From 3e1d84872da402884491d81981b559bddaaefb3e Mon Sep 17 00:00:00 2001 From: Musawer Ahmad Saqif Date: Tue, 30 Jun 2026 14:32:07 -0600 Subject: [PATCH] Validate BCs for Petrov-Galerkin matrix assembly --- cpp/dolfinx/fem/assembler.h | 53 +++++++++++++++++++++++++- python/dolfinx/fem/petsc.py | 8 +++- python/dolfinx/wrappers/petsc.cpp | 10 +++-- python/test/unit/fem/test_assembler.py | 30 +++++++++++++++ 4 files changed, 94 insertions(+), 7 deletions(-) diff --git a/cpp/dolfinx/fem/assembler.h b/cpp/dolfinx/fem/assembler.h index 0654e108e0d..82d8c87324b 100644 --- a/cpp/dolfinx/fem/assembler.h +++ b/cpp/dolfinx/fem/assembler.h @@ -40,6 +40,51 @@ class Form; template class FunctionSpace; +namespace impl +{ +/// @brief Check that Dirichlet boundary conditions are supported for matrix +/// assembly when test and trial spaces are distinct but share a dofmap. +/// @param[in] a Bilinear form. +/// @param[in] bcs Boundary conditions to check. +/// @throws std::runtime_error If a boundary condition is applied to either +/// space of a form whose distinct test and trial spaces share the same dofmap. +template +void check_bcs_for_dofmap_sharing_spaces( + const Form& a, + const std::vector>>& bcs) +{ + if (bcs.empty() or a.function_spaces().size() != 2) + return; + + const auto& V0 = a.function_spaces().at(0); + const auto& V1 = a.function_spaces().at(1); + assert(V0); + assert(V1); + if (V0 == V1) + return; + if (V0->mesh() != V1->mesh()) + return; + if (V0->dofmaps().size() != 1 or V1->dofmaps().size() != 1) + return; + if (V0->dofmap() != V1->dofmap()) + return; + if (V0->element()->signature() != V1->element()->signature()) + return; + + for (auto& bc : bcs) + { + assert(bc.get().function_space()); + if (V0->contains(*bc.get().function_space()) + or V1->contains(*bc.get().function_space())) + { + throw std::runtime_error( + "DirichletBCs on distinct test and trial spaces that share a " + "dofmap are not supported by this matrix assembly function."); + } + } +} +} // namespace impl + /// @brief Evaluate an Expression on cells or facets. /// /// This function accepts packed coefficient data, which allows it be @@ -497,13 +542,19 @@ void assemble_matrix( /// @param[in] coefficients Coefficients that appear in `a`. /// @param[in] bcs Boundary conditions to apply. For boundary condition /// dofs the row and column are zeroed. The diagonal entry is not set. +/// @param[in] check_bcs Check that the boundary conditions are supported by +/// this assembly function. template void assemble_matrix( auto mat_add, const Form& a, std::span constants, const std::map, std::pair, int>>& coefficients, - const std::vector>>& bcs) + const std::vector>>& bcs, + bool check_bcs = true) { + if (check_bcs) + impl::check_bcs_for_dofmap_sharing_spaces(a, bcs); + // Index maps for dof ranges // NOTE: For mixed-topology meshes, there will be multiple DOF maps, // but the index maps are the same. diff --git a/python/dolfinx/fem/petsc.py b/python/dolfinx/fem/petsc.py index e5dee8c5cf1..956c9205ccc 100644 --- a/python/dolfinx/fem/petsc.py +++ b/python/dolfinx/fem/petsc.py @@ -412,6 +412,8 @@ def _( | Sequence[Sequence[dict[tuple[IntegralType, int], npt.NDArray]]] | None ) = None, + *, + _check_bcs: bool = True, ) -> PETSc.Mat: # type: ignore[name-defined] """Assemble bilinear form into a matrix. @@ -422,6 +424,7 @@ def _( The returned matrix is not finalised, i.e. ghost values are not accumulated. """ + bcs = [] if bcs is None else bcs if A.getType() == PETSc.Mat.Type.NEST: # type: ignore[attr-defined] if not isinstance(a, Sequence): raise ValueError("Must provide a sequence of forms when assembling a nest matrix") @@ -431,7 +434,7 @@ def _( for j, (a_block, const, coeff) in enumerate(zip(a_row, const_row, coeff_row)): if a_block is not None: Asub = A.getNestSubMatrix(i, j) - assemble_matrix(Asub, a_block, bcs, diag, const, coeff) # type: ignore[arg-type] + assemble_matrix(Asub, a_block, bcs, diag, const, coeff, _check_bcs=False) # type: ignore[arg-type,call-arg] elif i == j: for bc in bcs: row_forms = [row_form for row_form in a_row if row_form is not None] @@ -474,6 +477,7 @@ def _( coeffs[i][j], # type: ignore[index] _bcs, True, + False, ) A.restoreLocalSubMatrix(is0[i], is1[j], Asub) elif i == j: @@ -502,7 +506,7 @@ def _( constants = pack_constants(a) if constants is None else constants # type: ignore[assignment] coeffs = pack_coefficients(a) if coeffs is None else coeffs # type: ignore[assignment] _bcs = [bc._cpp_object for bc in bcs] if bcs is not None else [] - _cpp.fem.petsc.assemble_matrix(A, a._cpp_object, constants, coeffs, _bcs) + _cpp.fem.petsc.assemble_matrix(A, a._cpp_object, constants, coeffs, _bcs, False, _check_bcs) if a.function_spaces[0] is a.function_spaces[1]: A.assemblyBegin(PETSc.Mat.AssemblyType.FLUSH) # type: ignore[attr-defined] A.assemblyEnd(PETSc.Mat.AssemblyType.FLUSH) # type: ignore[attr-defined] diff --git a/python/dolfinx/wrappers/petsc.cpp b/python/dolfinx/wrappers/petsc.cpp index c343b769889..4a60afed7fe 100644 --- a/python/dolfinx/wrappers/petsc.cpp +++ b/python/dolfinx/wrappers/petsc.cpp @@ -186,7 +186,7 @@ void petsc_fem_module(nb::module_& m) nb::c_contig>>& coefficients, const std::vector< const dolfinx::fem::DirichletBC*>& bcs, - bool unrolled) + bool unrolled, bool check_bcs) { std::vector>> @@ -204,18 +204,20 @@ void petsc_fem_module(nb::module_& m) a.function_spaces()[1]->dofmap()->bs(), ADD_VALUES); dolfinx::fem::assemble_matrix( set_fn, a, std::span(constants.data(), constants.size()), - dolfinx_wrappers::py_to_cpp_coeffs(coefficients), _bcs); + dolfinx_wrappers::py_to_cpp_coeffs(coefficients), _bcs, + check_bcs); } else { dolfinx::fem::assemble_matrix( dolfinx::la::petsc::Matrix::set_block_fn(A, ADD_VALUES), a, std::span(constants.data(), constants.size()), - dolfinx_wrappers::py_to_cpp_coeffs(coefficients), _bcs); + dolfinx_wrappers::py_to_cpp_coeffs(coefficients), _bcs, + check_bcs); } }, nb::arg("A"), nb::arg("a"), nb::arg("constants"), nb::arg("coeffs"), - nb::arg("bcs"), nb::arg("unrolled") = false, + nb::arg("bcs"), nb::arg("unrolled") = false, nb::arg("check_bcs") = true, "Assemble bilinear form into an existing PETSc matrix"); m.def( "assemble_matrix", diff --git a/python/test/unit/fem/test_assembler.py b/python/test/unit/fem/test_assembler.py index 55d77b88a51..213198080f1 100644 --- a/python/test/unit/fem/test_assembler.py +++ b/python/test/unit/fem/test_assembler.py @@ -178,6 +178,26 @@ def nest_matrix_norm(A): return math.sqrt(norm) +def _shared_dofmap_distinct_space_form_and_bc(): + msh = create_unit_square(MPI.COMM_WORLD, 4, 4) + U = functionspace(msh, ("Lagrange", 1)) + V = U.clone() + u = ufl.TrialFunction(U) + v = ufl.TestFunction(V) + a = form(inner(u, v) * dx) + + facets = locate_entities_boundary(msh, msh.topology.dim - 1, lambda x: np.isclose(x[0], 0.0)) + dofs = locate_dofs_topological(V, msh.topology.dim - 1, facets) + bc = dirichletbc(default_scalar_type(0), dofs, V) + return a, bc + + +def test_assembly_rejects_bcs_for_distinct_spaces_sharing_dofmap(): + a, bc = _shared_dofmap_distinct_space_form_and_bc() + with pytest.raises(RuntimeError, match="distinct test and trial spaces that share a dofmap"): + assemble_matrix(a, bcs=[bc]) + + @pytest.mark.petsc4py def test_vector_single_space_as_block(): from dolfinx.fem.petsc import create_vector as petsc_create_vector @@ -193,6 +213,16 @@ def test_vector_single_space_as_block(): class TestPETScAssemblers: """Test PETSc-based assemblers for matrices and vectors.""" + def test_assembly_rejects_bcs_for_distinct_spaces_sharing_dofmap(self): + """Test that PETSc assembly rejects ambiguous standalone BCs.""" + from dolfinx.fem.petsc import assemble_matrix as petsc_assemble_matrix + + a, bc = _shared_dofmap_distinct_space_form_and_bc() + with pytest.raises( + RuntimeError, match="distinct test and trial spaces that share a dofmap" + ): + petsc_assemble_matrix(a, bcs=[bc]) + @pytest.mark.parametrize("mode", [GhostMode.none, GhostMode.shared_facet]) def test_basic_assembly_petsc_matrixcsr(self, mode): """Test basic assembly of PETSc Mat and compare with MatrixCSR assembly."""