Fix domain_integral_type_map

tests/firedrake/submesh/test_submesh_assemble.py

@@ -532,3 +532,34 @@ def test_submesh_assemble_quad_triangle():
     a_ref = c_ref * inner(TrialFunction(V_q), TestFunction(V_t)) * ds_q(label_interf)
     A_ref = assemble(a_ref)
     assert np.allclose(A.M[0][1].values, A_ref.M.values)
+
+
+@pytest.mark.parallel(3)
+def test_assemble_parent_coefficient():
+    subdomain_id = 999
+    nx = 4
+    mesh = UnitSquareMesh(2*nx, nx, quadrilateral=True, reorder=False)
+    x, y = SpatialCoordinate(mesh)
+    M = FunctionSpace(mesh, "DG", 0)
+    marker = Function(M).interpolate(conditional(Or(x > 0.5, y > 0.5), 1, 0))
+    mesh = RelabeledMesh(mesh, [marker], [subdomain_id])
+    submesh = Submesh(mesh, mesh.topological_dimension, subdomain_id, ignore_halo=True, reorder=False)
+
+    def expr(m):
+        x = SpatialCoordinate(m)
+        return 1 + dot(x, x)
+
+    Vsub = FunctionSpace(submesh, "CG", 1)
+    vsub = TestFunction(Vsub)
+    usub = TrialFunction(Vsub)
+
+    Q = FunctionSpace(mesh, "DG", 0)
+    q = Function(Q).interpolate(expr(mesh))
+    A = assemble(inner(grad(usub) * q, grad(vsub))*dx(domain=submesh))
+
+    Qsub = FunctionSpace(submesh, "DG", 0)
+    qsub = Function(Qsub).interpolate(expr(submesh))
+    A_ref = assemble(inner(grad(usub) * qsub, grad(vsub))*dx)
+
+    A_ref.petscmat.axpy(-1, A.petscmat)
+    assert np.isclose(A_ref.petscmat.norm(PETSc.NormType.NORM_FROBENIUS), 0)

tsfc/driver.py

@@ -135,15 +135,22 @@ def compile_integral(integral_data, form_data, prefix, parameters, *, diagonal=F
             if coeff in form_data.coefficient_split:
                 coefficient_split[coeff] = form_data.coefficient_split[coeff]
             coefficient_numbers.append(form_data.original_coefficient_positions[i])
+
     mesh = integral_data.domain
     all_meshes = extract_domains(form_data.original_form)
     domain_number = all_meshes.index(mesh)
+    coefficient_meshes = chain.from_iterable(map(extract_domains, coefficients))
+
+    domain_integral_type_map = dict.fromkeys(all_meshes, None)
+    domain_integral_type_map.update(dict.fromkeys(coefficient_meshes, "cell"))
+    domain_integral_type_map.update(integral_data.domain_integral_type_map)
+
     integral_data_info = TSFCIntegralDataInfo(
         domain=integral_data.domain,
         integral_type=integral_data.integral_type,
         subdomain_id=integral_data.subdomain_id,
         domain_number=domain_number,
-        domain_integral_type_map={mesh: integral_data.domain_integral_type_map[mesh] if mesh in integral_data.domain_integral_type_map else None for mesh in all_meshes},
+        domain_integral_type_map=domain_integral_type_map,
         arguments=arguments,
         coefficients=coefficients,
         coefficient_split=coefficient_split,