assembly parts of an operator in Python

pyproject.toml

@@ -21,6 +21,7 @@ classifiers = [
 dependencies = [
     "maturin>=1.7",
     "numpy",
+    "scipy",
     "cffi",
     'patchelf; platform_system == "Linux"',
     "ndelement",

python/bempp/assembly/boundary.py

@@ -8,6 +8,7 @@
 from ndelement.reference_cell import ReferenceCellType
 from enum import Enum
 from _cffi_backend import _CDataBase
+from scipy.sparse import coo_matrix
 
 _dtypes = {
     0: np.float32,
@@ -27,6 +28,20 @@ def _dtype_value(dtype):
     raise TypeError("Invalid data type")
 
 
+def _convert_to_scipy(rs_sparse_mat: _CDataBase, dtype: typing.Type[np.floating]) -> coo_matrix:
+    """Convert a pointer to a sparse matrix in Rust to a SciPy COO matrix."""
+    shape = (_lib.sparse_shape0(rs_sparse_mat), _lib.sparse_shape1(rs_sparse_mat))
+    size = _lib.sparse_data_size(rs_sparse_mat)
+    data = np.empty(size, dtype=dtype)
+    rows = np.empty(size, dtype=np.uintp)
+    cols = np.empty(size, dtype=np.uintp)
+    _lib.sparse_data(rs_sparse_mat, _ffi.cast("void*", data.ctypes.data))
+    _lib.sparse_rows(rs_sparse_mat, _ffi.cast("uintptr_t*", rows.ctypes.data))
+    _lib.sparse_cols(rs_sparse_mat, _ffi.cast("uintptr_t*", cols.ctypes.data))
+    _lib.free_sparse_matrix(rs_sparse_mat)
+    return coo_matrix((data, (rows, cols)), shape=shape)
+
+
 class OperatorType(Enum):
     """Operator type."""
 
@@ -54,6 +69,7 @@ def __del__(self):
     def assemble_into_dense(
         self, trial_space: FunctionSpace, test_space: FunctionSpace
     ) -> npt.NDArray[np.floating]:
+        """Assemble operator into a dense matrix."""
         assert trial_space.dtype == test_space.dtype == self.dtype
         output = np.zeros(
             (test_space.global_size, trial_space.global_size), dtype=self.dtype, order="F"
@@ -66,6 +82,66 @@ def assemble_into_dense(
         )
         return output
 
+    def assemble_singular_into_dense(
+        self, trial_space: FunctionSpace, test_space: FunctionSpace
+    ) -> npt.NDArray[np.floating]:
+        """Assemble the singular part of an operator into a dense matrix."""
+        assert trial_space.dtype == test_space.dtype == self.dtype
+        output = np.zeros(
+            (test_space.global_size, trial_space.global_size), dtype=self.dtype, order="F"
+        )
+        _lib.boundary_assembler_assemble_singular_into_dense(
+            self._rs_assembler,
+            _ffi.cast("void*", output.ctypes.data),
+            trial_space._rs_space,
+            test_space._rs_space,
+        )
+        return output
+
+    def assemble_nonsingular_into_dense(
+        self, trial_space: FunctionSpace, test_space: FunctionSpace
+    ) -> npt.NDArray[np.floating]:
+        """Assemble the non-singular part of an operator into a dense matrix."""
+        assert trial_space.dtype == test_space.dtype == self.dtype
+        output = np.zeros(
+            (test_space.global_size, trial_space.global_size), dtype=self.dtype, order="F"
+        )
+        _lib.boundary_assembler_assemble_nonsingular_into_dense(
+            self._rs_assembler,
+            _ffi.cast("void*", output.ctypes.data),
+            trial_space._rs_space,
+            test_space._rs_space,
+        )
+        return output
+
+    def assemble_singular(
+        self, trial_space: FunctionSpace, test_space: FunctionSpace
+    ) -> coo_matrix:
+        """Assemble the singular part of an operator into a CSR matrix."""
+        assert trial_space.dtype == test_space.dtype == self.dtype
+        return _convert_to_scipy(
+            _lib.boundary_assembler_assemble_singular(
+                self._rs_assembler,
+                trial_space._rs_space,
+                test_space._rs_space,
+            ),
+            self.dtype,
+        )
+
+    def assemble_singular_correction(
+        self, trial_space: FunctionSpace, test_space: FunctionSpace
+    ) -> coo_matrix:
+        """Assemble the singular correction of an operator into a CSR matrix."""
+        assert trial_space.dtype == test_space.dtype == self.dtype
+        return _convert_to_scipy(
+            _lib.boundary_assembler_assemble_singular_correction(
+                self._rs_assembler,
+                trial_space._rs_space,
+                test_space._rs_space,
+            ),
+            self.dtype,
+        )
+
     def set_quadrature_degree(self, cell: ReferenceCellType, degree: int):
         """Set the (non-singular) quadrature degree for a cell type."""
         if not _lib.boundary_assembler_has_quadrature_degree(self._rs_assembler, cell.value):

python/test/test_boundary_assembler.py

@@ -135,3 +135,138 @@ def test_single_layer_sphere0_dp0():
     )
 
     assert np.allclose(mat, from_cl, rtol=1e-4)
+
+
+@pytest.mark.parametrize(
+    "operator",
+    [
+        OperatorType.SingleLayer,
+        OperatorType.DoubleLayer,
+        OperatorType.AdjointDoubleLayer,
+        OperatorType.Hypersingular,
+    ],
+)
+@pytest.mark.parametrize("test_degree", range(3))
+@pytest.mark.parametrize("trial_degree", range(3))
+def test_assemble_singular(operator, test_degree, trial_degree):
+    grid = regular_sphere(0)
+    test_element = create_family(Family.Lagrange, test_degree, Continuity.Discontinuous)
+    test_space = function_space(grid, test_element)
+    trial_element = create_family(Family.Lagrange, trial_degree, Continuity.Discontinuous)
+    trial_space = function_space(grid, trial_element)
+
+    a = create_laplace_assembler(operator)
+    mat = a.assemble_singular(trial_space, test_space)
+
+    dense = a.assemble_into_dense(trial_space, test_space)
+    for i, j, value in zip(mat.row, mat.col, mat.data):
+        assert np.isclose(dense[i, j], value)
+
+
+@pytest.mark.parametrize(
+    "operator",
+    [
+        OperatorType.SingleLayer,
+        OperatorType.DoubleLayer,
+        OperatorType.AdjointDoubleLayer,
+        OperatorType.Hypersingular,
+    ],
+)
+@pytest.mark.parametrize(
+    "test_degree, test_continuity",
+    [
+        (0, Continuity.Discontinuous),
+        (1, Continuity.Discontinuous),
+        (1, Continuity.Standard),
+        (2, Continuity.Standard),
+    ],
+)
+@pytest.mark.parametrize(
+    "trial_degree, trial_continuity",
+    [
+        (0, Continuity.Discontinuous),
+        (1, Continuity.Discontinuous),
+        (1, Continuity.Standard),
+        (2, Continuity.Standard),
+    ],
+)
+def test_assemble_singular_sparse_vs_dense(
+    operator, test_degree, test_continuity, trial_degree, trial_continuity
+):
+    grid = regular_sphere(0)
+    test_element = create_family(Family.Lagrange, test_degree, test_continuity)
+    test_space = function_space(grid, test_element)
+    trial_element = create_family(Family.Lagrange, trial_degree, trial_continuity)
+    trial_space = function_space(grid, trial_element)
+
+    a = create_laplace_assembler(operator)
+    mat = a.assemble_singular(trial_space, test_space)
+
+    dense = a.assemble_singular_into_dense(trial_space, test_space)
+    assert np.allclose(dense, mat.todense())
+
+
+@pytest.mark.parametrize(
+    "operator",
+    [
+        OperatorType.SingleLayer,
+        OperatorType.DoubleLayer,
+        OperatorType.AdjointDoubleLayer,
+        OperatorType.Hypersingular,
+    ],
+)
+@pytest.mark.parametrize("test_degree", range(3))
+@pytest.mark.parametrize("trial_degree", range(3))
+def test_assemble_singular_correction(operator, test_degree, trial_degree):
+    grid = regular_sphere(0)
+    test_element = create_family(Family.Lagrange, test_degree, Continuity.Discontinuous)
+    test_space = function_space(grid, test_element)
+    trial_element = create_family(Family.Lagrange, trial_degree, Continuity.Discontinuous)
+    trial_space = function_space(grid, trial_element)
+
+    a = create_laplace_assembler(operator)
+    a.assemble_singular_correction(trial_space, test_space)
+
+
+@pytest.mark.parametrize(
+    "operator",
+    [
+        OperatorType.SingleLayer,
+        OperatorType.DoubleLayer,
+        OperatorType.AdjointDoubleLayer,
+        OperatorType.Hypersingular,
+    ],
+)
+@pytest.mark.parametrize(
+    "test_degree, test_continuity",
+    [
+        (0, Continuity.Discontinuous),
+        (1, Continuity.Discontinuous),
+        (1, Continuity.Standard),
+        (2, Continuity.Standard),
+    ],
+)
+@pytest.mark.parametrize(
+    "trial_degree, trial_continuity",
+    [
+        (0, Continuity.Discontinuous),
+        (1, Continuity.Discontinuous),
+        (1, Continuity.Standard),
+        (2, Continuity.Standard),
+    ],
+)
+def test_assemble_singular_and_nonsingular(
+    operator, test_degree, test_continuity, trial_degree, trial_continuity
+):
+    grid = regular_sphere(0)
+    test_element = create_family(Family.Lagrange, test_degree, test_continuity)
+    test_space = function_space(grid, test_element)
+    trial_element = create_family(Family.Lagrange, trial_degree, trial_continuity)
+    trial_space = function_space(grid, trial_element)
+
+    a = create_laplace_assembler(operator)
+    mat = a.assemble_singular_into_dense(trial_space, test_space)
+    mat += a.assemble_nonsingular_into_dense(trial_space, test_space)
+
+    mat2 = a.assemble_into_dense(trial_space, test_space)
+    assert np.allclose(mat, mat2)

src/assembly/boundary/assemblers.rs

@@ -528,7 +528,7 @@ impl<
     }
 
     /// Assemble the singular contributions
-    fn assemble_singular<Space: FunctionSpace<T = T> + Sync>(
+    pub(crate) fn assemble_singular<Space: FunctionSpace<T = T> + Sync>(
         &self,
         shape: [usize; 2],
         trial_space: &Space,
@@ -692,7 +692,7 @@ impl<
     /// Assemble the singular correction
     ///
     /// The singular correction is the contribution is the terms for adjacent cells are assembled using an (incorrect) non-singular quadrature rule
-    fn assemble_singular_correction<Space: FunctionSpace<T = T> + Sync>(
+    pub(crate) fn assemble_singular_correction<Space: FunctionSpace<T = T> + Sync>(
         &self,
         shape: [usize; 2],
         trial_space: &Space,
@@ -823,7 +823,7 @@ impl<
             )
     }
     /// Assemble the non-singular contributions into a dense matrix
-    fn assemble_nonsingular<Space: FunctionSpace<T = T> + Sync>(
+    pub(crate) fn assemble_nonsingular<Space: FunctionSpace<T = T> + Sync>(
         &self,
         output: &RawData2D<T>,
         trial_space: &Space,