Assembly of combination of facet and cell functions

dolfinx
1

I’m assembling interface contributions to two (or more) submeshes, closely following the helpful example from Jorgen Dokken:

problem.py for working with submeshes

There, the forms containing the interface contributions were in terms of restrictions on the boundary of the “bulk” functions. Now I want to also have them depend on a function defined directly on that interface.

As an example, I created a unit square mesh, divided into two submeshes along x=0.5 and defined a parameter function along the interface, p = y. Now I want to assemble ∫ v_Γ p (1 - u) dS. Creating the functionspace for the interface and defining this functions worked fine, but its contribution are not correctly assembled. Here is my MWE:

#!/usr/bin/env python3

from mpi4py import MPI
import numpy as np

from dolfinx.mesh import (
    create_unit_square,
    locate_entities,
    locate_entities_boundary,
    create_submesh,
    compute_incident_entities,
    meshtags
)
from dolfinx.fem import (
    form,
    functionspace,
    Function,
)

from dolfinx.fem.petsc import assemble_vector

from ufl import (
    Measure,
    TestFunction,
)

# Create main mesh
mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)

# cell and facet indices
tdim = mesh.topology.dim
fdim = mesh.topology.dim - 1
left_cells = locate_entities(mesh, tdim, lambda x: x[0] <= 0.5)
right_cells = locate_entities(mesh, tdim, lambda x: x[0] >= 0.5)
mesh.topology.create_connectivity(fdim, tdim)
left_facets = compute_incident_entities(mesh.topology, left_cells, tdim, fdim)
right_facets = compute_incident_entities(mesh.topology, right_cells, tdim, fdim)
interface = np.intersect1d(left_facets, right_facets)

# create submeshes
submesh_left, cell_map_left, _, _ = create_submesh(mesh, tdim, left_cells)

# facet tag = 3 at interface
num_facets_local = (
    mesh.topology.index_map(fdim).size_local + mesh.topology.index_map(fdim).num_ghosts
)
facet_indices = np.arange(num_facets_local, dtype=np.int32)
values = np.zeros(num_facets_local, dtype=np.int32)
values[interface] = 3
facet_tags = meshtags(mesh, fdim, facet_indices, values)

# create entity maps
num_cells_local = (
    mesh.topology.index_map(tdim).size_local + mesh.topology.index_map(tdim).num_ghosts
)
mesh_to_left = np.full(num_cells_local, -1, dtype=np.int32)
mesh_to_left[cell_map_left] = np.arange(len(cell_map_left), dtype=np.int32)

f_to_c = mesh.topology.connectivity(fdim, tdim)
for facet in facet_tags.find(3):
    cells = f_to_c.links(facet)
    assert len(cells) == 2
    a_map = mesh_to_left[cells]
    mesh_to_left[cells] = max(a_map)

interface_mesh, cell_map_interface, _, _ = create_submesh(mesh, fdim, interface)
inverse_facet_map = np.full(num_facets_local, -1, dtype=np.int32)
inverse_facet_map[cell_map_interface] = np.arange(len(cell_map_interface), dtype=np.int32)

entity_maps = {submesh_left: mesh_to_left, interface_mesh: inverse_facet_map}

# Create measures
cell_tags_left = meshtags(mesh, tdim, left_cells, 1)
dx_l = Measure("dx", domain=mesh, subdomain_data=cell_tags_left, subdomain_id=1)
dS = Measure("dS", domain=mesh, subdomain_data=facet_tags, subdomain_id=3)

# create functions
V_l = functionspace(submesh_left, ("Lagrange", 1))
V_i = functionspace(interface_mesh, ("Lagrange", 1))

p = Function(V_i)
p.interpolate(lambda x: x[1])

v_l = TestFunction(V_l)
u_l = Function(V_l)
u_l.vector.array.fill(0.0)

# forms
F_l = v_l("+") * (1.0 - u_l("+")) * dS
F_form = form(F_l, entity_maps=entity_maps)
F_vec = assemble_vector(F_form)

F_l2 = v_l("+") * p * (1.0 - u_l("+")) * dS
F_form2 = form(F_l2, entity_maps=entity_maps)
F_vec2 = assemble_vector(F_form2)

np.set_printoptions(precision=3)
print(f"Values of facet fct (as expected):\n{p.vector.array}")
print(f"Values of facet integral without p (as expected):\n{F_vec.array}")
print(f"Values of facet integral with p (unexpected):\n{F_vec2.array}")

Any hint on what I’m doing wrong is greatly appreciated.