Extracting 2D boundary mesh from 3D mesh

mesh
1

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.

2

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.

3

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)


in_mesh = meshio.read("bot_mesh.xdmf")

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:
    xdmf.read(mesh2D)
print(mesh2D.geometry().dim())
2 Likes
4

Thank you, that seems to work!

5

Hi,
I am able to extract boundary mesh using above MWE, but, unable to obtain the subdomains from original 3D mesh. The application is cylinder inscribed in a cube. (sub-domain1=cylinder, subdomain2=cube). Like using mesh0 = SubMesh(mesh, subdomains, 0) gives cylinder mesh.

How can we get the subdomain mesh at boundary, bot_boundary which would be a 2D circle?

Like, the bot_boundary=SubMesh(bmesh, part_of_bot, 1) giving submesh at y=-1 boundary.
, can we extract original subdomain data from bot_boundary.

I tried using,
bot_boundary_0=SubMesh(bot_boundary, subdomains,0) and creating bmesh using only mesh0, but, all causing entire page shut. Any help is greatly appreciated…