Tag entities between two subdomains

mesh
1

Hi everybody,

I am looking for a way to tag the interface between two subdomains in dolfinx (if this is not already done in, e.g., gmsh).
The cell tags are defined for each cell. My idea was to write all nodes that belong to each cell tag to a list and find duplicate nodes as they belong to the interface. However, I do not know how to proceed. Does anyone have any idea?
Thank you.

2

The solution is similar to this post: Plotting subdomains/internal boundaries - #8 by dokken

3

You can use the facet_to_cell connectivity to do this.

Consider the following:

import numba
import numpy.typing as npt
import dolfinx
from mpi4py import MPI
import numpy as np
import time

N = 100
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, N, N)


def some_cells(x):
    return (abs(x[0]-0.5) < 0.2) & (abs(x[1]-0.5) < 0.4)


@numba.njit(cache=True)
def create_interface(facet_marker: npt.NDArray[np.int32],
                     tag: int,
                     fa: npt.NDArray[np.int32],
                     fo: npt.NDArray[np.int32],
                     cell_marker: npt.NDArray[np.int32]):
    for f in range(len(facet_marker)):
        cells = fa[fo[f]:fo[f+1]]
        c_val = cell_marker[cells]
        if len(np.unique(c_val)) == 2:
            facet_marker[f] = tag


# Cold start of function
create_interface(np.empty(0, dtype=np.int32), 0, np.empty(
    0, dtype=np.int32), np.empty(0, dtype=np.int32), np.empty(0, dtype=np.int32))

cell_tag_0 = 1
cell_tag_1 = 3
cells = dolfinx.mesh.locate_entities(mesh, mesh.topology.dim, some_cells)
ct = dolfinx.mesh.meshtags(mesh, mesh.topology.dim,
                           cells, np.full_like(cells, cell_tag_0, dtype=np.int32))

mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
f_to_c = mesh.topology.connectivity(mesh.topology.dim-1, mesh.topology.dim)
fa = f_to_c.array
fo = f_to_c.offsets

facet_map = mesh.topology.index_map(mesh.topology.dim-1)
num_facets = facet_map.size_local + facet_map.num_ghosts
facet_marker = np.zeros(num_facets, dtype=np.int32)
cell_map = mesh.topology.index_map(mesh.topology.dim)
num_cells = cell_map.size_local + cell_map.num_ghosts
cell_marker = np.full(num_cells, cell_tag_1, dtype=np.int32)
cell_marker[ct.indices] = ct.values

interface_marker = 7
start = time.perf_counter()
create_interface(facet_marker, interface_marker, fa, fo, cell_marker)
end = time.perf_counter()
print(f"Find interface {end-start:.2e}")
ft = dolfinx.mesh.meshtags(mesh, mesh.topology.dim-1,
                           np.arange(num_facets), facet_marker)
with dolfinx.io.XDMFFile(mesh.comm, "facet_marker.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(ft)
1 Like
4

Okay great! I will also consider this new solution. Thank you!