From fa7037e2a9423bf9c48e27f8d34aa671381884c6 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Fri, 12 Dec 2025 13:55:07 +0000 Subject: [PATCH] Fix domain_integral_type_map --- .../submesh/test_submesh_assemble.py | 31 +++++++++++++++++++ tsfc/driver.py | 9 +++++- 2 files changed, 39 insertions(+), 1 deletion(-) diff --git a/tests/firedrake/submesh/test_submesh_assemble.py b/tests/firedrake/submesh/test_submesh_assemble.py index 28f244581f..87f4db35dd 100644 --- a/tests/firedrake/submesh/test_submesh_assemble.py +++ b/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) diff --git a/tsfc/driver.py b/tsfc/driver.py index c82a4eb5ed..3f30038db2 100644 --- a/tsfc/driver.py +++ b/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,