Interior facet with consistent orientation

[jorgensd]

Reference implementation:

# 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)

made for https://fenicsproject.discourse.group/t/consistent-side-of-interior-boundaries-with-dolfinx/15968/2