Guard cloned-space BC assembly in C++

cpp/dolfinx/fem/assembler.h

@@ -40,6 +40,51 @@ class Form;
 template <std::floating_point T>
 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 <dolfinx::scalar T, std::floating_point U>
+void check_bcs_for_dofmap_sharing_spaces(
+    const Form<T, U>& a,
+    const std::vector<std::reference_wrapper<const DirichletBC<T, U>>>& 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 <dolfinx::scalar T, std::floating_point U>
 void assemble_matrix(
     auto mat_add, const Form<T, U>& a, std::span<const T> constants,
     const std::map<std::pair<IntegralType, int>,
                    std::pair<std::span<const T>, int>>& coefficients,
-    const std::vector<std::reference_wrapper<const DirichletBC<T, U>>>& bcs)
+    const std::vector<std::reference_wrapper<const DirichletBC<T, U>>>& 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.

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]

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<PetscScalar, PetscReal>*>& bcs,
-         bool unrolled)
+         bool unrolled, bool check_bcs)
       {
         std::vector<std::reference_wrapper<
             const dolfinx::fem::DirichletBC<PetscScalar, PetscReal>>>
@@ -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",

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."""