MeshView submesh vertex map to parent mesh

1

Hello,

I am getting started in using the cecile/mixed-dimensional branch. I have a 3D cube geometry with a cylindric hole at the center. I have defined the inner surface of the cylinder as a 2D subdomain using

> class Inner(SubDomain):
>     def inside(self,x,on_boundary):
>         tol = 1e-1
>         return (on_boundary and np.linalg.norm(x[0:2]) <= ri + tol)
>
> marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
> 
> marker.set_all(0)
> sm = Inner()
> sm.mark(marker, 1)
> marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
> submesh = MeshView.create(marker, 1)

My question is how do I map the vertex of the submesh to the parent mesh and vice-versa? Is there an easy way of doing this? I tried with some functions of MeshView but I can’t make it work.

Thank you!

1 Like
2

Hello,

The mapping between parent and child mesh vertices (and cells) is built when building the MeshView object. You can access these from the topology of your submesh :

submesh.topology().mapping()[mesh.id()].vertex_map()
submesh.topology().mapping()[mesh.id()].cell_map()

Note : In your code you redefine the variable marker after having marked your entities, making it zero for all entities so I suspect your submesh to be empty.

2 Likes
3

Thank you for the answer!

4

Hello, does there exist a similar mapping between the parent and child mesh facets ?

5

Hello,
There is no such mapping for mesh facets. Note that if your submesh is of codimension 1, then the submesh cells are facets in the parent mesh that you can access through the corresponding cell_map()

2 Likes
6

Thank you very much for your answer !

In fact, I was asking this question in relation with this post. Such a mapping would have been useful to recover facet labels on a submesh (of codimension 0) after using MeshView, but I will have to find another way then.

8

Hi @cdaversin. Nice approach. My question is that is it also possible to get the vertex_map() and cell_map() from the entity submesh if it was not created in same .py script rather was created by same MeshView but in some other .py code and wrote in .h5 file and then read it in separate code file by reading submesh from that .h5 file?
I am currently working with the scripts where mesh generation using MeshView is done in some dedicated code file and they are read to other scripts using .h5 files. Thanks!

9

Hi,

When creating a submesh with MeshView, the new Mesh object which is generated embeds a map of the mappings it has with other meshes (parent or sibling meshes), as indicated in my previous post :

submesh.topology().mapping()

Unfortunately this info is not stored when writing your submesh as a .h5 (or .xdmf) file. Here is the code sample I wrote to check this :

from fenics import *

# Create mesh and its submesh
mesh = UnitSquareMesh(10,10)
marker = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
for f in facets(mesh):
   marker[f] = 0.5 - DOLFIN_EPS < f.midpoint().y() < 0.5 + DOLFIN_EPS
submesh = MeshView.create(marker, 1)

# Write meshes to .h5
mesh_hfile_ = "XDMF/mesh.h5"
mesh_hfile = HDF5File(MPI.comm_world, mesh_hfile_, "w")
mesh_hfile.write(mesh, "/mesh")
mesh_hfile.close()

submesh_hfile_ = "XDMF/submesh.h5"
submesh_hfile = HDF5File(MPI.comm_world, submesh_hfile_, "w")
submesh_hfile.write(submesh, "/mesh")
submesh_hfile.close()

# Load meshes in new Mesh objects
loaded_mesh = Mesh()
mesh_hfile = HDF5File(MPI.comm_world, mesh_hfile_, "r")
mesh_hfile.read(loaded_mesh, "/mesh", False)
mesh_hfile.close()

loaded_submesh = Mesh()
submesh_hfile = HDF5File(MPI.comm_world, submesh_hfile_, "r")
submesh_hfile.read(loaded_submesh, "/mesh", False)
submesh_hfile.close()

# Does loaded submesh have a mapping ?
mapping = submesh.topology().mapping()
print("mapping (initial submesh) = ", mapping)
loaded_mapping = loaded_submesh.topology().mapping()
print("mapping (loaded submesh) = ", loaded_mapping)

Which returns :

mapping (initial submesh) =  {0: <dolfin.cpp.mesh.MeshView object at 0x7efdcb7b2f70>}
mapping (loaded submesh) =  {}

The submesh generated in this script is mapped with its parent of index 0, while the same submesh loaded from a .h5 file doesn’t have any mapping. So I am afraid you will have to do that in the same script, or find a way to store the mappings in a separate file.

1 Like