Consider the following minimal example:
import dolfinx
import numpy as np
from mpi4py import MPI
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
V = dolfinx.fem.FunctionSpace(mesh, ("DG", 1))
u = dolfinx.fem.Function(V)
mesh.topology.create_entities(mesh.topology.dim -1)
f_map = mesh.topology.index_map(mesh.topology.dim -1)
num_facets = f_map.size_local + f_map.num_ghosts
facets = np.arange(num_facets, dtype=np.int32)
values = np.zeros_like(facets, dtype=np.int32)
def interface(x):
return np.isclose(x[0], 0.5)
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
values[dolfinx.mesh.locate_entities(mesh, mesh.topology.dim-1, interface)] = 1
ft = dolfinx.mesh.meshtags(mesh, mesh.topology.dim-1, facets, values)
c_map = mesh.topology.index_map(mesh.topology.dim)
num_cells = c_map.size_local + c_map.num_ghosts
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim-1)
c_to_f = mesh.topology.connectivity(mesh.topology.dim, mesh.topology.dim-1)
for cell in range(num_cells):
for facet in c_to_f.links(cell):
if ft.values[facet] == 1:
print(cell, facet, V.dofmap.cell_dofs(cell))
u.x.array[V.dofmap.cell_dofs(cell)] = 2
with dolfinx.io.VTXWriter(mesh.comm, "u.bp", [u], engine="BP4") as bp:
bp.write(0.0)
giving:
