Confusion on mapping BoundaryMesh to parent mesh

mesh
1

Hello, I am trying to figure out how to relate face indices on a BoundaryMesh with the relevant indices on the parent mesh. In the minimal code below I attempt to extract the faces on the bottom of a cube using MeshFunction but I must not be understanding how BoundaryMesh.entity_map() works. Any help as to why the first example doesn’t work but the second one does would be very much appreciated… Thank you!

from dolfin import *

# volumetric mesh and its boundary mesh
vmesh = UnitCubeMesh(10, 10, 10)
bmesh = BoundaryMesh(vmesh, "exterior")
# meshfunctions of dimension 2 for the volumetric/boundary mesh
mfv = MeshFunction("size_t", vmesh, 2)
mfb = MeshFunction("size_t", vmesh, 2)

# maps face indices on the boundary mesh to indices on the full mesh
facemap = bmesh.entity_map(2)

# expect smv to extract all faces on the boundary pointing downward from vmesh
for f in faces(bmesh):
    if Facet(vmesh, cellmap[f.index()]).normal().z() < -0.5:
        mfv[facemap[f.index()]] = 1

smv = SubMesh(vmesh,mfv,1)

# meanwhile somehow this code works
# expect mf to have value "1" for each face on boundary mesh pointing downward
for f in faces(bmesh):
    if Facet(vmesh, cellmap[f.index()]).normal().z() < -0.5:
        mfb[f] = 1

smb = SubMesh(bmesh,mfb,1)


vtkfile("smv.pvd") << smv
vtkfile("smb.pvd") << smb

SubMesh from mfv:

mfv
mfv.png928×875 128 KB

[since I am a new user it won’t let me upload a screenshot of the SubMesh from mfb but the result is as expected (just the bottom surface of the cube)]

Edit: This is on the latest stable version of fenics using docker

2

Hi, consider

from dolfin import *
import numpy as np

# volumetric mesh and its boundary mesh
mesh = UnitCubeMesh(10, 10, 10)
bdry_mesh = BoundaryMesh(mesh, 'exterior')

# Map CELLS of bdry mesh to FACETS of mesh
c2f = bdry_mesh.entity_map(2)

# Select CELLS of bdry_mesh that are of interest
bottom_cells = np.fromiter((cell.index() for cell in cells(bdry_mesh) if near(cell.midpoint().z(), 0)),
                           dtype='uintp')

# Mark the corrsponding FACETS of mesh
facet_f = MeshFunction('size_t', mesh, mesh.topology().dim()-1, 0)
facet_f.array()[c2f.array()[bottom_cells]] = 1

File('foo.pvd') << facet_f