Fix consistent restriction of one-sided interior facet integrals

[jorgensd]

Currently on main, there is no way of doing consistent one-sided integration on interior facets.
In old FEniCS, the workaround for this was to supply a cell-marker with values such that "+" would always be the side with the largest marker value. This PR reintroduces this fix, as it is crucial for post-processing in many applications, such as torque calculations in electromagnetics.

MWE:

from petsc4py import PETSc
import ufl
import dolfinx
from mpi4py import MPI
import numpy as np

mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, 10, 10, 10)
tdim = mesh.topology.dim
fdim = tdim - 1

plus_restriction = 5
minus_restriction = 3
facet_marker = 1
minus_cells = dolfinx.mesh.locate_entities(mesh, tdim, lambda x: x[1] <= 0.5 + 1e-16)
cell_map = mesh.topology.index_map(tdim)
all_cells = np.arange(cell_map.size_local + cell_map.num_ghosts, dtype=np.int32)
# All cells with x[1]>0.5 is tagged with plus restriction
cell_marker = np.full(len(all_cells), plus_restriction, dtype=np.int32)
# All cells with x[1]<=0.5 tagged with minus restriction
cell_marker[minus_cells] = 1
ct = dolfinx.mesh.meshtags(mesh, tdim, all_cells, cell_marker)

# Tag all facets on surface between tags with facet marker
midline_facets = dolfinx.mesh.locate_entities(mesh, fdim, lambda x: np.isclose(x[1], 0.5))
argsort_f = np.argsort(midline_facets)
ft = dolfinx.mesh.meshtags(mesh, fdim, midline_facets[argsort_f],
                           np.full(len(midline_facets), facet_marker, dtype=np.int32))


# Interpolate non-zero vlaue on minus restriction
V = dolfinx.fem.FunctionSpace(mesh, ("DG", 0))
u = dolfinx.fem.Function(V)
u.x.array[:] = 0
u.interpolate(lambda x: np.ones(x.shape[1]), minus_cells)
u.x.scatter_forward()

with dolfinx.io.XDMFFile(mesh.comm, "func.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_function(u, 0)

dS = ufl.Measure("dS", domain=mesh, subdomain_data=ft)
dx = ufl.Measure("dx", domain=mesh, subdomain_data=ct)
# Arbitrary orientation of +
form_0 = u("+") * dS(1)
# Arbitrary orientation of -
form_1 = u("-") * dS(1)
zero = dolfinx.fem.Constant(mesh, PETSc.ScalarType(0))
# + should be on side with larger value, i.e. exact value should be 0
form_2 = u("+") * dS(1) + zero * dx(88)
# - should point out of domain with smaller value, exact value is 1
form_3 = u("-") * dS(1) + zero * dx(88)
forms = [form_0, form_1, form_2, form_3]

for form in forms:
    scalar_form = dolfinx.fem.form(form)
    local_val = dolfinx.fem.assemble_scalar(scalar_form)
    global_val = mesh.comm.allreduce(local_val, op=MPI.SUM)
    if mesh.comm.rank == 0:
        print(str(form), global_val, flush=True)

On main:

mpirun -n 3 python3 inner_facets.py 
{ (f)('+') } * dS(<Mesh #0>[1], {}) 0.6250000000000003
{ (f)('-') } * dS(<Mesh #0>[1], {}) 0.3750000000000003
{ c_0 } * dx(<Mesh #0>[88], {})
  +  { (f)('+') } * dS(<Mesh #0>[1], {}) 0.6250000000000003
{ c_0 } * dx(<Mesh #0>[88], {})
  +  { (f)('-') } * dS(<Mesh #0>[1], {}) 0.3750000000000003

on this branch

mpirun -n 3 python3 inner_facets.py 
{ (f)('+') } * dS(<Mesh #0>[1], {}) 0.6250000000000003
{ (f)('-') } * dS(<Mesh #0>[1], {}) 0.3750000000000003
{ c_0 } * dx(<Mesh #0>[88], {})
  +  { (f)('+') } * dS(<Mesh #0>[1], {}) 0.0
{ c_0 } * dx(<Mesh #0>[88], {})
  +  { (f)('-') } * dS(<Mesh #0>[1], {}) 1.0000000000000004

[garth-wells]

This feels too implicit to me. We already have 'one-sided' integration using ds. Would a generalisation of ds cover the use cases that this PR aims to address?

[jorgensd]

Replaced by: #2269