Plotting subdomains/internal boundaries

mesh
8

There seems to be a bug when using subdomains to mark cells in parallel with shared_vertex/shared_facet.
You can submit a bug at the bitbucket page. However, it is not likely to be fixed as DOLFIN is no longer developed.

A workaround would be to write the tags to XDMF in serial, and then read then in again when running in parallel.

In the development version of dolfin, dolfinx this works out of the box:

import numpy as np
import dolfinx
from mpi4py import MPI

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 25, 25,
                                       ghost_mode=dolfinx.mesh.GhostMode.shared_facet)


def Omega(x):
    return (x[0] - 0.5)**2 + (x[1]-0.5)**2 - 0.2**2 > 0


tdim = mesh.topology.dim
indices = dolfinx.mesh.locate_entities(mesh, tdim, Omega)

num_cells = mesh.topology.index_map(tdim).size_local + \
    mesh.topology.index_map(tdim).num_ghosts
local_cell_marker = np.zeros(num_cells, dtype=np.int32)
local_cell_marker[indices] = 1
ct = dolfinx.mesh.MeshTags(mesh, tdim, np.arange(
    num_cells, dtype=np.int32), local_cell_marker)


mesh.topology.create_connectivity(tdim-1, tdim)
f_to_c = mesh.topology.connectivity(tdim-1, tdim)


num_facets = mesh.topology.index_map(tdim-1).size_local + \
    mesh.topology.index_map(tdim-1).num_ghosts

facets = []

for i in range(num_facets):
    cells = f_to_c.links(i)
    if len(cells) == 2:
        values = np.unique(local_cell_marker[cells])
        if len(values) == 2:
            facets.append(i)

all_facets = np.zeros(num_facets, dtype=np.int32)
all_facets[facets] = 1
ft = dolfinx.mesh.MeshTags(mesh, tdim-1, np.arange(num_facets, dtype=np.int32),
                           all_facets)

with dolfinx.io.XDMFFile(mesh.comm, "mf.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    ct.name = "Cell tag"
    xdmf.write_meshtags(ct)
    ft.name = "Facet tag"
    xdmf.write_meshtags(ft)


print(MPI.COMM_WORLD.rank, facets)
1 Like