Problem transferring facet tags to submesh

mesh
1

Hello,

I have generated a 3D mesh using Gmsh, which includes several physical volumes (103, 104, 105, 106, 107) and several physical surfaces (1, 2) to impose boundary conditions.

Next, I would like to solve my variational form only on a part of the mesh (on volumes 106 and 107), while imposing boundary conditions with the physical surfaces from the main mesh (1, 2), which are also surfaces of my volumes 106 and 107.

To achieve this, I created a submesh from the main mesh and then tried to transfer the facet tags from the parent mesh to the submesh, following the procedure provided in the following post: Transfer meshtags to submesh in dolfinx - #7 by dokken

However, I encountered an error. Could you help me identify the problem?

Thank you very much!

The full code is as follows:

import dolfinx
from dolfinx.io import XDMFFile
from mpi4py import MPI
import numpy as np

# Read the mesh
with XDMFFile(MPI.COMM_WORLD, "3Dmesh/3Dblock_additionalMaterials_mesh.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh(name="Grid")
    cell_tags = xdmf.read_meshtags(mesh, name="Grid")
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)

with XDMFFile(MPI.COMM_WORLD, "3Dmesh/3Dblock_additionalMaterials_facet_mesh.xdmf", "r") as xdmf:
    facet_tags = xdmf.read_meshtags(mesh, name="Grid")

# Check physical volumes and surfaces of the mesh
print(np.unique(cell_tags.values))
print(np.unique(facet_tags.values))

# Create submesh which is the union of physical volumes 106 and 107
submesh, entity_map, _ = dolfinx.mesh.create_submesh(mesh, mesh.topology.dim, cell_tags.indices[(cell_tags.values==106)|(cell_tags.values==107)])

# Transfer facet tags from parent mesh to submesh 
tdim = mesh.topology.dim
fdim = tdim - 1
c_to_f = mesh.topology.connectivity(tdim, fdim)
f_map = mesh.topology.index_map(fdim)
all_facets = f_map.size_local + f_map.num_ghosts
all_values = np.zeros(all_facets, dtype=np.int32)
all_values[facet_tags.indices] = facet_tags.values
print(np.unique(all_values))

submesh.topology.create_entities(fdim)
subf_map = submesh.topology.index_map(fdim)
submesh.topology.create_connectivity(tdim, fdim)
c_to_f_sub = submesh.topology.connectivity(tdim, fdim)
num_sub_facets = subf_map.size_local + subf_map.size_global
sub_values = np.empty(num_sub_facets, dtype=np.int32)
for i, entity in enumerate(entity_map):
    parent_facets = c_to_f.links(entity)
    child_facets = c_to_f_sub.links(i)
    for child, parent in zip(child_facets, parent_facets):
        sub_values[child] = all_values[parent]
sub_meshtag = dolfinx.mesh.meshtags(submesh, submesh.topology.dim-1, np.arange(
    num_sub_facets, dtype=np.int32), sub_values)

The error is:

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Input In [1], in <cell line: 38>()
     40     child_facets = c_to_f_sub.links(i)
     41     for child, parent in zip(child_facets, parent_facets):
---> 42         sub_values[child] = all_values[parent]
     43 sub_meshtag = dolfinx.mesh.meshtags(submesh, submesh.topology.dim-1, np.arange(
     44     num_sub_facets, dtype=np.int32), sub_values)

IndexError: index 128753 is out of bounds for axis 0 with size 64368

The mesh can be accessed at the following link:

2

I cannot reproduce this with DOLFINx v0.6.0:

import dolfinx
from dolfinx.io.gmshio import read_from_msh
from mpi4py import MPI
import numpy as np

mesh, cell_tags, facet_tags = dolfinx.io.gmshio.read_from_msh("mesh.msh", MPI.COMM_WORLD, 0)

# Check physical volumes and surfaces of the mesh
print(np.unique(cell_tags.values))
print(np.unique(facet_tags.values))

# Create submesh which is the union of physical volumes 106 and 107
submesh, entity_map, vertex_map, geom_map = dolfinx.mesh.create_submesh(mesh, mesh.topology.dim, cell_tags.indices[(cell_tags.values==106)|(cell_tags.values==107)])

# Transfer facet tags from parent mesh to submesh 
tdim = mesh.topology.dim
fdim = tdim - 1
c_to_f = mesh.topology.connectivity(tdim, fdim)
f_map = mesh.topology.index_map(fdim)
all_facets = f_map.size_local + f_map.num_ghosts
all_values = np.zeros(all_facets, dtype=np.int32)
all_values[facet_tags.indices] = facet_tags.values
print(np.unique(all_values))

submesh.topology.create_entities(fdim)
subf_map = submesh.topology.index_map(fdim)
submesh.topology.create_connectivity(tdim, fdim)
c_to_f_sub = submesh.topology.connectivity(tdim, fdim)
num_sub_facets = subf_map.size_local + subf_map.num_ghosts
sub_values = np.empty(num_sub_facets, dtype=np.int32)
for i, entity in enumerate(entity_map):
    parent_facets = c_to_f.links(entity)
    child_facets = c_to_f_sub.links(i)
    for child, parent in zip(child_facets, parent_facets):
        sub_values[child] = all_values[parent]
sub_meshtag = dolfinx.mesh.meshtags(submesh, submesh.topology.dim-1, np.arange(
    num_sub_facets, dtype=np.int32), sub_values)

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "submesh.xdmf", "w") as xdmf:
    xdmf.write_mesh(submesh)
    submesh.topology.create_connectivity(submesh.topology.dim-1, submesh.topology.dim)
    xdmf.write_meshtags(sub_meshtag)
2 Likes
3

Thank you for your response, I was using the previous version of dolfinx, hence the error.