Create corresponding cell tags from existing facet tags

dolfinx
1

Hello all,

I was wondering if there is a preferred way in dolfinx to create the corresponding cell tags from existing facet tags. So basically mark every element that belongs to a tagged facet with the same numerical tag as a cell tag.

Version:
fenics-dolfinx 0.7.2
gmsh 4.12.2

MWE:

from dolfinx import mesh, io
from dolfinx.mesh import compute_midpoints, locate_entities
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.fem import Function, FunctionSpace
from ufl import FiniteElement, Mesh
from dolfinx.io import XDMFFile

import gmsh

gmsh.initialize()
gmsh.model.add("box_with_tags")

# Create box (x, y, z, Lx, Ly, Lz)
box = gmsh.model.occ.addBox(0, 0, 0, 10, 10, 10)
gmsh.model.occ.synchronize()

# tags for top and bottom
bottom_tag = 1
top_tag = 2

# Create PhysicalGroups
gmsh.model.addPhysicalGroup(3, [1], tag=3)
gmsh.model.addPhysicalGroup(2, [1], bottom_tag)  
gmsh.model.addPhysicalGroup(2, [2], top_tag)  

gmsh.model.occ.synchronize()
gmsh.model.mesh.generate(3)

gmsh.write("box_with_tags.msh")

gmsh.finalize()

def create_celltags(msh, facetags, tag1, tag2):

    #Create connectivity
    msh.topology.create_connectivity(msh.topology.dim, msh.topology.dim - 1)
    facet_to_cell_conn = msh.topology.connectivity(msh.topology.dim, msh.topology.dim - 1)
    
    #Get facets with tag1 or tag2
    facet_indices_tag1 = np.where(facetags.values == tag1)[0]
    facet_indices_tag2 = np.where(facetags.values == tag2)[0]

    #Find the related element_indices
    elements_tag1 = np.unique(np.hstack([facet_to_cell_conn.links(f) for f in facet_indices_tag1]))
    elements_tag2 = np.unique(np.hstack([facet_to_cell_conn.links(f) for f in facet_indices_tag2]))

    #Combine tag1 and tag2
    adjacent_elements = np.concatenate((elements_tag1, elements_tag2))

    
    cell_marker_values = np.concatenate((
        np.full(len(elements_tag1), tag1, dtype=np.int32),  # Tag 1 → 1
        np.full(len(elements_tag2), tag2, dtype=np.int32)   # Tag 2 → 2
    ))

    #create cell_tags
    cell_tags = mesh.meshtags(domain, domain.topology.dim, adjacent_elements, cell_marker_values)
    
    print(f"Number of marked cells: {len(adjacent_elements)}")
    return cell_tags




domain, markers, facets = io.gmshio.read_from_msh("box_with_tags.msh", MPI.COMM_WORLD)

cell_tags = create_celltags(domain, facets, tag1=1, tag2=2)
cell_tags.name = "Cell tags"  

with XDMFFile(MPI.COMM_WORLD, "box_with_tags.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_meshtags(facets, domain.geometry)
    xdmf.write_meshtags(cell_tags, domain.geometry)

I dont get an error but i see in the visualization, that the wrong cells are marked.
Edit: I cant post a picture of the tagged facets, since new users can only upload one picture, but its top and bottom.

image
image831×654 47.9 KB