You can follow the logic in: Wrong FacetNormal vector on internal boundaries - #2 by dokken
i.e. you pack your integration data by hand, and choose consistently that a cell marked with a value (say a will be the first in the integration entities array).
Please note that is probably better to use a dg-0 function to indicate right and left value.
Here is a reference implementation that does this handling consistently:
# Consistent interior integral
# Author: Jørgen S. Dokken
# SPDX-License-Identifier: MIT
from mpi4py import MPI
import dolfinx
import numpy as np
import packaging.version
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10, dolfinx.cpp.mesh.CellType.triangle, ghost_mode=dolfinx.mesh.GhostMode.shared_facet)
# Divide mesh in two parts
def left(x):
return x[0] < 0.5 + 1e-14
num_cells = mesh.topology.index_map(mesh.topology.dim).size_local+mesh.topology.index_map(mesh.topology.dim).num_ghosts
cells = np.arange(num_cells, dtype=np.int32)
values = np.full_like(cells, 3, dtype=np.int32) # Set all cells to 3
values[dolfinx.mesh.locate_entities(mesh, mesh.topology.dim, left)] = 5 # Set left side to 5
# Create meshtag
ct = dolfinx.mesh.meshtags(mesh, mesh.topology.dim, cells, values)
# Find common facets
mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim - 1)
cell_left = ct.find(5)
cell_right = ct.find(3)
left_facets= dolfinx.mesh.compute_incident_entities(mesh.topology, cell_left, mesh.topology.dim, mesh.topology.dim-1)
right_facets= dolfinx.mesh.compute_incident_entities(mesh.topology, cell_right, mesh.topology.dim, mesh.topology.dim-1)
common_facets = np.intersect1d(left_facets, right_facets)
if packaging.version.parse(dolfinx.__version__) < packaging.version.parse("0.9"):
integration_entities = []
c_to_f = mesh.topology.connectivity(mesh.topology.dim, mesh.topology.dim - 1)
f_to_c = mesh.topology.connectivity(mesh.topology.dim-1, mesh.topology.dim)
for f in common_facets:
cells = f_to_c.links(f)
if len(cells) == 2:
cell_values = ct.values[cells]
if len(np.unique(cell_values)) < 2:
# Only integrate over common boundary
continue
# We set all highest values first ("+") restriction
argsort = np.argsort(cell_values)
for cell in cells[argsort[::-1]]:
facets = c_to_f.links(cell)
pos = np.argwhere(facets == f)
integration_entities.append(cell)
integration_entities.append(pos[0, 0])
else:
pass
else:
mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim - 1)
integration_entities = dolfinx.fem.compute_integration_domains(dolfinx.fem.IntegralType.interior_facet, mesh.topology, common_facets,
mesh.topology.dim-1)
connected_cells = integration_entities.reshape(-1, 4)[:, [0, 2]]
switch = ct.values[connected_cells[:,0]] < ct.values[connected_cells[:,1]]
if len(switch)> 0 and any(switch):
tmp_entities = integration_entities.reshape(-1, 4)
tmp_entities[switch] = tmp_entities[switch][:, [2, 3, 0, 1]]
integration_entities = tmp_entities.flatten()
import ufl
dS_custom = ufl.Measure("dS", domain=mesh, subdomain_data=[(12, np.array(integration_entities, dtype=np.int32))])
# Reference implementation, which doesn't work
# facet_marker = dolfinx.mesh.meshtags(mesh, mesh.topology.dim-1, common_facets, np.full_like(common_facets, 12))
# dS_custom = ufl.Measure("dS", subdomain_data=facet_marker)
V = dolfinx.fem.functionspace(mesh, ("DG", 0, (2, )))
u = dolfinx.fem.Function(V)
# Set all values on right side to (2.5, 2.5)
if packaging.version.parse(dolfinx.__version__) < packaging.version.parse("0.9"):
u.interpolate(lambda x: np.full((2, x.shape[1]), 2.5, dtype=np.float64), cells=ct.find(3))
else:
u.interpolate(lambda x: np.full((2, x.shape[1]), 2.5, dtype=np.float64), cells0=ct.find(3))
u.x.scatter_forward()
# Set all values on left side to (0.3, 0)
def g(x):
values = np.zeros((2, x.shape[1]), dtype=np.float64)
values[0] = 0.3
return values
if packaging.version.parse(dolfinx.__version__) < packaging.version.parse("0.9"):
u.interpolate(g, cells=ct.find(5))
else:
u.interpolate(g, cells0=ct.find(5))
u.x.scatter_forward()
# Evaluate integral from either side
n = ufl.FacetNormal(mesh)
f_left = dolfinx.fem.form(ufl.dot(u("+"),n("+")) * dS_custom(12))
left_val = mesh.comm.allreduce(dolfinx.fem.assemble_scalar(f_left), op=MPI.SUM)
f_right = dolfinx.fem.form(ufl.dot(u("-"),n("-")) * dS_custom(12))
right_val = mesh.comm.allreduce(dolfinx.fem.assemble_scalar(f_right), op=MPI.SUM)
print(left_val, right_val)
assert np.isclose(left_val, 0.3)
assert np.isclose(right_val, -2.5)