method `global_index()` crashes

Question

The code crashes if the method global_index() is called to face entity, see the marked line.

How to get the global_index of a mesh entity? I need this the create face to face mapping between the sub mesh and the parent mesh. Any ideas how to do this?

import dolfin
import numpy as np
dolfin.set_log_level(dolfin.ERROR)
class LeftHalf(dolfin.SubDomain):
    def inside(self, x, on_boundary):
        return dolfin.between(x[0], (0.0, 0.5))
class LeftSide(dolfin.SubDomain):
    def inside(self, x, on_boundary):
        return ( dolfin.near(x[0], 0.0) and on_boundary )
fem_mesh = dolfin.BoxMesh(dolfin.Point(0.,0.,0.), dolfin.Point(1.,1.,1.),10,10,10)
fem_spacedim = fem_mesh.geometry().dim()

left_domain = LeftHalf()
left_boundary = LeftSide()

cell_ids = dolfin.CellFunction("size_t", fem_mesh)
cell_ids.set_all(0)
left_domain.mark(cell_ids,1)


face_ids = dolfin.FaceFunction("size_t", fem_mesh)
face_ids.set_all(0)
left_boundary.mark(face_ids,1)

sub_mesh = dolfin.SubMesh(fem_mesh, cell_ids, 1)
# creating face mapping between meshes
parent_cells = list(dolfin.cells(fem_mesh))

child_faces = list(dolfin.faces(sub_mesh))
child_cells = list(dolfin.cells(sub_mesh))

n_child_faces = sub_mesh.num_faces()
assert n_child_faces > 0, "Number of faces must be larger zero."        

child_face_touched = np.zeros((n_child_faces, ), dtype = np.bool)

parent_face_indices = np.zeros((n_child_faces, ), dtype = np.uintp)

parent_cell_indices = sub_mesh.data().array("parent_cell_indices", fem_spacedim)
assert parent_cell_indices.size == sub_mesh.num_cells(), "Dimensions mismatch!"

parent_vertex_indices = sub_mesh.data().array("parent_vertex_indices", 0)
assert parent_vertex_indices.size == sub_mesh.num_vertices(), "Dimensions mismatch!"


child_face2parent_vertex = {}
parent_face2parent_vertex = {}

# loop over cells
for child_cell_ind, parent_cell_ind in enumerate(parent_cell_indices):
    parent_cell = parent_cells[parent_cell_ind]
    child_cell = child_cells[child_cell_ind]

    child_face2parent_vertex.clear()

    for child_face in dolfin.faces(child_cell):
        # crash in next line
        if not(child_face_touched[child_face.global_index()]): 
            child_face2parent_vertex[child_face.global_index()] = parent_vertex_indices[child_face.entities(0)]

    if len(child_face2parent_vertex) > 0:
        parent_face2parent_vertex.clear()

        for parent_face in dolfin.faces(parent_cell):
            parent_face2parent_vertex[parent_face.index()] = parent_face.entities(0)

        for sub_key in child_face2parent_vertex:
            for parent_key in parent_face2parent_vertex:
                if np.all(child_face2parent_vertex[sub_key] == parent_face2parent_vertex[parent_key]):
                    parent_face_indices[sub_key] = parent_key
                    child_face_touched[sub_key] = True
                    break
    assert np.all(child_face_touched), "Not all faces touched."

Comments

I was able to re-program the code using the method enitities() which returns the global indices. I ask myself what global_index is good for?

import dolfin
import numpy as np
dolfin.set_log_level(dolfin.ERROR)
class LeftHalf(dolfin.SubDomain):
    def inside(self, x, on_boundary):
        return dolfin.between(x[0], (0.0, 0.5))
class LeftSide(dolfin.SubDomain):
    def inside(self, x, on_boundary):
        return ( dolfin.near(x[0], 0.0) and on_boundary )
fem_mesh = dolfin.BoxMesh(dolfin.Point(0.,0.,0.), dolfin.Point(1.,1.,1.),10,10,10)
fem_spacedim = fem_mesh.geometry().dim()

left_domain = LeftHalf()
left_boundary = LeftSide()

cell_ids = dolfin.CellFunction("size_t", fem_mesh)
cell_ids.set_all(0)
left_domain.mark(cell_ids,1)

face_ids = dolfin.FaceFunction("size_t", fem_mesh)
face_ids.set_all(0)
left_boundary.mark(face_ids,1)

dolfin.plot(face_ids)
dolfin.interactive()

sub_mesh = dolfin.SubMesh(fem_mesh, cell_ids, 1)

parent_faces = list(dolfin.faces(fem_mesh))
parent_cells = list(dolfin.cells(fem_mesh))

child_faces = list(dolfin.faces(sub_mesh))
child_cells = list(dolfin.cells(sub_mesh))

n_child_faces = sub_mesh.num_faces()
assert n_child_faces > 0, "Number of faces must be larger zero."        

child_face_touched = np.zeros((n_child_faces, ), dtype = np.bool)

parent_face_indices = np.zeros((n_child_faces, ), dtype = np.uintp)

parent_cell_indices = sub_mesh.data().array("parent_cell_indices", fem_spacedim)
assert parent_cell_indices.size == sub_mesh.num_cells(), "Dimensions mismatch!"

parent_vertex_indices = sub_mesh.data().array("parent_vertex_indices", 0)
assert parent_vertex_indices.size == sub_mesh.num_vertices(), "Dimensions mismatch!"

child_face2parent_vertex = {}
parent_face2parent_vertex = {}

# loop over cells
for child_cell_ind, parent_cell_ind in enumerate(parent_cell_indices):
    parent_cell = parent_cells[parent_cell_ind]
    child_cell = child_cells[child_cell_ind]

    child_face2parent_vertex.clear()

    for child_face_ind in child_cell.entities(fem_spacedim - 1):
        if not(child_face_touched[child_face_ind]):
            child_face = child_faces[child_face_ind]
            child_face2parent_vertex[child_face_ind] = np.sort(parent_vertex_indices[child_face.entities(0)])
    if len(child_face2parent_vertex) > 0:
        parent_face2parent_vertex.clear()

        for parent_face_ind in parent_cell.entities(fem_spacedim - 1):
            parent_face = parent_faces[parent_face_ind]
            parent_face2parent_vertex[parent_face_ind] = np.sort(parent_face.entities(0))

        for sub_key in child_face2parent_vertex:
            for parent_key in parent_face2parent_vertex:
                if np.all(child_face2parent_vertex[sub_key] == parent_face2parent_vertex[parent_key]):
                    parent_face_indices[sub_key] = parent_key
                    child_face_touched[sub_key] = True
                    break
assert np.all(child_face_touched), "Not all faces touched."

sub_face_ids = dolfin.FaceFunction("size_t", sub_mesh)
sub_face_ids.set_all(0)
sub_face_ids.array()[:] = face_ids.array()[parent_face_indices]

dolfin.plot(sub_face_ids)
dolfin.interactive()