Identify facets on the (outer) boundary of a 3D mesh

mesh
,
1

I have a complicated 3D mesh (loaded from a xdmf-file with cell data available, see the screenshot below for a clipped representation).

To set the boundary conditions, I’m wondering if there is there a direct way (or a workaround) to automatically identify the facets or dofs of the mesh that sit at outer boundaries?

(I cannot use, say, the dolfinx.fem.locate_dofs_geometrical function because I don’t have a parametrization of the boundaries)

clipped view of the geometry

2

Would something like dolfinx/python/demo/demo_biharmonic.py at 5f2b765005b8bd28b138a5959b9a46ad76e277d8 · FEniCS/dolfinx · GitHub work for you?

1 Like
3

Worked for me – after I set entity_dim=2 in the call to fem.locate_dofs_topological in line 174

tdim = mesh.topology.dim
mesh.topology.create_connectivity(tdim - 1, tdim)
facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
os_dofs = locate_dofs_topological(V=V, entity_dim=2, entities=facets)

# ## write out the facet field
fctfnm = 'facet_function.bp'
fctpltfun = dolfinx.fem.Function(V)
fctpltfun.name = 'facets_idf'
fct_cnv = dolfinx.io.VTXWriter(V.mesh.comm, fctfnm, fctpltfun, engine="BP4")
fctpltfun.x.array[os_dofs] = 1.  # set function to 1 at the surface dofs
fct_cnv.write(-1)  # write the function

Only now I have the problem that also interior surfaces (between the different subregions) are detected. Any quick idea on that? Otherwise, I will open another topic.

facet function (including interior facets where subregions meet)

4

That’s odd. If you look at the implementation in dolfinx/cpp/dolfinx/mesh/utils.cpp at 5f2b765005b8bd28b138a5959b9a46ad76e277d8 · FEniCS/dolfinx · GitHub exterior_facet_indices is marking only faces on the outer boundary (i.e., with only one connected cell), not between subdomains. Can you double check if somehow the mesh you have is “detached” at the interfaces between subdomains?

5

Thank you for this insight. In fact, when I filtered out the dofs that belonged to a subregion, it appeared that “the other half of the facet” remains (see the second image below). :face_with_monocle:

both facets
one facet

Looks like I have to talk back to the engineers. The more that I expect that this will be problematic when I’m going to solve the PDEs on the mesh.