Getting correct boundary submesh

I want to be able to define various boundary SubMesh's but am having a hard time
determining the correct inputs to MeshFunction in order that I can obtain a SubMesh with an expected number of vertices.

Given the discretization and subdomains defined by Top and Bottom,
I expect to obtain SubMesh's, submesh_bottom, submesh_top with
100 vertices each. Here is the output

It is true that when I use markers = dl.MeshFunction('size_t', bmesh, 2) the coordinates of
the SubMesh's are as expected. Point for clarification: can someone point out why I can’t use markers = dl.MeshFunction('size_t', mesh, 2)?

The second issue is that I can not use these subsets to define boundary conditions over the original mesh.

Question: how can I extract a correct SubMesh and also be able to use the MeshFunction data to define a DirichletBC as well as define a Measure to include a surface integral term over the correct sub boundary domain in the weak form that defines my PDE?

Hello,

As far as I know SubMesh is intended to create submeshes composed of a subset of cells from the parent mesh (i.e. the initial mesh and the submesh have the same topological dimension). That is probably why it doesn’t work as expected when to try to build a submesh composed of a set of lower dimensional entities (and why it works with bmesh, top_dim = 2).
Note : you can export the markers and your meshes to XDMF files to visualize them (with Paraview for instance) and check what you get :

with XDMFFile(mesh.mpi_comm(), filename) as file:
    file.write(boundary_markers)
with XDMFFile(mesh.mpi_comm(), filename) as file:
    file.write(submesh_bottom)

In the context of mixed-dimensional problems, we have developed a new MeshView class allowing to build submeshes of lower dimension. It should soon be merged into the master branch but meanwhile, you can give it a try using the branch cecile/mixed-dimensional or the dedicated Docker container ceciledc/fenics_mixed_dimensional. The syntax would be :

submesh_bottom = MeshView.create(boundary_markers, 1)

P.S. : For next posts, please provide the code directly in the post and not as an image. It helps people who want to reproduce your issue

2 Likes

Fantastic, thank you for the clarification!