[BUG]: Packing of codim-0 form not working for mixture of measures

[jorgensd]

Summarize the issue

Issue is that coefficient packing is not restricted to the integral type. Thus if you have a coefficient on a subdomain not present in the cell integral (dx), it is still packed for this domain, where the map is undefined.

How to reproduce the bug

See mwe below. Will add more info when I got time

Minimal Example (Python)

from mpi4py import MPI
import dolfinx
import dolfinx.fem.petsc
import ufl
import numpy as np

def half(x):
    return x[1] <= 0.5 + 1e-14

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 2, 2, dolfinx.mesh.CellType.triangle)

# Split domain in half and set an interface tag of 5
gdim = mesh.geometry.dim
tdim = mesh.topology.dim
fdim = tdim - 1
mesh.topology.create_entities(fdim)
num_facets_local = (
    mesh.topology.index_map(fdim).size_local + mesh.topology.index_map(fdim).num_ghosts
)
facets = np.arange(num_facets_local, dtype=np.int32)
values = np.full_like(facets, 0, dtype=np.int32)

bottom_cells = dolfinx.mesh.locate_entities(mesh, tdim, half)
num_cells_local = (
    mesh.topology.index_map(tdim).size_local + mesh.topology.index_map(tdim).num_ghosts
)
cells = np.full(num_cells_local, 4, dtype=np.int32)
cells[bottom_cells] = 3
ct = dolfinx.mesh.meshtags(
    mesh, tdim, np.arange(num_cells_local, dtype=np.int32), cells
)
all_b_facets = dolfinx.mesh.compute_incident_entities(
    mesh.topology, ct.find(3), tdim, fdim
)
all_t_facets = dolfinx.mesh.compute_incident_entities(
    mesh.topology, ct.find(4), tdim, fdim
)
interface = np.intersect1d(all_b_facets, all_t_facets)
values[interface] = 5

mt = dolfinx.mesh.meshtags(mesh, mesh.topology.dim - 1, facets, values)

submesh_b, submesh_b_to_mesh, b_v_map = dolfinx.mesh.create_submesh(
    mesh, tdim, ct.find(3)
)[0:3]
submesh_t, submesh_t_to_mesh, t_v_map = dolfinx.mesh.create_submesh(
    mesh, tdim, ct.find(4)
)[0:3]
parent_to_sub_b = np.full(num_facets_local, -1, dtype=np.int32)
parent_to_sub_b[submesh_b_to_mesh] = np.arange(len(submesh_b_to_mesh), dtype=np.int32)
parent_to_sub_t = np.full(num_facets_local, -1, dtype=np.int32)
parent_to_sub_t[submesh_t_to_mesh] = np.arange(len(submesh_t_to_mesh), dtype=np.int32)    

# We need to modify the cell maps, as for `dS` integrals of interfaces between submeshes, there is no entity to map to.
# We use the entity on the same side to fix this (as all restrictions are one-sided)

# Hack, as we use one-sided restrictions, pad dS integral with the same entity from the same cell on both sides
mesh.topology.create_connectivity(fdim, tdim)
f_to_c = mesh.topology.connectivity(fdim, tdim)
for facet in mt.find(5):
    cells = f_to_c.links(facet)
    assert len(cells) == 2
    b_map = parent_to_sub_b[cells]
    t_map = parent_to_sub_t[cells]
    assert max(b_map) > -1
    assert max(t_map) > -1
    parent_to_sub_b[cells] = max(b_map)
    parent_to_sub_t[cells] = max(t_map)

entity_maps = {submesh_b._cpp_object: parent_to_sub_b, submesh_t._cpp_object: parent_to_sub_t}

V_0 = dolfinx.fem.functionspace(submesh_b, ("Lagrange", 3))
V_1 = dolfinx.fem.functionspace(submesh_t, ("Lagrange", 2))

c_0 = dolfinx.fem.Function(V_0)
c_1 = dolfinx.fem.Function(V_1)
c_0.name = "c_0"
c_0.x.array[:] = 2
c_1.name = "c_1"
c_1.x.array[:] = 5

v_0 = ufl.TestFunction(V_0)
v_1 = ufl.TestFunction(V_1)

ct = dolfinx.mesh.meshtags(mesh, tdim, submesh_b_to_mesh, np.full_like(submesh_b_to_mesh, 1,dtype=np.int32))
dx = ufl.Measure("dx", domain=mesh, subdomain_data=ct, subdomain_id=1)
dInterface = ufl.Measure("dS", domain=mesh, subdomain_data=mt, subdomain_id=5)
b_res = "+"
t_res = "-"

c_1_r = c_1(t_res)
v_0_r = v_0(b_res)
v_1_r = v_1(t_res)

F = ufl.inner(c_0, v_0) * dx + ufl.inner(c_1_r, v_0_r) * dInterface
F_c = dolfinx.fem.form(F, entity_maps=entity_maps)

print(dolfinx.fem.assemble.pack_coefficients(F_c._cpp_object))
print(dolfinx.fem.petsc.assemble_vector(F_c))

Output (Python)

No response

Version

main branch

DOLFINx git commit

No response

Installation

No response

Additional information

No response

[jpdean]

I think this deals with a similar problem but for the codim 1 case, so we can probably generalise it to deal with this too

[jorgensd]

Slightly updated example to avoid segfault due to over-specifying dx