 Extracting 2D boundary mesh from 3D mesh

Hello all,

I’m currently trying to extract a certain 2D portion of the boundary of a 3D mesh. As a minimum working example, consider trying to extract the bottom surface of a unit cube.

from dolfin import *

dim = 3
N = 2**5
mesh = UnitCubeMesh(N,N,N)

bdim = dim-1
bmesh = BoundaryMesh(mesh, "exterior")
mapping = bmesh.entity_map(bdim)
part_of_bot = MeshFunction("size_t", bmesh, bdim)
for cell in cells(bmesh):
curr_facet_normal = Facet(mesh, mapping[cell.index()]).normal()
if near(curr_facet_normal.y(), -1.0): #On bot boundary
part_of_bot[cell] = 1

bot_boundary = SubMesh(bmesh, part_of_bot, 1)
File('bot_boundary.pvd') << bot_boundary

When directly examining the output mesh files, I can see that it preserved the y coordinate (which makes sense). However, I need the 2D mesh, with the zero y coordinate pruned. How might I go about getting rid of that coordinate from bot_boundary, or better yet, is there a way to build bot_boundary (say an option that I can pass to SubMesh) that doesn’t preserve the y coordinate in the first place? I can’t find a relevant function in the docs (Restriction seemed relevant, but that’s been removed since).

Thanks for any input.

If it helps, my goal is eventually to pass the resulting 2d mesh to VectorFunctionSpace. There’s a dimension option for VectorFunctionSpace which I believe will let it ignore the 3rd coordinate, however that won’t work my application and the MWE as they lie in the X-Z plane. I suppose I could rotate bot_boundary into the X-Y plane and proceed like that.

There must be a better way to do it than that though.

Consider the following MWE using meshio (which can be installed with pip), for instance
pip3 install --no-binary=h5py h5py meshio --user

import meshio
import numpy as np
from dolfin import *

dim = 3
N = 2**5
mesh = UnitCubeMesh(N, N, N)

bdim = dim-1
bmesh = BoundaryMesh(mesh, "exterior")
mapping = bmesh.entity_map(bdim)
part_of_bot = MeshFunction("size_t", bmesh, bdim)
for cell in cells(bmesh):
curr_facet_normal = Facet(mesh, mapping[cell.index()]).normal()
if near(curr_facet_normal.y(), -1.0):  # On bot boundary
part_of_bot[cell] = 1

bot_boundary = SubMesh(bmesh, part_of_bot, 1)
#File('bot_boundary.pvd') << bot_boundary
with XDMFFile("bot_mesh.xdmf") as xdmf:
xdmf.write(bot_boundary)

cells = in_mesh.get_cells_type("triangle")
points = np.delete(in_mesh.points, 1, axis=1)
out_mesh = meshio.Mesh(points=points, cells={"triangle": cells})
meshio.write("pruned_mesh.xdmf", out_mesh)

mesh2D = Mesh()
with XDMFFile("pruned_mesh.xdmf") as xdmf: