Repository: bempp/bempp-rs (language: python)

Issue:
Add Python wrapper
- [x] FunctionSpace
    - [x] `grid`
    - [x] `element`
    - [x] `is_serial`
    - [x] `get_local_dof_numbers`
    - [x] `local_size`
    - [x] `global_size`
    - [x] `cell_dofs`
    - ~[ ] `cell_colouring`~
    - [x] `global_dof_index`
    - [x] `ownership`
- [x] Assembly
    - [x] Potential
        - [x] Options
            - [x] `quadrature_degrees`
            - [x] `batch_size`
        - [x] `assemble_into_dense`
    - [x] Boundary
        - [x] Options
            - [x] `quadrature_degrees`
            - [x] `singular_quadrature_degrees`
            - [x] `batch_size`
        - [x] `assemble_into_dense`
        - [x] `assemble_into_csr`


Discussion / hints:

Fix patch:
diff --git a/.github/workflows/style-checks.yml b/.github/workflows/style-checks.yml
index 45dc3788..c008f71a 100644
--- a/.github/workflows/style-checks.yml
+++ b/.github/workflows/style-checks.yml
@@ -37,3 +37,17 @@ jobs:
           cargo clippy ${{ matrix.feature-flags }} -- -D warnings
           cargo clippy --tests ${{ matrix.feature-flags }} -- -D warnings
           cargo clippy --examples ${{ matrix.feature-flags }} -- -D warnings
+
+      - name: Set up Python
+        uses: actions/setup-python@v4
+        with:
+          python-version: "3.12"
+      - run: |
+          python3 -m pip install ruff
+          python3 -m ruff check python
+          python3 -m ruff format --check python
+        name: Run ruff checks
+      - run: |
+          python3 -m pip install mypy
+          python3 -m mypy python
+        name: Run mypy checks
diff --git a/benches/assembly_benchmark.rs b/benches/assembly_benchmark.rs
index 26ea9d2c..6c3f4c05 100644
--- a/benches/assembly_benchmark.rs
+++ b/benches/assembly_benchmark.rs
@@ -20,12 +20,12 @@ pub fn assembly_parts_benchmark(c: &mut Criterion) {
 
         let colouring = space.cell_colouring();
         let mut a = BoundaryAssembler::<f64, _, _>::new_laplace_single_layer();
-        a.quadrature_degree(ReferenceCellType::Triangle, 16);
-        a.singular_quadrature_degree(
+        a.set_quadrature_degree(ReferenceCellType::Triangle, 16);
+        a.set_singular_quadrature_degree(
             (ReferenceCellType::Triangle, ReferenceCellType::Triangle),
             4,
         );
-        a.batch_size(128);
+        a.set_batch_size(128);
 
         group.bench_function(
             &format!(
diff --git a/examples/assembly.rs b/examples/assembly.rs
index 586eaa09..0f94e8d4 100644
--- a/examples/assembly.rs
+++ b/examples/assembly.rs
@@ -21,11 +21,11 @@ fn main() {
 
     // Adjust the quadrature degree for non-singular integrals on a triangle.
     // This makes the integrals use a quadrature rule with 16 points
-    a.quadrature_degree(ReferenceCellType::Triangle, 16);
+    a.set_quadrature_degree(ReferenceCellType::Triangle, 16);
 
     // Adjust the quadrature degree for singular integrals on a pair ortriangles.
     // This makes the integrals use a quadrature rule based on a rule on an interval with 4 points
-    a.singular_quadrature_degree(
+    a.set_singular_quadrature_degree(
         (ReferenceCellType::Triangle, ReferenceCellType::Triangle),
         4,
     );
diff --git a/pyproject.toml b/pyproject.toml
index 318e317c..5c1b192b 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -21,12 +21,13 @@ classifiers = [
 dependencies = [
     "maturin>=1.7",
     "numpy",
+    "scipy",
     "cffi",
     'patchelf; platform_system == "Linux"',
     "ndelement",
     "ndgrid"
 ]
-packages = ["bempp"]
+packages = ["bempp", "bempp.assembly"]
 
 [project.urls]
 homepage = "https://github.com/bempp/bempp-rs"
diff --git a/python/bempp/__init__.py b/python/bempp/__init__.py
index 1cffbb27..6b92b118 100644
--- a/python/bempp/__init__.py
+++ b/python/bempp/__init__.py
@@ -1,3 +1,3 @@
 """Bempp."""
 
-from bempp import function_space
+from bempp import function_space, assembly
diff --git a/python/bempp/assembly/__init__.py b/python/bempp/assembly/__init__.py
new file mode 100644
index 00000000..205f2856
--- /dev/null
+++ b/python/bempp/assembly/__init__.py
@@ -0,0 +1,3 @@
+"""Assembly."""
+
+from bempp.assembly import boundary
diff --git a/python/bempp/assembly/boundary.py b/python/bempp/assembly/boundary.py
new file mode 100644
index 00000000..301d4278
--- /dev/null
+++ b/python/bempp/assembly/boundary.py
@@ -0,0 +1,205 @@
+"""Boundary assembly."""
+
+import typing
+import numpy as np
+import numpy.typing as npt
+from bempp._bempprs import lib as _lib, ffi as _ffi
+from bempp.function_space import FunctionSpace
+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,
+    1: np.float64,
+}
+_ctypes = {
+    np.float32: "float",
+    np.float64: "double",
+}
+
+
+def _dtype_value(dtype):
+    """Get the u8 identifier for a data type."""
+    for i, j in _dtypes.items():
+        if j == dtype:
+            return i
+    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."""
+
+    SingleLayer = 0
+    DoubleLayer = 1
+    AdjointDoubleLayer = 2
+    Hypersingular = 3
+    ElectricField = 4
+    MagneticField = 5
+
+
+class BoundaryAssembler(object):
+    """Boundary assembler."""
+
+    def __init__(self, rs_assembler: _CDataBase, owned: bool = True):
+        """Initialise."""
+        self._rs_assembler = rs_assembler
+        self._owned = owned
+
+    def __del__(self):
+        """Delete."""
+        if self._owned:
+            _lib.free_boundary_assembler(self._rs_assembler)
+
+    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"
+        )
+        _lib.boundary_assembler_assemble_into_dense(
+            self._rs_assembler,
+            _ffi.cast("void*", output.ctypes.data),
+            trial_space._rs_space,
+            test_space._rs_space,
+        )
+        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):
+            raise ValueError(f"Invalid cell type: {cell}")
+        _lib.boundary_assembler_set_quadrature_degree(self._rs_assembler, cell.value, degree)
+
+    def quadrature_degree(self, cell: ReferenceCellType) -> int:
+        """Get the (non-singular) quadrature degree for a cell type."""
+        if not _lib.boundary_assembler_has_quadrature_degree(self._rs_assembler, cell.value):
+            raise ValueError(f"Invalid cell type: {cell}")
+        return _lib.boundary_assembler_quadrature_degree(self._rs_assembler, cell.value)
+
+    def set_singular_quadrature_degree(
+        self, cell0: ReferenceCellType, cell1: ReferenceCellType, degree: int
+    ):
+        """Set the singular quadrature degree for a cell type."""
+        if not _lib.boundary_assembler_has_singular_quadrature_degree(
+            self._rs_assembler, cell0.value, cell1.value
+        ):
+            raise ValueError(f"Invalid cell pair: {cell0} and {cell1}")
+        _lib.boundary_assembler_set_singular_quadrature_degree(
+            self._rs_assembler, cell0.value, cell1.value, degree
+        )
+
+    def singular_quadrature_degree(self, cell0: ReferenceCellType, cell1: ReferenceCellType) -> int:
+        """Get the singular_quadrature degree for a cell type."""
+        if not _lib.boundary_assembler_has_singular_quadrature_degree(
+            self._rs_assembler, cell0.value, cell1.value
+        ):
+            raise ValueError(f"Invalid cell pair: {cell0} and {cell1}")
+        return _lib.boundary_assembler_singular_quadrature_degree(
+            self._rs_assembler, cell0.value, cell1.value
+        )
+
+    def set_batch_size(self, batch_size: int):
+        """Set the batch size."""
+        _lib.boundary_assembler_set_batch_size(self._rs_assembler, batch_size)
+
+    @property
+    def batch_size(self) -> int:
+        """Get the batch size."""
+        return _lib.boundary_assembler_batch_size(self._rs_assembler)
+
+    @property
+    def dtype(self):
+        """Data type."""
+        return _dtypes[_lib.boundary_assembler_dtype(self._rs_assembler)]
+
+    @property
+    def _ctype(self):
+        """C data type."""
+        return _ctypes[self.dtype]
+
+
+def create_laplace_assembler(
+    operator_type: OperatorType, dtype: typing.Type[np.floating] = np.float64
+):
+    """Create a Laplace assembler."""
+    return BoundaryAssembler(
+        _lib.laplace_boundary_assembler_new(operator_type.value, _dtype_value(dtype))
+    )
diff --git a/python/bempp/function_space.py b/python/bempp/function_space.py
index 9fe091e7..272237ed 100644
--- a/python/bempp/function_space.py
+++ b/python/bempp/function_space.py
@@ -9,6 +9,7 @@
 from ndgrid.ownership import Owned, Ghost
 from ndelement.ciarlet import ElementFamily, CiarletElement
 from ndelement.reference_cell import ReferenceCellType
+from _cffi_backend import _CDataBase
 
 _dtypes = {
     0: np.float32,
@@ -23,13 +24,15 @@
 class FunctionSpace(object):
     """Function space."""
 
-    def __init__(self, rs_space):
+    def __init__(self, rs_space: _CDataBase, owned: bool = True):
         """Initialise."""
         self._rs_space = rs_space
+        self._owned = True
 
     def __del__(self):
         """Delete."""
-        _lib.free_space(self._rs_space)
+        if self._owned:
+            _lib.free_space(self._rs_space)
 
     def element(self, entity: ReferenceCellType) -> CiarletElement:
         """Get the grid that this space is defined on."""
diff --git a/src/assembly/boundary.rs b/src/assembly/boundary.rs
index 6162d01b..e3ab7598 100644
--- a/src/assembly/boundary.rs
+++ b/src/assembly/boundary.rs
@@ -1,7 +1,7 @@
 //! Assembly of boundary operators
 mod assemblers;
 pub(crate) mod cell_pair_assemblers;
-mod integrands;
+pub mod integrands;
 
 pub use assemblers::BoundaryAssembler;
 
diff --git a/src/assembly/boundary/assemblers.rs b/src/assembly/boundary/assemblers.rs
index 7453a0c7..78f298df 100644
--- a/src/assembly/boundary/assemblers.rs
+++ b/src/assembly/boundary/assemblers.rs
@@ -484,12 +484,17 @@ impl<
     }
 
     /// Set (non-singular) quadrature degree for a cell type
-    pub fn quadrature_degree(&mut self, cell: ReferenceCellType, degree: usize) {
+    pub fn set_quadrature_degree(&mut self, cell: ReferenceCellType, degree: usize) {
         *self.options.quadrature_degrees.get_mut(&cell).unwrap() = degree;
     }
 
+    /// Get (non-singular) quadrature degree for a cell type
+    pub fn quadrature_degree(&self, cell: ReferenceCellType) -> Option<usize> {
+        self.options.quadrature_degrees.get(&cell).copied()
+    }
+
     /// Set singular quadrature degree for a pair of cell types
-    pub fn singular_quadrature_degree(
+    pub fn set_singular_quadrature_degree(
         &mut self,
         cells: (ReferenceCellType, ReferenceCellType),
         degree: usize,
@@ -501,13 +506,29 @@ impl<
             .unwrap() = degree;
     }
 
+    /// Get singular quadrature degree for a pair of cell types
+    pub fn singular_quadrature_degree(
+        &self,
+        cells: (ReferenceCellType, ReferenceCellType),
+    ) -> Option<usize> {
+        self.options
+            .singular_quadrature_degrees
+            .get(&cells)
+            .copied()
+    }
+
     /// Set the maximum size of a batch of cells to send to an assembly function
-    pub fn batch_size(&mut self, size: usize) {
+    pub fn set_batch_size(&mut self, size: usize) {
         self.options.batch_size = size;
     }
 
+    /// Get the maximum size of a batch of cells to send to an assembly function
+    pub fn batch_size(&self) -> usize {
+        self.options.batch_size
+    }
+
     /// 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,
@@ -671,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,
@@ -802,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,
@@ -937,9 +958,12 @@ impl<
     > BoundaryAssembly for BoundaryAssembler<T, Integrand, Kernel>
 {
     type T = T;
-    fn assemble_singular_into_dense<Space: FunctionSpace<T = T> + Sync>(
+    fn assemble_singular_into_dense<
+        Space: FunctionSpace<T = T> + Sync,
+        Array2: RandomAccessMut<2, Item = T> + Shape<2> + RawAccessMut<Item = T>,
+    >(
         &self,
-        output: &mut RlstArray<T, 2>,
+        output: &mut Array2,
         trial_space: &Space,
         test_space: &Space,
     ) {
@@ -969,9 +993,12 @@ impl<
         .unwrap()
     }
 
-    fn assemble_singular_correction_into_dense<Space: FunctionSpace<T = T> + Sync>(
+    fn assemble_singular_correction_into_dense<
+        Space: FunctionSpace<T = T> + Sync,
+        Array2: RandomAccessMut<2, Item = T> + Shape<2> + RawAccessMut<Item = T>,
+    >(
         &self,
-        output: &mut RlstArray<T, 2>,
+        output: &mut Array2,
         trial_space: &Space,
         test_space: &Space,
     ) {
@@ -1002,9 +1029,12 @@ impl<
         .unwrap()
     }
 
-    fn assemble_into_dense<Space: FunctionSpace<T = T> + Sync>(
+    fn assemble_into_dense<
+        Space: FunctionSpace<T = T> + Sync,
+        Array2: RandomAccessMut<2, Item = T> + Shape<2> + RawAccessMut<Item = T>,
+    >(
         &self,
-        output: &mut RlstArray<T, 2>,
+        output: &mut Array2,
         trial_space: &Space,
         test_space: &Space,
     ) {
@@ -1021,9 +1051,12 @@ impl<
         self.assemble_singular_into_dense(output, trial_space, test_space);
     }
 
-    fn assemble_nonsingular_into_dense<Space: FunctionSpace<T = T> + Sync>(
+    fn assemble_nonsingular_into_dense<
+        Space: FunctionSpace<T = T> + Sync,
+        Array2: RandomAccessMut<2, Item = T> + Shape<2> + RawAccessMut<Item = T>,
+    >(
         &self,
-        output: &mut RlstArray<T, 2>,
+        output: &mut Array2,
         trial_space: &Space,
         test_space: &Space,
         trial_colouring: &HashMap<ReferenceCellType, Vec<Vec<usize>>>,
diff --git a/src/assembly/boundary/integrands.rs b/src/assembly/boundary/integrands.rs
index 8b552761..18356253 100644
--- a/src/assembly/boundary/integrands.rs
+++ b/src/assembly/boundary/integrands.rs
@@ -27,6 +27,7 @@ pub struct BoundaryIntegrandSum<
 impl<T: RlstScalar, I0: BoundaryIntegrand<T = T>, I1: BoundaryIntegrand<T = T>>
     BoundaryIntegrandSum<T, I0, I1>
 {
+    /// Create new
     pub fn new(integrand0: I0, integrand1: I1) -> Self {
         Self {
             integrand0,
@@ -63,6 +64,7 @@ pub struct BoundaryIntegrandScalarProduct<T: RlstScalar, I: BoundaryIntegrand<T
 }
 
 impl<T: RlstScalar, I: BoundaryIntegrand<T = T>> BoundaryIntegrandScalarProduct<T, I> {
+    /// Create new
     pub fn new(scalar: T, integrand: I) -> Self {
         Self { scalar, integrand }
     }
diff --git a/src/assembly/boundary/integrands/adjoint_double_layer.rs b/src/assembly/boundary/integrands/adjoint_double_layer.rs
index 1cd6d0c9..ea2fcd97 100644
--- a/src/assembly/boundary/integrands/adjoint_double_layer.rs
+++ b/src/assembly/boundary/integrands/adjoint_double_layer.rs
@@ -2,11 +2,13 @@
 use crate::traits::{Access1D, Access2D, BoundaryIntegrand, GeometryAccess};
 use rlst::RlstScalar;
 
+/// Integrand for an adjoint double layer boundary operator
 pub struct AdjointDoubleLayerBoundaryIntegrand<T: RlstScalar> {
     _t: std::marker::PhantomData<T>,
 }
 
 impl<T: RlstScalar> AdjointDoubleLayerBoundaryIntegrand<T> {
+    /// Create new
     pub fn new() -> Self {
         Self {
             _t: std::marker::PhantomData,
diff --git a/src/assembly/boundary/integrands/double_layer.rs b/src/assembly/boundary/integrands/double_layer.rs
index 01979a43..cc1b2fed 100644
--- a/src/assembly/boundary/integrands/double_layer.rs
+++ b/src/assembly/boundary/integrands/double_layer.rs
@@ -2,11 +2,13 @@
 use crate::traits::{Access1D, Access2D, BoundaryIntegrand, GeometryAccess};
 use rlst::RlstScalar;
 
+/// Integrand for a double layer boundary operator
 pub struct DoubleLayerBoundaryIntegrand<T: RlstScalar> {
     _t: std::marker::PhantomData<T>,
 }
 
 impl<T: RlstScalar> DoubleLayerBoundaryIntegrand<T> {
+    /// Create new
     pub fn new() -> Self {
         Self {
             _t: std::marker::PhantomData,
diff --git a/src/assembly/boundary/integrands/hypersingular.rs b/src/assembly/boundary/integrands/hypersingular.rs
index ae77f61c..478bdf2f 100644
--- a/src/assembly/boundary/integrands/hypersingular.rs
+++ b/src/assembly/boundary/integrands/hypersingular.rs
@@ -2,11 +2,13 @@
 use crate::traits::{Access1D, Access2D, BoundaryIntegrand, GeometryAccess};
 use rlst::RlstScalar;
 
+/// Integrand for the curl curl term of a hypersingular boundary operator
 pub struct HypersingularCurlCurlBoundaryIntegrand<T: RlstScalar> {
     _t: std::marker::PhantomData<T>,
 }
 
 impl<T: RlstScalar> HypersingularCurlCurlBoundaryIntegrand<T> {
+    /// Create new
     pub fn new() -> Self {
         Self {
             _t: std::marker::PhantomData,
@@ -52,11 +54,13 @@ unsafe impl<T: RlstScalar> BoundaryIntegrand for HypersingularCurlCurlBoundaryIn
     }
 }
 
+/// Integrand for the normal normal term of a hypersingular boundary operator
 pub struct HypersingularNormalNormalBoundaryIntegrand<T: RlstScalar> {
     _t: std::marker::PhantomData<T>,
 }
 
 impl<T: RlstScalar> HypersingularNormalNormalBoundaryIntegrand<T> {
+    /// Create new
     pub fn new() -> Self {
         Self {
             _t: std::marker::PhantomData,
diff --git a/src/assembly/boundary/integrands/single_layer.rs b/src/assembly/boundary/integrands/single_layer.rs
index 600c1b05..458477ac 100644
--- a/src/assembly/boundary/integrands/single_layer.rs
+++ b/src/assembly/boundary/integrands/single_layer.rs
@@ -2,11 +2,13 @@
 use crate::traits::{Access1D, Access2D, BoundaryIntegrand, GeometryAccess};
 use rlst::RlstScalar;
 
+/// Integrand for a single layer boundary operator
 pub struct SingleLayerBoundaryIntegrand<T: RlstScalar> {
     _t: std::marker::PhantomData<T>,
 }
 
 impl<T: RlstScalar> SingleLayerBoundaryIntegrand<T> {
+    /// Create new
     pub fn new() -> Self {
         Self {
             _t: std::marker::PhantomData,
diff --git a/src/bindings.rs b/src/bindings.rs
index dbcabf12..ab48e835 100644
--- a/src/bindings.rs
+++ b/src/bindings.rs
@@ -12,14 +12,19 @@ pub enum DType {
     C64 = 3,
 }
 
-#[derive(Debug, PartialEq, Clone, Copy)]
-#[repr(u8)]
-pub enum RealDType {
-    F32 = 0,
-    F64 = 1,
+impl DType {
+    fn from(value: u8) -> Option<Self> {
+        match value {
+            0 => Some(DType::F32),
+            1 => Some(DType::F64),
+            2 => Some(DType::C32),
+            3 => Some(DType::C64),
+            _ => None,
+        }
+    }
 }
 
-mod function {
+pub mod function {
     use super::DType;
     use crate::{function::SerialFunctionSpace, traits::FunctionSpace};
     use ndelement::{
@@ -663,3 +668,3770 @@ mod function {
         }
     }
 }
+
+pub mod assembly {
+    use super::function::{extract_space, FunctionSpaceWrapper, GridType, SpaceType};
+    use super::DType;
+    use crate::assembly::common::SparseMatrixData;
+    use crate::{
+        assembly::boundary::integrands::{
+            AdjointDoubleLayerBoundaryIntegrand, BoundaryIntegrandScalarProduct,
+            BoundaryIntegrandSum, DoubleLayerBoundaryIntegrand,
+            HypersingularCurlCurlBoundaryIntegrand, HypersingularNormalNormalBoundaryIntegrand,
+            SingleLayerBoundaryIntegrand,
+        },
+        assembly::boundary::BoundaryAssembler,
+        assembly::kernels::KernelEvaluator,
+        function::SerialFunctionSpace,
+        traits::{
+            BoundaryAssembly, BoundaryIntegrand, FunctionSpace,
+            KernelEvaluator as KernelEvaluatorTrait,
+        },
+    };
+    use green_kernels::{helmholtz_3d::Helmholtz3dKernel, laplace_3d::Laplace3dKernel};
+    use ndelement::{ciarlet::CiarletElement, types::ReferenceCellType};
+    use ndgrid::SingleElementGrid;
+    use rlst::{c32, c64, rlst_array_from_slice_mut2, MatrixInverse, RlstScalar};
+    use std::ffi::c_void;
+    use std::slice::from_raw_parts_mut;
+
+    type LaplaceHypersingularBoundaryIntegrand<T> = HypersingularCurlCurlBoundaryIntegrand<T>;
+    type HelmholtzHypersingularBoundaryIntegrand<T> = BoundaryIntegrandSum<
+        T,
+        HypersingularCurlCurlBoundaryIntegrand<T>,
+        BoundaryIntegrandScalarProduct<T, HypersingularNormalNormalBoundaryIntegrand<T>>,
+    >;
+
+    #[derive(Debug, PartialEq, Clone, Copy)]
+    #[repr(u8)]
+    pub enum BoundaryOperator {
+        SingleLayer = 0,
+        DoubleLayer = 1,
+        AdjointDoubleLayer = 2,
+        Hypersingular = 3,
+        ElectricField = 4,
+        MagneticField = 5,
+    }
+
+    impl BoundaryOperator {
+        fn from(value: u8) -> Option<Self> {
+            match value {
+                0 => Some(BoundaryOperator::SingleLayer),
+                1 => Some(BoundaryOperator::DoubleLayer),
+                2 => Some(BoundaryOperator::AdjointDoubleLayer),
+                3 => Some(BoundaryOperator::Hypersingular),
+                4 => Some(BoundaryOperator::ElectricField),
+                5 => Some(BoundaryOperator::MagneticField),
+                _ => None,
+            }
+        }
+    }
+
+    #[derive(Debug, PartialEq, Clone, Copy)]
+    #[repr(u8)]
+    pub enum KernelType {
+        Laplace = 0,
+        Helmholtz = 1,
+    }
+
+    #[repr(C)]
+    pub struct BoundaryAssemblerWrapper {
+        pub assembler: *const c_void,
+        pub itype: BoundaryOperator,
+        pub ktype: KernelType,
+        pub dtype: DType,
+    }
+    impl Drop for BoundaryAssemblerWrapper {
+        fn drop(&mut self) {
+            let Self {
+                assembler,
+                itype,
+                ktype,
+                dtype,
+            } = self;
+            match ktype {
+                KernelType::Laplace => match itype {
+                    BoundaryOperator::SingleLayer => match dtype {
+                        DType::F32 => drop(unsafe {
+                            Box::from_raw(
+                                *assembler
+                                    as *mut BoundaryAssembler<
+                                        f32,
+                                        SingleLayerBoundaryIntegrand<f32>,
+                                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                                    >,
+                            )
+                        }),
+                        DType::F64 => drop(unsafe {
+                            Box::from_raw(
+                                *assembler
+                                    as *mut BoundaryAssembler<
+                                        f64,
+                                        SingleLayerBoundaryIntegrand<f64>,
+                                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                                    >,
+                            )
+                        }),
+                        _ => {
+                            panic!("Invalid data type");
+                        }
+                    },
+                    BoundaryOperator::DoubleLayer => match dtype {
+                        DType::F32 => drop(unsafe {
+                            Box::from_raw(
+                                *assembler
+                                    as *mut BoundaryAssembler<
+                                        f32,
+                                        DoubleLayerBoundaryIntegrand<f32>,
+                                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                                    >,
+                            )
+                        }),
+                        DType::F64 => drop(unsafe {
+                            Box::from_raw(
+                                *assembler
+                                    as *mut BoundaryAssembler<
+                                        f64,
+                                        DoubleLayerBoundaryIntegrand<f64>,
+                                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                                    >,
+                            )
+                        }),
+                        _ => {
+                            panic!("Invalid data type");
+                        }
+                    },
+                    BoundaryOperator::AdjointDoubleLayer => match dtype {
+                        DType::F32 => drop(unsafe {
+                            Box::from_raw(
+                                *assembler
+                                    as *mut BoundaryAssembler<
+                                        f32,
+                                        AdjointDoubleLayerBoundaryIntegrand<f32>,
+                                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                                    >,
+                            )
+                        }),
+                        DType::F64 => drop(unsafe {
+                            Box::from_raw(
+                                *assembler
+                                    as *mut BoundaryAssembler<
+                                        f64,
+                                        AdjointDoubleLayerBoundaryIntegrand<f64>,
+                                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                                    >,
+                            )
+                        }),
+                        _ => {
+                            panic!("Invalid data type");
+                        }
+                    },
+                    BoundaryOperator::Hypersingular => match dtype {
+                        DType::F32 => drop(unsafe {
+                            Box::from_raw(
+                                *assembler
+                                    as *mut BoundaryAssembler<
+                                        f32,
+                                        LaplaceHypersingularBoundaryIntegrand<f32>,
+                                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                                    >,
+                            )
+                        }),
+                        DType::F64 => drop(unsafe {
+                            Box::from_raw(
+                                *assembler
+                                    as *mut BoundaryAssembler<
+                                        f64,
+                                        LaplaceHypersingularBoundaryIntegrand<f64>,
+                                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                                    >,
+                            )
+                        }),
+                        _ => {
+                            panic!("Invalid data type");
+                        }
+                    },
+                    _ => {
+                        panic!("Invalid operator");
+                    }
+                },
+                KernelType::Helmholtz => match itype {
+                    BoundaryOperator::SingleLayer => match dtype {
+                        DType::C32 => drop(unsafe {
+                            Box::from_raw(
+                                *assembler
+                                    as *mut BoundaryAssembler<
+                                        c32,
+                                        SingleLayerBoundaryIntegrand<c32>,
+                                        KernelEvaluator<c32, Laplace3dKernel<c32>>,
+                                    >,
+                            )
+                        }),
+                        DType::C64 => drop(unsafe {
+                            Box::from_raw(
+                                *assembler
+                                    as *mut BoundaryAssembler<
+                                        c64,
+                                        SingleLayerBoundaryIntegrand<c64>,
+                                        KernelEvaluator<c64, Laplace3dKernel<c64>>,
+                                    >,
+                            )
+                        }),
+                        _ => {
+                            panic!("Invalid data type");
+                        }
+                    },
+                    BoundaryOperator::DoubleLayer => match dtype {
+                        DType::C32 => drop(unsafe {
+                            Box::from_raw(
+                                *assembler
+                                    as *mut BoundaryAssembler<
+                                        c32,
+                                        DoubleLayerBoundaryIntegrand<c32>,
+                                        KernelEvaluator<c32, Laplace3dKernel<c32>>,
+                                    >,
+                            )
+                        }),
+                        DType::C64 => drop(unsafe {
+                            Box::from_raw(
+                                *assembler
+                                    as *mut BoundaryAssembler<
+                                        c64,
+                                        DoubleLayerBoundaryIntegrand<c64>,
+                                        KernelEvaluator<c64, Laplace3dKernel<c64>>,
+                                    >,
+                            )
+                        }),
+                        _ => {
+                            panic!("Invalid data type");
+                        }
+                    },
+                    BoundaryOperator::AdjointDoubleLayer => match dtype {
+                        DType::C32 => drop(unsafe {
+                            Box::from_raw(
+                                *assembler
+                                    as *mut BoundaryAssembler<
+                                        c32,
+                                        AdjointDoubleLayerBoundaryIntegrand<c32>,
+                                        KernelEvaluator<c32, Laplace3dKernel<c32>>,
+                                    >,
+                            )
+                        }),
+                        DType::C64 => drop(unsafe {
+                            Box::from_raw(
+                                *assembler
+                                    as *mut BoundaryAssembler<
+                                        c64,
+                                        AdjointDoubleLayerBoundaryIntegrand<c64>,
+                                        KernelEvaluator<c64, Laplace3dKernel<c64>>,
+                                    >,
+                            )
+                        }),
+                        _ => {
+                            panic!("Invalid data type");
+                        }
+                    },
+                    BoundaryOperator::Hypersingular => match dtype {
+                        DType::C32 => drop(unsafe {
+                            Box::from_raw(
+                                *assembler
+                                    as *mut BoundaryAssembler<
+                                        c32,
+                                        HelmholtzHypersingularBoundaryIntegrand<c32>,
+                                        KernelEvaluator<c32, Laplace3dKernel<c32>>,
+                                    >,
+                            )
+                        }),
+                        DType::C64 => drop(unsafe {
+                            Box::from_raw(
+                                *assembler
+                                    as *mut BoundaryAssembler<
+                                        c64,
+                                        HelmholtzHypersingularBoundaryIntegrand<c64>,
+                                        KernelEvaluator<c64, Laplace3dKernel<c64>>,
+                                    >,
+                            )
+                        }),
+                        _ => {
+                            panic!("Invalid data type");
+                        }
+                    },
+                    _ => {
+                        panic!("Invalid operator");
+                    }
+                },
+            }
+        }
+    }
+
+    #[no_mangle]
+    pub unsafe extern "C" fn free_boundary_assembler(a: *mut BoundaryAssemblerWrapper) {
+        assert!(!a.is_null());
+        unsafe { drop(Box::from_raw(a)) }
+    }
+
+    #[repr(C)]
+    pub struct SparseMatrixWrapper {
+        pub matrix: *const c_void,
+        pub dtype: DType,
+    }
+    impl Drop for SparseMatrixWrapper {
+        fn drop(&mut self) {
+            let Self { matrix, dtype } = self;
+            match dtype {
+                DType::F32 => drop(unsafe { Box::from_raw(*matrix as *mut SparseMatrixData<f32>) }),
+                DType::F64 => drop(unsafe { Box::from_raw(*matrix as *mut SparseMatrixData<f64>) }),
+                DType::C32 => drop(unsafe { Box::from_raw(*matrix as *mut SparseMatrixData<c32>) }),
+                DType::C64 => drop(unsafe { Box::from_raw(*matrix as *mut SparseMatrixData<c64>) }),
+            }
+        }
+    }
+
+    #[no_mangle]
+    pub unsafe extern "C" fn free_sparse_matrix(m: *mut SparseMatrixWrapper) {
+        assert!(!m.is_null());
+        unsafe { drop(Box::from_raw(m)) }
+    }
+
+    pub(crate) unsafe fn extract_sparse_matrix<T: RlstScalar + MatrixInverse>(
+        matrix: *const SparseMatrixWrapper,
+    ) -> *const SparseMatrixData<T> {
+        (*matrix).matrix as *const SparseMatrixData<T>
+    }
+
+    #[no_mangle]
+    pub unsafe extern "C" fn sparse_data_size(matrix: *mut SparseMatrixWrapper) -> usize {
+        match (*matrix).dtype {
+            DType::F32 => (*extract_sparse_matrix::<f32>(matrix)).data.len(),
+            DType::F64 => (*extract_sparse_matrix::<f64>(matrix)).data.len(),
+            DType::C32 => (*extract_sparse_matrix::<c32>(matrix)).data.len(),
+            DType::C64 => (*extract_sparse_matrix::<c64>(matrix)).data.len(),
+        }
+    }
+
+    unsafe fn sparse_data_internal<T: RlstScalar + MatrixInverse>(
+        matrix: *mut SparseMatrixWrapper,
+        data: *mut c_void,
+    ) {
+        let data = data as *mut T;
+        for (i, j) in (*extract_sparse_matrix::<T>(matrix))
+            .data
+            .iter()
+            .enumerate()
+        {
+            *data.add(i) = *j;
+        }
+    }
+    #[no_mangle]
+    pub unsafe extern "C" fn sparse_data(matrix: *mut SparseMatrixWrapper, data: *mut c_void) {
+        match (*matrix).dtype {
+            DType::F32 => sparse_data_internal::<f32>(matrix, data),
+            DType::F64 => sparse_data_internal::<f64>(matrix, data),
+            DType::C32 => sparse_data_internal::<c32>(matrix, data),
+            DType::C64 => sparse_data_internal::<c64>(matrix, data),
+        }
+    }
+
+    unsafe fn sparse_rows_internal<T: RlstScalar + MatrixInverse>(
+        matrix: *mut SparseMatrixWrapper,
+        rows: *mut usize,
+    ) {
+        for (i, j) in (*extract_sparse_matrix::<T>(matrix))
+            .rows
+            .iter()
+            .enumerate()
+        {
+            *rows.add(i) = *j;
+        }
+    }
+    #[no_mangle]
+    pub unsafe extern "C" fn sparse_rows(matrix: *mut SparseMatrixWrapper, rows: *mut usize) {
+        match (*matrix).dtype {
+            DType::F32 => sparse_rows_internal::<f32>(matrix, rows),
+            DType::F64 => sparse_rows_internal::<f64>(matrix, rows),
+            DType::C32 => sparse_rows_internal::<c32>(matrix, rows),
+            DType::C64 => sparse_rows_internal::<c64>(matrix, rows),
+        }
+    }
+
+    unsafe fn sparse_cols_internal<T: RlstScalar + MatrixInverse>(
+        matrix: *mut SparseMatrixWrapper,
+        cols: *mut usize,
+    ) {
+        for (i, j) in (*extract_sparse_matrix::<T>(matrix))
+            .cols
+            .iter()
+            .enumerate()
+        {
+            *cols.add(i) = *j;
+        }
+    }
+    #[no_mangle]
+    pub unsafe extern "C" fn sparse_cols(matrix: *mut SparseMatrixWrapper, cols: *mut usize) {
+        match (*matrix).dtype {
+            DType::F32 => sparse_cols_internal::<f32>(matrix, cols),
+            DType::F64 => sparse_cols_internal::<f64>(matrix, cols),
+            DType::C32 => sparse_cols_internal::<c32>(matrix, cols),
+            DType::C64 => sparse_cols_internal::<c64>(matrix, cols),
+        }
+    }
+
+    #[no_mangle]
+    pub unsafe extern "C" fn sparse_shape0(matrix: *mut SparseMatrixWrapper) -> usize {
+        match (*matrix).dtype {
+            DType::F32 => (*extract_sparse_matrix::<f32>(matrix)).shape[0],
+            DType::F64 => (*extract_sparse_matrix::<f64>(matrix)).shape[0],
+            DType::C32 => (*extract_sparse_matrix::<c32>(matrix)).shape[0],
+            DType::C64 => (*extract_sparse_matrix::<c64>(matrix)).shape[0],
+        }
+    }
+
+    #[no_mangle]
+    pub unsafe extern "C" fn sparse_shape1(matrix: *mut SparseMatrixWrapper) -> usize {
+        match (*matrix).dtype {
+            DType::F32 => (*extract_sparse_matrix::<f32>(matrix)).shape[1],
+            DType::F64 => (*extract_sparse_matrix::<f64>(matrix)).shape[1],
+            DType::C32 => (*extract_sparse_matrix::<c32>(matrix)).shape[1],
+            DType::C64 => (*extract_sparse_matrix::<c64>(matrix)).shape[1],
+        }
+    }
+
+    pub(crate) unsafe fn extract_boundary_assembler<
+        T: RlstScalar + MatrixInverse,
+        Integrand: BoundaryIntegrand<T = T>,
+        Kernel: KernelEvaluatorTrait<T = T>,
+    >(
+        assembler: *const BoundaryAssemblerWrapper,
+    ) -> *const BoundaryAssembler<T, Integrand, Kernel> {
+        (*assembler).assembler as *const BoundaryAssembler<T, Integrand, Kernel>
+    }
+
+    pub(crate) unsafe fn extract_boundary_assembler_mut<
+        T: RlstScalar + MatrixInverse,
+        Integrand: BoundaryIntegrand<T = T>,
+        Kernel: KernelEvaluatorTrait<T = T>,
+    >(
+        assembler: *const BoundaryAssemblerWrapper,
+    ) -> *mut BoundaryAssembler<T, Integrand, Kernel> {
+        (*assembler).assembler as *mut BoundaryAssembler<T, Integrand, Kernel>
+    }
+
+    #[no_mangle]
+    pub unsafe extern "C" fn boundary_assembler_has_quadrature_degree(
+        assembler: *mut BoundaryAssemblerWrapper,
+        cell: u8,
+    ) -> bool {
+        let cell = ReferenceCellType::from(cell).unwrap();
+        match (*assembler).ktype {
+            KernelType::Laplace => match (*assembler).itype {
+                BoundaryOperator::SingleLayer => match (*assembler).dtype {
+                    DType::F32 => (*extract_boundary_assembler::<
+                        f32,
+                        SingleLayerBoundaryIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    DType::F64 => (*extract_boundary_assembler::<
+                        f64,
+                        SingleLayerBoundaryIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::F32 => (*extract_boundary_assembler::<
+                        f32,
+                        DoubleLayerBoundaryIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    DType::F64 => (*extract_boundary_assembler::<
+                        f64,
+                        DoubleLayerBoundaryIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::AdjointDoubleLayer => match (*assembler).dtype {
+                    DType::F32 => (*extract_boundary_assembler::<
+                        f32,
+                        AdjointDoubleLayerBoundaryIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    DType::F64 => (*extract_boundary_assembler::<
+                        f64,
+                        AdjointDoubleLayerBoundaryIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::Hypersingular => match (*assembler).dtype {
+                    DType::F32 => (*extract_boundary_assembler::<
+                        f32,
+                        LaplaceHypersingularBoundaryIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    DType::F64 => (*extract_boundary_assembler::<
+                        f64,
+                        LaplaceHypersingularBoundaryIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+            KernelType::Helmholtz => match (*assembler).itype {
+                BoundaryOperator::SingleLayer => match (*assembler).dtype {
+                    DType::C32 => (*extract_boundary_assembler::<
+                        c32,
+                        SingleLayerBoundaryIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    DType::C64 => (*extract_boundary_assembler::<
+                        c64,
+                        SingleLayerBoundaryIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::C32 => (*extract_boundary_assembler::<
+                        c32,
+                        DoubleLayerBoundaryIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    DType::C64 => (*extract_boundary_assembler::<
+                        c64,
+                        DoubleLayerBoundaryIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::AdjointDoubleLayer => match (*assembler).dtype {
+                    DType::C32 => (*extract_boundary_assembler::<
+                        c32,
+                        AdjointDoubleLayerBoundaryIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    DType::C64 => (*extract_boundary_assembler::<
+                        c64,
+                        AdjointDoubleLayerBoundaryIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::Hypersingular => match (*assembler).dtype {
+                    DType::C32 => (*extract_boundary_assembler::<
+                        c32,
+                        HelmholtzHypersingularBoundaryIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    DType::C64 => (*extract_boundary_assembler::<
+                        c64,
+                        HelmholtzHypersingularBoundaryIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+        }
+        .is_some()
+    }
+
+    #[no_mangle]
+    pub unsafe extern "C" fn boundary_assembler_set_quadrature_degree(
+        assembler: *mut BoundaryAssemblerWrapper,
+        cell: u8,
+        degree: usize,
+    ) {
+        let cell = ReferenceCellType::from(cell).unwrap();
+        match (*assembler).ktype {
+            KernelType::Laplace => match (*assembler).itype {
+                BoundaryOperator::SingleLayer => match (*assembler).dtype {
+                    DType::F32 => (*extract_boundary_assembler_mut::<
+                        f32,
+                        SingleLayerBoundaryIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .set_quadrature_degree(cell, degree),
+                    DType::F64 => (*extract_boundary_assembler_mut::<
+                        f64,
+                        SingleLayerBoundaryIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .set_quadrature_degree(cell, degree),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::F32 => (*extract_boundary_assembler_mut::<
+                        f32,
+                        DoubleLayerBoundaryIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .set_quadrature_degree(cell, degree),
+                    DType::F64 => (*extract_boundary_assembler_mut::<
+                        f64,
+                        DoubleLayerBoundaryIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .set_quadrature_degree(cell, degree),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::AdjointDoubleLayer => match (*assembler).dtype {
+                    DType::F32 => (*extract_boundary_assembler_mut::<
+                        f32,
+                        AdjointDoubleLayerBoundaryIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .set_quadrature_degree(cell, degree),
+                    DType::F64 => (*extract_boundary_assembler_mut::<
+                        f64,
+                        AdjointDoubleLayerBoundaryIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .set_quadrature_degree(cell, degree),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::Hypersingular => match (*assembler).dtype {
+                    DType::F32 => (*extract_boundary_assembler_mut::<
+                        f32,
+                        LaplaceHypersingularBoundaryIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .set_quadrature_degree(cell, degree),
+                    DType::F64 => (*extract_boundary_assembler_mut::<
+                        f64,
+                        LaplaceHypersingularBoundaryIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .set_quadrature_degree(cell, degree),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+            KernelType::Helmholtz => match (*assembler).itype {
+                BoundaryOperator::SingleLayer => match (*assembler).dtype {
+                    DType::C32 => (*extract_boundary_assembler_mut::<
+                        c32,
+                        SingleLayerBoundaryIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .set_quadrature_degree(cell, degree),
+                    DType::C64 => (*extract_boundary_assembler_mut::<
+                        c64,
+                        SingleLayerBoundaryIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .set_quadrature_degree(cell, degree),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::C32 => (*extract_boundary_assembler_mut::<
+                        c32,
+                        DoubleLayerBoundaryIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .set_quadrature_degree(cell, degree),
+                    DType::C64 => (*extract_boundary_assembler_mut::<
+                        c64,
+                        DoubleLayerBoundaryIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .set_quadrature_degree(cell, degree),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::AdjointDoubleLayer => match (*assembler).dtype {
+                    DType::C32 => (*extract_boundary_assembler_mut::<
+                        c32,
+                        AdjointDoubleLayerBoundaryIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .set_quadrature_degree(cell, degree),
+                    DType::C64 => (*extract_boundary_assembler_mut::<
+                        c64,
+                        AdjointDoubleLayerBoundaryIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .set_quadrature_degree(cell, degree),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::Hypersingular => match (*assembler).dtype {
+                    DType::C32 => (*extract_boundary_assembler_mut::<
+                        c32,
+                        HelmholtzHypersingularBoundaryIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .set_quadrature_degree(cell, degree),
+                    DType::C64 => (*extract_boundary_assembler_mut::<
+                        c64,
+                        HelmholtzHypersingularBoundaryIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .set_quadrature_degree(cell, degree),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+        }
+    }
+
+    #[no_mangle]
+    pub unsafe extern "C" fn boundary_assembler_quadrature_degree(
+        assembler: *mut BoundaryAssemblerWrapper,
+        cell: u8,
+    ) -> usize {
+        let cell = ReferenceCellType::from(cell).unwrap();
+        match (*assembler).ktype {
+            KernelType::Laplace => match (*assembler).itype {
+                BoundaryOperator::SingleLayer => match (*assembler).dtype {
+                    DType::F32 => (*extract_boundary_assembler::<
+                        f32,
+                        SingleLayerBoundaryIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    DType::F64 => (*extract_boundary_assembler::<
+                        f64,
+                        SingleLayerBoundaryIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::F32 => (*extract_boundary_assembler::<
+                        f32,
+                        DoubleLayerBoundaryIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    DType::F64 => (*extract_boundary_assembler::<
+                        f64,
+                        DoubleLayerBoundaryIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::AdjointDoubleLayer => match (*assembler).dtype {
+                    DType::F32 => (*extract_boundary_assembler::<
+                        f32,
+                        AdjointDoubleLayerBoundaryIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    DType::F64 => (*extract_boundary_assembler::<
+                        f64,
+                        AdjointDoubleLayerBoundaryIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::Hypersingular => match (*assembler).dtype {
+                    DType::F32 => (*extract_boundary_assembler::<
+                        f32,
+                        LaplaceHypersingularBoundaryIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    DType::F64 => (*extract_boundary_assembler::<
+                        f64,
+                        LaplaceHypersingularBoundaryIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+            KernelType::Helmholtz => match (*assembler).itype {
+                BoundaryOperator::SingleLayer => match (*assembler).dtype {
+                    DType::C32 => (*extract_boundary_assembler::<
+                        c32,
+                        SingleLayerBoundaryIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    DType::C64 => (*extract_boundary_assembler::<
+                        c64,
+                        SingleLayerBoundaryIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::C32 => (*extract_boundary_assembler::<
+                        c32,
+                        DoubleLayerBoundaryIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    DType::C64 => (*extract_boundary_assembler::<
+                        c64,
+                        DoubleLayerBoundaryIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::AdjointDoubleLayer => match (*assembler).dtype {
+                    DType::C32 => (*extract_boundary_assembler::<
+                        c32,
+                        AdjointDoubleLayerBoundaryIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    DType::C64 => (*extract_boundary_assembler::<
+                        c64,
+                        AdjointDoubleLayerBoundaryIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::Hypersingular => match (*assembler).dtype {
+                    DType::C32 => (*extract_boundary_assembler::<
+                        c32,
+                        HelmholtzHypersingularBoundaryIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    DType::C64 => (*extract_boundary_assembler::<
+                        c64,
+                        HelmholtzHypersingularBoundaryIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+        }
+        .unwrap()
+    }
+
+    #[no_mangle]
+    pub unsafe extern "C" fn boundary_assembler_has_singular_quadrature_degree(
+        assembler: *mut BoundaryAssemblerWrapper,
+        cell0: u8,
+        cell1: u8,
+    ) -> bool {
+        let cells = (
+            ReferenceCellType::from(cell0).unwrap(),
+            ReferenceCellType::from(cell1).unwrap(),
+        );
+        match (*assembler).ktype {
+            KernelType::Laplace => match (*assembler).itype {
+                BoundaryOperator::SingleLayer => match (*assembler).dtype {
+                    DType::F32 => (*extract_boundary_assembler::<
+                        f32,
+                        SingleLayerBoundaryIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .singular_quadrature_degree(cells),
+                    DType::F64 => (*extract_boundary_assembler::<
+                        f64,
+                        SingleLayerBoundaryIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .singular_quadrature_degree(cells),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::F32 => (*extract_boundary_assembler::<
+                        f32,
+                        DoubleLayerBoundaryIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .singular_quadrature_degree(cells),
+                    DType::F64 => (*extract_boundary_assembler::<
+                        f64,
+                        DoubleLayerBoundaryIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .singular_quadrature_degree(cells),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::AdjointDoubleLayer => match (*assembler).dtype {
+                    DType::F32 => (*extract_boundary_assembler::<
+                        f32,
+                        AdjointDoubleLayerBoundaryIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .singular_quadrature_degree(cells),
+                    DType::F64 => (*extract_boundary_assembler::<
+                        f64,
+                        AdjointDoubleLayerBoundaryIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .singular_quadrature_degree(cells),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::Hypersingular => match (*assembler).dtype {
+                    DType::F32 => (*extract_boundary_assembler::<
+                        f32,
+                        LaplaceHypersingularBoundaryIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .singular_quadrature_degree(cells),
+                    DType::F64 => (*extract_boundary_assembler::<
+                        f64,
+                        LaplaceHypersingularBoundaryIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .singular_quadrature_degree(cells),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+            KernelType::Helmholtz => match (*assembler).itype {
+                BoundaryOperator::SingleLayer => match (*assembler).dtype {
+                    DType::C32 => (*extract_boundary_assembler::<
+                        c32,
+                        SingleLayerBoundaryIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .singular_quadrature_degree(cells),
+                    DType::C64 => (*extract_boundary_assembler::<
+                        c64,
+                        SingleLayerBoundaryIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .singular_quadrature_degree(cells),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::C32 => (*extract_boundary_assembler::<
+                        c32,
+                        DoubleLayerBoundaryIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .singular_quadrature_degree(cells),
+                    DType::C64 => (*extract_boundary_assembler::<
+                        c64,
+                        DoubleLayerBoundaryIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .singular_quadrature_degree(cells),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::AdjointDoubleLayer => match (*assembler).dtype {
+                    DType::C32 => (*extract_boundary_assembler::<
+                        c32,
+                        AdjointDoubleLayerBoundaryIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .singular_quadrature_degree(cells),
+                    DType::C64 => (*extract_boundary_assembler::<
+                        c64,
+                        AdjointDoubleLayerBoundaryIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .singular_quadrature_degree(cells),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::Hypersingular => match (*assembler).dtype {
+                    DType::C32 => (*extract_boundary_assembler::<
+                        c32,
+                        HelmholtzHypersingularBoundaryIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .singular_quadrature_degree(cells),
+                    DType::C64 => (*extract_boundary_assembler::<
+                        c64,
+                        HelmholtzHypersingularBoundaryIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .singular_quadrature_degree(cells),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+        }
+        .is_some()
+    }
+
+    #[no_mangle]
+    pub unsafe extern "C" fn boundary_assembler_set_singular_quadrature_degree(
+        assembler: *mut BoundaryAssemblerWrapper,
+        cell0: u8,
+        cell1: u8,
+        degree: usize,
+    ) {
+        let cells = (
+            ReferenceCellType::from(cell0).unwrap(),
+            ReferenceCellType::from(cell1).unwrap(),
+        );
+        match (*assembler).ktype {
+            KernelType::Laplace => match (*assembler).itype {
+                BoundaryOperator::SingleLayer => match (*assembler).dtype {
+                    DType::F32 => (*extract_boundary_assembler_mut::<
+                        f32,
+                        SingleLayerBoundaryIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .set_singular_quadrature_degree(cells, degree),
+                    DType::F64 => (*extract_boundary_assembler_mut::<
+                        f64,
+                        SingleLayerBoundaryIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .set_singular_quadrature_degree(cells, degree),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::F32 => (*extract_boundary_assembler_mut::<
+                        f32,
+                        DoubleLayerBoundaryIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .set_singular_quadrature_degree(cells, degree),
+                    DType::F64 => (*extract_boundary_assembler_mut::<
+                        f64,
+                        DoubleLayerBoundaryIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .set_singular_quadrature_degree(cells, degree),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::AdjointDoubleLayer => match (*assembler).dtype {
+                    DType::F32 => (*extract_boundary_assembler_mut::<
+                        f32,
+                        AdjointDoubleLayerBoundaryIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .set_singular_quadrature_degree(cells, degree),
+                    DType::F64 => (*extract_boundary_assembler_mut::<
+                        f64,
+                        AdjointDoubleLayerBoundaryIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .set_singular_quadrature_degree(cells, degree),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::Hypersingular => match (*assembler).dtype {
+                    DType::F32 => (*extract_boundary_assembler_mut::<
+                        f32,
+                        LaplaceHypersingularBoundaryIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .set_singular_quadrature_degree(cells, degree),
+                    DType::F64 => (*extract_boundary_assembler_mut::<
+                        f64,
+                        LaplaceHypersingularBoundaryIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .set_singular_quadrature_degree(cells, degree),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+            KernelType::Helmholtz => match (*assembler).itype {
+                BoundaryOperator::SingleLayer => match (*assembler).dtype {
+                    DType::C32 => (*extract_boundary_assembler_mut::<
+                        c32,
+                        SingleLayerBoundaryIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .set_singular_quadrature_degree(cells, degree),
+                    DType::C64 => (*extract_boundary_assembler_mut::<
+                        c64,
+                        SingleLayerBoundaryIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .set_singular_quadrature_degree(cells, degree),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::C32 => (*extract_boundary_assembler_mut::<
+                        c32,
+                        DoubleLayerBoundaryIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .set_singular_quadrature_degree(cells, degree),
+                    DType::C64 => (*extract_boundary_assembler_mut::<
+                        c64,
+                        DoubleLayerBoundaryIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .set_singular_quadrature_degree(cells, degree),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::AdjointDoubleLayer => match (*assembler).dtype {
+                    DType::C32 => (*extract_boundary_assembler_mut::<
+                        c32,
+                        AdjointDoubleLayerBoundaryIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .set_singular_quadrature_degree(cells, degree),
+                    DType::C64 => (*extract_boundary_assembler_mut::<
+                        c64,
+                        AdjointDoubleLayerBoundaryIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .set_singular_quadrature_degree(cells, degree),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::Hypersingular => match (*assembler).dtype {
+                    DType::C32 => (*extract_boundary_assembler_mut::<
+                        c32,
+                        HelmholtzHypersingularBoundaryIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .set_singular_quadrature_degree(cells, degree),
+                    DType::C64 => (*extract_boundary_assembler_mut::<
+                        c64,
+                        HelmholtzHypersingularBoundaryIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .set_singular_quadrature_degree(cells, degree),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+        }
+    }
+
+    #[no_mangle]
+    pub unsafe extern "C" fn boundary_assembler_singular_quadrature_degree(
+        assembler: *mut BoundaryAssemblerWrapper,
+        cell0: u8,
+        cell1: u8,
+    ) -> usize {
+        let cells = (
+            ReferenceCellType::from(cell0).unwrap(),
+            ReferenceCellType::from(cell1).unwrap(),
+        );
+        match (*assembler).ktype {
+            KernelType::Laplace => match (*assembler).itype {
+                BoundaryOperator::SingleLayer => match (*assembler).dtype {
+                    DType::F32 => (*extract_boundary_assembler::<
+                        f32,
+                        SingleLayerBoundaryIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .singular_quadrature_degree(cells),
+                    DType::F64 => (*extract_boundary_assembler::<
+                        f64,
+                        SingleLayerBoundaryIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .singular_quadrature_degree(cells),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::F32 => (*extract_boundary_assembler::<
+                        f32,
+                        DoubleLayerBoundaryIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .singular_quadrature_degree(cells),
+                    DType::F64 => (*extract_boundary_assembler::<
+                        f64,
+                        DoubleLayerBoundaryIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .singular_quadrature_degree(cells),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::AdjointDoubleLayer => match (*assembler).dtype {
+                    DType::F32 => (*extract_boundary_assembler::<
+                        f32,
+                        AdjointDoubleLayerBoundaryIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .singular_quadrature_degree(cells),
+                    DType::F64 => (*extract_boundary_assembler::<
+                        f64,
+                        AdjointDoubleLayerBoundaryIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .singular_quadrature_degree(cells),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::Hypersingular => match (*assembler).dtype {
+                    DType::F32 => (*extract_boundary_assembler::<
+                        f32,
+                        LaplaceHypersingularBoundaryIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .singular_quadrature_degree(cells),
+                    DType::F64 => (*extract_boundary_assembler::<
+                        f64,
+                        LaplaceHypersingularBoundaryIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .singular_quadrature_degree(cells),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+            KernelType::Helmholtz => match (*assembler).itype {
+                BoundaryOperator::SingleLayer => match (*assembler).dtype {
+                    DType::C32 => (*extract_boundary_assembler::<
+                        c32,
+                        SingleLayerBoundaryIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .singular_quadrature_degree(cells),
+                    DType::C64 => (*extract_boundary_assembler::<
+                        c64,
+                        SingleLayerBoundaryIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .singular_quadrature_degree(cells),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::C32 => (*extract_boundary_assembler::<
+                        c32,
+                        DoubleLayerBoundaryIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .singular_quadrature_degree(cells),
+                    DType::C64 => (*extract_boundary_assembler::<
+                        c64,
+                        DoubleLayerBoundaryIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .singular_quadrature_degree(cells),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::AdjointDoubleLayer => match (*assembler).dtype {
+                    DType::C32 => (*extract_boundary_assembler::<
+                        c32,
+                        AdjointDoubleLayerBoundaryIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .singular_quadrature_degree(cells),
+                    DType::C64 => (*extract_boundary_assembler::<
+                        c64,
+                        AdjointDoubleLayerBoundaryIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .singular_quadrature_degree(cells),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::Hypersingular => match (*assembler).dtype {
+                    DType::C32 => (*extract_boundary_assembler::<
+                        c32,
+                        HelmholtzHypersingularBoundaryIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .singular_quadrature_degree(cells),
+                    DType::C64 => (*extract_boundary_assembler::<
+                        c64,
+                        HelmholtzHypersingularBoundaryIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .singular_quadrature_degree(cells),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+        }
+        .unwrap()
+    }
+
+    #[no_mangle]
+    pub unsafe extern "C" fn boundary_assembler_set_batch_size(
+        assembler: *mut BoundaryAssemblerWrapper,
+        batch_size: usize,
+    ) {
+        match (*assembler).ktype {
+            KernelType::Laplace => match (*assembler).itype {
+                BoundaryOperator::SingleLayer => match (*assembler).dtype {
+                    DType::F32 => (*extract_boundary_assembler_mut::<
+                        f32,
+                        SingleLayerBoundaryIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .set_batch_size(batch_size),
+                    DType::F64 => (*extract_boundary_assembler_mut::<
+                        f64,
+                        SingleLayerBoundaryIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .set_batch_size(batch_size),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::F32 => (*extract_boundary_assembler_mut::<
+                        f32,
+                        DoubleLayerBoundaryIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .set_batch_size(batch_size),
+                    DType::F64 => (*extract_boundary_assembler_mut::<
+                        f64,
+                        DoubleLayerBoundaryIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .set_batch_size(batch_size),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::AdjointDoubleLayer => match (*assembler).dtype {
+                    DType::F32 => (*extract_boundary_assembler_mut::<
+                        f32,
+                        AdjointDoubleLayerBoundaryIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .set_batch_size(batch_size),
+                    DType::F64 => (*extract_boundary_assembler_mut::<
+                        f64,
+                        AdjointDoubleLayerBoundaryIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .set_batch_size(batch_size),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::Hypersingular => match (*assembler).dtype {
+                    DType::F32 => (*extract_boundary_assembler_mut::<
+                        f32,
+                        LaplaceHypersingularBoundaryIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .set_batch_size(batch_size),
+                    DType::F64 => (*extract_boundary_assembler_mut::<
+                        f64,
+                        LaplaceHypersingularBoundaryIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .set_batch_size(batch_size),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+            KernelType::Helmholtz => match (*assembler).itype {
+                BoundaryOperator::SingleLayer => match (*assembler).dtype {
+                    DType::C32 => (*extract_boundary_assembler_mut::<
+                        c32,
+                        SingleLayerBoundaryIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .set_batch_size(batch_size),
+                    DType::C64 => (*extract_boundary_assembler_mut::<
+                        c64,
+                        SingleLayerBoundaryIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .set_batch_size(batch_size),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::C32 => (*extract_boundary_assembler_mut::<
+                        c32,
+                        DoubleLayerBoundaryIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .set_batch_size(batch_size),
+                    DType::C64 => (*extract_boundary_assembler_mut::<
+                        c64,
+                        DoubleLayerBoundaryIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .set_batch_size(batch_size),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::AdjointDoubleLayer => match (*assembler).dtype {
+                    DType::C32 => (*extract_boundary_assembler_mut::<
+                        c32,
+                        AdjointDoubleLayerBoundaryIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .set_batch_size(batch_size),
+                    DType::C64 => (*extract_boundary_assembler_mut::<
+                        c64,
+                        AdjointDoubleLayerBoundaryIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .set_batch_size(batch_size),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::Hypersingular => match (*assembler).dtype {
+                    DType::C32 => (*extract_boundary_assembler_mut::<
+                        c32,
+                        HelmholtzHypersingularBoundaryIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .set_batch_size(batch_size),
+                    DType::C64 => (*extract_boundary_assembler_mut::<
+                        c64,
+                        HelmholtzHypersingularBoundaryIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .set_batch_size(batch_size),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+        }
+    }
+
+    #[no_mangle]
+    pub unsafe extern "C" fn boundary_assembler_batch_size(
+        assembler: *mut BoundaryAssemblerWrapper,
+    ) -> usize {
+        match (*assembler).ktype {
+            KernelType::Laplace => match (*assembler).itype {
+                BoundaryOperator::SingleLayer => match (*assembler).dtype {
+                    DType::F32 => (*extract_boundary_assembler::<
+                        f32,
+                        SingleLayerBoundaryIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .batch_size(),
+                    DType::F64 => (*extract_boundary_assembler::<
+                        f64,
+                        SingleLayerBoundaryIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .batch_size(),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::F32 => (*extract_boundary_assembler::<
+                        f32,
+                        DoubleLayerBoundaryIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .batch_size(),
+                    DType::F64 => (*extract_boundary_assembler::<
+                        f64,
+                        DoubleLayerBoundaryIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .batch_size(),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::AdjointDoubleLayer => match (*assembler).dtype {
+                    DType::F32 => (*extract_boundary_assembler::<
+                        f32,
+                        AdjointDoubleLayerBoundaryIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .batch_size(),
+                    DType::F64 => (*extract_boundary_assembler::<
+                        f64,
+                        AdjointDoubleLayerBoundaryIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .batch_size(),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::Hypersingular => match (*assembler).dtype {
+                    DType::F32 => (*extract_boundary_assembler::<
+                        f32,
+                        LaplaceHypersingularBoundaryIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .batch_size(),
+                    DType::F64 => (*extract_boundary_assembler::<
+                        f64,
+                        LaplaceHypersingularBoundaryIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .batch_size(),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+            KernelType::Helmholtz => match (*assembler).itype {
+                BoundaryOperator::SingleLayer => match (*assembler).dtype {
+                    DType::C32 => (*extract_boundary_assembler::<
+                        c32,
+                        SingleLayerBoundaryIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .batch_size(),
+                    DType::C64 => (*extract_boundary_assembler::<
+                        c64,
+                        SingleLayerBoundaryIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .batch_size(),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::C32 => (*extract_boundary_assembler::<
+                        c32,
+                        DoubleLayerBoundaryIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .batch_size(),
+                    DType::C64 => (*extract_boundary_assembler::<
+                        c64,
+                        DoubleLayerBoundaryIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .batch_size(),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::AdjointDoubleLayer => match (*assembler).dtype {
+                    DType::C32 => (*extract_boundary_assembler::<
+                        c32,
+                        AdjointDoubleLayerBoundaryIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .batch_size(),
+                    DType::C64 => (*extract_boundary_assembler::<
+                        c64,
+                        AdjointDoubleLayerBoundaryIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .batch_size(),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::Hypersingular => match (*assembler).dtype {
+                    DType::C32 => (*extract_boundary_assembler::<
+                        c32,
+                        HelmholtzHypersingularBoundaryIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .batch_size(),
+                    DType::C64 => (*extract_boundary_assembler::<
+                        c64,
+                        HelmholtzHypersingularBoundaryIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .batch_size(),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+        }
+    }
+
+    unsafe fn boundary_assembler_assemble_into_dense_internal_real<
+        T: RlstScalar<Real = T> + MatrixInverse,
+        Integrand: BoundaryIntegrand<T = T>,
+        Kernel: KernelEvaluatorTrait<T = T>,
+        Space: FunctionSpace<T = T> + Sync,
+    >(
+        assembler: *mut BoundaryAssemblerWrapper,
+        output: *mut c_void,
+        trial_space: *const FunctionSpaceWrapper,
+        test_space: *const FunctionSpaceWrapper,
+    ) {
+        let trial_dim = (*extract_space::<Space>(trial_space)).global_size();
+        let test_dim = (*extract_space::<Space>(test_space)).global_size();
+        let mut output = rlst_array_from_slice_mut2!(
+            from_raw_parts_mut(output as *mut T, trial_dim * test_dim),
+            [test_dim, trial_dim]
+        );
+
+        (*extract_boundary_assembler::<T, Integrand, Kernel>(assembler)).assemble_into_dense(
+            &mut output,
+            &*extract_space::<Space>(trial_space),
+            &*extract_space::<Space>(test_space),
+        )
+    }
+    unsafe fn boundary_assembler_assemble_into_dense_internal_complex<
+        T: RlstScalar<Complex = T> + MatrixInverse,
+        Integrand: BoundaryIntegrand<T = T>,
+        Kernel: KernelEvaluatorTrait<T = T>,
+        Space: FunctionSpace<T = T> + Sync,
+    >(
+        assembler: *mut BoundaryAssemblerWrapper,
+        output: *mut c_void,
+        trial_space: *const FunctionSpaceWrapper,
+        test_space: *const FunctionSpaceWrapper,
+    ) {
+        let trial_dim = (*extract_space::<Space>(trial_space)).global_size();
+        let test_dim = (*extract_space::<Space>(test_space)).global_size();
+        let mut output = rlst_array_from_slice_mut2!(
+            from_raw_parts_mut(output as *mut T, trial_dim * test_dim),
+            [test_dim, trial_dim]
+        );
+
+        (*extract_boundary_assembler::<T, Integrand, Kernel>(assembler)).assemble_into_dense(
+            &mut output,
+            &*extract_space::<Space>(trial_space),
+            &*extract_space::<Space>(test_space),
+        )
+    }
+    #[no_mangle]
+    pub unsafe extern "C" fn boundary_assembler_assemble_into_dense(
+        assembler: *mut BoundaryAssemblerWrapper,
+        output: *mut c_void,
+        trial_space: *const FunctionSpaceWrapper,
+        test_space: *const FunctionSpaceWrapper,
+    ) {
+        match (*assembler).dtype {
+            DType::F32 => {
+                assert_eq!((*trial_space).dtype, DType::F32);
+                assert_eq!((*test_space).dtype, DType::F32);
+            }
+            DType::F64 => {
+                assert_eq!((*trial_space).dtype, DType::F64);
+                assert_eq!((*test_space).dtype, DType::F64);
+            }
+            DType::C32 => {
+                assert_eq!((*trial_space).dtype, DType::F32);
+                assert_eq!((*test_space).dtype, DType::F32);
+            }
+            DType::C64 => {
+                assert_eq!((*trial_space).dtype, DType::F64);
+                assert_eq!((*test_space).dtype, DType::F64);
+            }
+        }
+        assert_eq!((*test_space).stype, (*trial_space).stype);
+        assert_eq!((*test_space).gtype, (*trial_space).gtype);
+        match (*assembler).ktype {
+            KernelType::Laplace => match (*assembler).itype {
+                BoundaryOperator::SingleLayer => match (*assembler).dtype {
+                    DType::F32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_into_dense_internal_real::<
+                                    f32,
+                                    SingleLayerBoundaryIntegrand<f32>,
+                                    KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                                    SerialFunctionSpace<
+                                        f32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    DType::F64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_into_dense_internal_real::<
+                                    f64,
+                                    SingleLayerBoundaryIntegrand<f64>,
+                                    KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                                    SerialFunctionSpace<
+                                        f64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::F32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_into_dense_internal_real::<
+                                    f32,
+                                    DoubleLayerBoundaryIntegrand<f32>,
+                                    KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                                    SerialFunctionSpace<
+                                        f32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    DType::F64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_into_dense_internal_real::<
+                                    f64,
+                                    DoubleLayerBoundaryIntegrand<f64>,
+                                    KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                                    SerialFunctionSpace<
+                                        f64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::AdjointDoubleLayer => match (*assembler).dtype {
+                    DType::F32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_into_dense_internal_real::<
+                                    f32,
+                                    AdjointDoubleLayerBoundaryIntegrand<f32>,
+                                    KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                                    SerialFunctionSpace<
+                                        f32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    DType::F64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_into_dense_internal_real::<
+                                    f64,
+                                    AdjointDoubleLayerBoundaryIntegrand<f64>,
+                                    KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                                    SerialFunctionSpace<
+                                        f64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::Hypersingular => match (*assembler).dtype {
+                    DType::F32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_into_dense_internal_real::<
+                                    f32,
+                                    LaplaceHypersingularBoundaryIntegrand<f32>,
+                                    KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                                    SerialFunctionSpace<
+                                        f32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    DType::F64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_into_dense_internal_real::<
+                                    f64,
+                                    LaplaceHypersingularBoundaryIntegrand<f64>,
+                                    KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                                    SerialFunctionSpace<
+                                        f64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+            KernelType::Helmholtz => match (*assembler).itype {
+                BoundaryOperator::SingleLayer => match (*assembler).dtype {
+                    DType::C32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_into_dense_internal_complex::<
+                                    c32,
+                                    SingleLayerBoundaryIntegrand<c32>,
+                                    KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                                    SerialFunctionSpace<
+                                        c32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    DType::C64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_into_dense_internal_complex::<
+                                    c64,
+                                    SingleLayerBoundaryIntegrand<c64>,
+                                    KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                                    SerialFunctionSpace<
+                                        c64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::C32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_into_dense_internal_complex::<
+                                    c32,
+                                    DoubleLayerBoundaryIntegrand<c32>,
+                                    KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                                    SerialFunctionSpace<
+                                        c32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    DType::C64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_into_dense_internal_complex::<
+                                    c64,
+                                    DoubleLayerBoundaryIntegrand<c64>,
+                                    KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                                    SerialFunctionSpace<
+                                        c64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::AdjointDoubleLayer => match (*assembler).dtype {
+                    DType::C32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_into_dense_internal_complex::<
+                                    c32,
+                                    AdjointDoubleLayerBoundaryIntegrand<c32>,
+                                    KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                                    SerialFunctionSpace<
+                                        c32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    DType::C64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_into_dense_internal_complex::<
+                                    c64,
+                                    AdjointDoubleLayerBoundaryIntegrand<c64>,
+                                    KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                                    SerialFunctionSpace<
+                                        c64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::Hypersingular => match (*assembler).dtype {
+                    DType::C32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_into_dense_internal_complex::<
+                                    c32,
+                                    HelmholtzHypersingularBoundaryIntegrand<c32>,
+                                    KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                                    SerialFunctionSpace<
+                                        c32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    DType::C64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_into_dense_internal_complex::<
+                                    c64,
+                                    HelmholtzHypersingularBoundaryIntegrand<c64>,
+                                    KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                                    SerialFunctionSpace<
+                                        c64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+        }
+    }
+
+    unsafe fn boundary_assembler_assemble_singular_into_dense_internal_real<
+        T: RlstScalar<Real = T> + MatrixInverse,
+        Integrand: BoundaryIntegrand<T = T>,
+        Kernel: KernelEvaluatorTrait<T = T>,
+        Space: FunctionSpace<T = T> + Sync,
+    >(
+        assembler: *mut BoundaryAssemblerWrapper,
+        output: *mut c_void,
+        trial_space: *const FunctionSpaceWrapper,
+        test_space: *const FunctionSpaceWrapper,
+    ) {
+        let trial_dim = (*extract_space::<Space>(trial_space)).global_size();
+        let test_dim = (*extract_space::<Space>(test_space)).global_size();
+        let mut output = rlst_array_from_slice_mut2!(
+            from_raw_parts_mut(output as *mut T, trial_dim * test_dim),
+            [test_dim, trial_dim]
+        );
+
+        (*extract_boundary_assembler::<T, Integrand, Kernel>(assembler))
+            .assemble_singular_into_dense(
+                &mut output,
+                &*extract_space::<Space>(trial_space),
+                &*extract_space::<Space>(test_space),
+            )
+    }
+    unsafe fn boundary_assembler_assemble_singular_into_dense_internal_complex<
+        T: RlstScalar<Complex = T> + MatrixInverse,
+        Integrand: BoundaryIntegrand<T = T>,
+        Kernel: KernelEvaluatorTrait<T = T>,
+        Space: FunctionSpace<T = T> + Sync,
+    >(
+        assembler: *mut BoundaryAssemblerWrapper,
+        output: *mut c_void,
+        trial_space: *const FunctionSpaceWrapper,
+        test_space: *const FunctionSpaceWrapper,
+    ) {
+        let trial_dim = (*extract_space::<Space>(trial_space)).global_size();
+        let test_dim = (*extract_space::<Space>(test_space)).global_size();
+        let mut output = rlst_array_from_slice_mut2!(
+            from_raw_parts_mut(output as *mut T, trial_dim * test_dim),
+            [test_dim, trial_dim]
+        );
+
+        (*extract_boundary_assembler::<T, Integrand, Kernel>(assembler))
+            .assemble_singular_into_dense(
+                &mut output,
+                &*extract_space::<Space>(trial_space),
+                &*extract_space::<Space>(test_space),
+            )
+    }
+    #[no_mangle]
+    pub unsafe extern "C" fn boundary_assembler_assemble_singular_into_dense(
+        assembler: *mut BoundaryAssemblerWrapper,
+        output: *mut c_void,
+        trial_space: *const FunctionSpaceWrapper,
+        test_space: *const FunctionSpaceWrapper,
+    ) {
+        match (*assembler).dtype {
+            DType::F32 => {
+                assert_eq!((*trial_space).dtype, DType::F32);
+                assert_eq!((*test_space).dtype, DType::F32);
+            }
+            DType::F64 => {
+                assert_eq!((*trial_space).dtype, DType::F64);
+                assert_eq!((*test_space).dtype, DType::F64);
+            }
+            DType::C32 => {
+                assert_eq!((*trial_space).dtype, DType::F32);
+                assert_eq!((*test_space).dtype, DType::F32);
+            }
+            DType::C64 => {
+                assert_eq!((*trial_space).dtype, DType::F64);
+                assert_eq!((*test_space).dtype, DType::F64);
+            }
+        }
+        assert_eq!((*test_space).stype, (*trial_space).stype);
+        assert_eq!((*test_space).gtype, (*trial_space).gtype);
+        match (*assembler).ktype {
+            KernelType::Laplace => match (*assembler).itype {
+                BoundaryOperator::SingleLayer => match (*assembler).dtype {
+                    DType::F32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_into_dense_internal_real::<
+                                    f32,
+                                    SingleLayerBoundaryIntegrand<f32>,
+                                    KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                                    SerialFunctionSpace<
+                                        f32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    DType::F64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_into_dense_internal_real::<
+                                    f64,
+                                    SingleLayerBoundaryIntegrand<f64>,
+                                    KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                                    SerialFunctionSpace<
+                                        f64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::F32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_into_dense_internal_real::<
+                                    f32,
+                                    DoubleLayerBoundaryIntegrand<f32>,
+                                    KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                                    SerialFunctionSpace<
+                                        f32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    DType::F64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_into_dense_internal_real::<
+                                    f64,
+                                    DoubleLayerBoundaryIntegrand<f64>,
+                                    KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                                    SerialFunctionSpace<
+                                        f64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::AdjointDoubleLayer => match (*assembler).dtype {
+                    DType::F32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_into_dense_internal_real::<
+                                    f32,
+                                    AdjointDoubleLayerBoundaryIntegrand<f32>,
+                                    KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                                    SerialFunctionSpace<
+                                        f32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    DType::F64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_into_dense_internal_real::<
+                                    f64,
+                                    AdjointDoubleLayerBoundaryIntegrand<f64>,
+                                    KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                                    SerialFunctionSpace<
+                                        f64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::Hypersingular => match (*assembler).dtype {
+                    DType::F32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_into_dense_internal_real::<
+                                    f32,
+                                    LaplaceHypersingularBoundaryIntegrand<f32>,
+                                    KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                                    SerialFunctionSpace<
+                                        f32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    DType::F64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_into_dense_internal_real::<
+                                    f64,
+                                    LaplaceHypersingularBoundaryIntegrand<f64>,
+                                    KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                                    SerialFunctionSpace<
+                                        f64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+            KernelType::Helmholtz => match (*assembler).itype {
+                BoundaryOperator::SingleLayer => match (*assembler).dtype {
+                    DType::C32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_into_dense_internal_complex::<
+                                    c32,
+                                    SingleLayerBoundaryIntegrand<c32>,
+                                    KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                                    SerialFunctionSpace<
+                                        c32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    DType::C64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_into_dense_internal_complex::<
+                                    c64,
+                                    SingleLayerBoundaryIntegrand<c64>,
+                                    KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                                    SerialFunctionSpace<
+                                        c64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::C32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_into_dense_internal_complex::<
+                                    c32,
+                                    DoubleLayerBoundaryIntegrand<c32>,
+                                    KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                                    SerialFunctionSpace<
+                                        c32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    DType::C64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_into_dense_internal_complex::<
+                                    c64,
+                                    DoubleLayerBoundaryIntegrand<c64>,
+                                    KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                                    SerialFunctionSpace<
+                                        c64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::AdjointDoubleLayer => match (*assembler).dtype {
+                    DType::C32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_into_dense_internal_complex::<
+                                    c32,
+                                    AdjointDoubleLayerBoundaryIntegrand<c32>,
+                                    KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                                    SerialFunctionSpace<
+                                        c32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    DType::C64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_into_dense_internal_complex::<
+                                    c64,
+                                    AdjointDoubleLayerBoundaryIntegrand<c64>,
+                                    KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                                    SerialFunctionSpace<
+                                        c64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::Hypersingular => match (*assembler).dtype {
+                    DType::C32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_into_dense_internal_complex::<
+                                    c32,
+                                    HelmholtzHypersingularBoundaryIntegrand<c32>,
+                                    KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                                    SerialFunctionSpace<
+                                        c32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    DType::C64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_into_dense_internal_complex::<
+                                    c64,
+                                    HelmholtzHypersingularBoundaryIntegrand<c64>,
+                                    KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                                    SerialFunctionSpace<
+                                        c64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+        }
+    }
+
+    unsafe fn boundary_assembler_assemble_nonsingular_into_dense_internal_real<
+        T: RlstScalar<Real = T> + MatrixInverse,
+        Integrand: BoundaryIntegrand<T = T>,
+        Kernel: KernelEvaluatorTrait<T = T>,
+        Space: FunctionSpace<T = T> + Sync,
+    >(
+        assembler: *mut BoundaryAssemblerWrapper,
+        output: *mut c_void,
+        trial_space: *const FunctionSpaceWrapper,
+        test_space: *const FunctionSpaceWrapper,
+    ) {
+        let trial_dim = (*extract_space::<Space>(trial_space)).global_size();
+        let test_dim = (*extract_space::<Space>(test_space)).global_size();
+        let mut output = rlst_array_from_slice_mut2!(
+            from_raw_parts_mut(output as *mut T, trial_dim * test_dim),
+            [test_dim, trial_dim]
+        );
+        let trial_colouring = (*extract_space::<Space>(trial_space)).cell_colouring();
+        let test_colouring = (*extract_space::<Space>(test_space)).cell_colouring();
+
+        (*extract_boundary_assembler::<T, Integrand, Kernel>(assembler))
+            .assemble_nonsingular_into_dense(
+                &mut output,
+                &*extract_space::<Space>(trial_space),
+                &*extract_space::<Space>(test_space),
+                &trial_colouring,
+                &test_colouring,
+            )
+    }
+    unsafe fn boundary_assembler_assemble_nonsingular_into_dense_internal_complex<
+        T: RlstScalar<Complex = T> + MatrixInverse,
+        Integrand: BoundaryIntegrand<T = T>,
+        Kernel: KernelEvaluatorTrait<T = T>,
+        Space: FunctionSpace<T = T> + Sync,
+    >(
+        assembler: *mut BoundaryAssemblerWrapper,
+        output: *mut c_void,
+        trial_space: *const FunctionSpaceWrapper,
+        test_space: *const FunctionSpaceWrapper,
+    ) {
+        let trial_dim = (*extract_space::<Space>(trial_space)).global_size();
+        let test_dim = (*extract_space::<Space>(test_space)).global_size();
+        let mut output = rlst_array_from_slice_mut2!(
+            from_raw_parts_mut(output as *mut T, trial_dim * test_dim),
+            [test_dim, trial_dim]
+        );
+        let trial_colouring = (*extract_space::<Space>(trial_space)).cell_colouring();
+        let test_colouring = (*extract_space::<Space>(test_space)).cell_colouring();
+
+        (*extract_boundary_assembler::<T, Integrand, Kernel>(assembler))
+            .assemble_nonsingular_into_dense(
+                &mut output,
+                &*extract_space::<Space>(trial_space),
+                &*extract_space::<Space>(test_space),
+                &trial_colouring,
+                &test_colouring,
+            )
+    }
+    #[no_mangle]
+    pub unsafe extern "C" fn boundary_assembler_assemble_nonsingular_into_dense(
+        assembler: *mut BoundaryAssemblerWrapper,
+        output: *mut c_void,
+        trial_space: *const FunctionSpaceWrapper,
+        test_space: *const FunctionSpaceWrapper,
+    ) {
+        match (*assembler).dtype {
+            DType::F32 => {
+                assert_eq!((*trial_space).dtype, DType::F32);
+                assert_eq!((*test_space).dtype, DType::F32);
+            }
+            DType::F64 => {
+                assert_eq!((*trial_space).dtype, DType::F64);
+                assert_eq!((*test_space).dtype, DType::F64);
+            }
+            DType::C32 => {
+                assert_eq!((*trial_space).dtype, DType::F32);
+                assert_eq!((*test_space).dtype, DType::F32);
+            }
+            DType::C64 => {
+                assert_eq!((*trial_space).dtype, DType::F64);
+                assert_eq!((*test_space).dtype, DType::F64);
+            }
+        }
+        assert_eq!((*test_space).stype, (*trial_space).stype);
+        assert_eq!((*test_space).gtype, (*trial_space).gtype);
+        match (*assembler).ktype {
+            KernelType::Laplace => match (*assembler).itype {
+                BoundaryOperator::SingleLayer => match (*assembler).dtype {
+                    DType::F32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_nonsingular_into_dense_internal_real::<
+                                    f32,
+                                    SingleLayerBoundaryIntegrand<f32>,
+                                    KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                                    SerialFunctionSpace<
+                                        f32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    DType::F64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_nonsingular_into_dense_internal_real::<
+                                    f64,
+                                    SingleLayerBoundaryIntegrand<f64>,
+                                    KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                                    SerialFunctionSpace<
+                                        f64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::F32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_nonsingular_into_dense_internal_real::<
+                                    f32,
+                                    DoubleLayerBoundaryIntegrand<f32>,
+                                    KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                                    SerialFunctionSpace<
+                                        f32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    DType::F64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_nonsingular_into_dense_internal_real::<
+                                    f64,
+                                    DoubleLayerBoundaryIntegrand<f64>,
+                                    KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                                    SerialFunctionSpace<
+                                        f64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::AdjointDoubleLayer => match (*assembler).dtype {
+                    DType::F32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_nonsingular_into_dense_internal_real::<
+                                    f32,
+                                    AdjointDoubleLayerBoundaryIntegrand<f32>,
+                                    KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                                    SerialFunctionSpace<
+                                        f32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    DType::F64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_nonsingular_into_dense_internal_real::<
+                                    f64,
+                                    AdjointDoubleLayerBoundaryIntegrand<f64>,
+                                    KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                                    SerialFunctionSpace<
+                                        f64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::Hypersingular => match (*assembler).dtype {
+                    DType::F32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_nonsingular_into_dense_internal_real::<
+                                    f32,
+                                    LaplaceHypersingularBoundaryIntegrand<f32>,
+                                    KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                                    SerialFunctionSpace<
+                                        f32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    DType::F64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_nonsingular_into_dense_internal_real::<
+                                    f64,
+                                    LaplaceHypersingularBoundaryIntegrand<f64>,
+                                    KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                                    SerialFunctionSpace<
+                                        f64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+            KernelType::Helmholtz => match (*assembler).itype {
+                BoundaryOperator::SingleLayer => match (*assembler).dtype {
+                    DType::C32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_nonsingular_into_dense_internal_complex::<
+                                    c32,
+                                    SingleLayerBoundaryIntegrand<c32>,
+                                    KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                                    SerialFunctionSpace<
+                                        c32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    DType::C64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_nonsingular_into_dense_internal_complex::<
+                                    c64,
+                                    SingleLayerBoundaryIntegrand<c64>,
+                                    KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                                    SerialFunctionSpace<
+                                        c64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::C32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_nonsingular_into_dense_internal_complex::<
+                                    c32,
+                                    DoubleLayerBoundaryIntegrand<c32>,
+                                    KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                                    SerialFunctionSpace<
+                                        c32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    DType::C64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_nonsingular_into_dense_internal_complex::<
+                                    c64,
+                                    DoubleLayerBoundaryIntegrand<c64>,
+                                    KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                                    SerialFunctionSpace<
+                                        c64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::AdjointDoubleLayer => match (*assembler).dtype {
+                    DType::C32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_nonsingular_into_dense_internal_complex::<
+                                    c32,
+                                    AdjointDoubleLayerBoundaryIntegrand<c32>,
+                                    KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                                    SerialFunctionSpace<
+                                        c32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    DType::C64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_nonsingular_into_dense_internal_complex::<
+                                    c64,
+                                    AdjointDoubleLayerBoundaryIntegrand<c64>,
+                                    KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                                    SerialFunctionSpace<
+                                        c64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::Hypersingular => match (*assembler).dtype {
+                    DType::C32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_nonsingular_into_dense_internal_complex::<
+                                    c32,
+                                    HelmholtzHypersingularBoundaryIntegrand<c32>,
+                                    KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                                    SerialFunctionSpace<
+                                        c32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    DType::C64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_nonsingular_into_dense_internal_complex::<
+                                    c64,
+                                    HelmholtzHypersingularBoundaryIntegrand<c64>,
+                                    KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                                    SerialFunctionSpace<
+                                        c64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(
+                                    assembler, output, trial_space, test_space
+                                );
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+        }
+    }
+
+    unsafe fn boundary_assembler_assemble_singular_internal_real<
+        T: RlstScalar<Real = T> + MatrixInverse,
+        Integrand: BoundaryIntegrand<T = T>,
+        Kernel: KernelEvaluatorTrait<T = T>,
+        Space: FunctionSpace<T = T> + Sync,
+    >(
+        assembler: *mut BoundaryAssemblerWrapper,
+        trial_space: *const FunctionSpaceWrapper,
+        test_space: *const FunctionSpaceWrapper,
+    ) -> *const SparseMatrixWrapper {
+        let trial_dim = (*extract_space::<Space>(trial_space)).global_size();
+        let test_dim = (*extract_space::<Space>(test_space)).global_size();
+        Box::into_raw(Box::new(SparseMatrixWrapper {
+            matrix: Box::into_raw(Box::new(
+                (*extract_boundary_assembler::<T, Integrand, Kernel>(assembler)).assemble_singular(
+                    [test_dim, trial_dim],
+                    &*extract_space::<Space>(trial_space),
+                    &*extract_space::<Space>(test_space),
+                ),
+            )) as *const c_void,
+            dtype: (*assembler).dtype,
+        }))
+    }
+    unsafe fn boundary_assembler_assemble_singular_internal_complex<
+        T: RlstScalar<Complex = T> + MatrixInverse,
+        Integrand: BoundaryIntegrand<T = T>,
+        Kernel: KernelEvaluatorTrait<T = T>,
+        Space: FunctionSpace<T = T> + Sync,
+    >(
+        assembler: *mut BoundaryAssemblerWrapper,
+        trial_space: *const FunctionSpaceWrapper,
+        test_space: *const FunctionSpaceWrapper,
+    ) -> *const SparseMatrixWrapper {
+        let trial_dim = (*extract_space::<Space>(trial_space)).global_size();
+        let test_dim = (*extract_space::<Space>(test_space)).global_size();
+        Box::into_raw(Box::new(SparseMatrixWrapper {
+            matrix: Box::into_raw(Box::new(
+                (*extract_boundary_assembler::<T, Integrand, Kernel>(assembler)).assemble_singular(
+                    [test_dim, trial_dim],
+                    &*extract_space::<Space>(trial_space),
+                    &*extract_space::<Space>(test_space),
+                ),
+            )) as *const c_void,
+            dtype: (*assembler).dtype,
+        }))
+    }
+    #[no_mangle]
+    pub unsafe extern "C" fn boundary_assembler_assemble_singular(
+        assembler: *mut BoundaryAssemblerWrapper,
+        trial_space: *const FunctionSpaceWrapper,
+        test_space: *const FunctionSpaceWrapper,
+    ) -> *const SparseMatrixWrapper {
+        match (*assembler).dtype {
+            DType::F32 => {
+                assert_eq!((*trial_space).dtype, DType::F32);
+                assert_eq!((*test_space).dtype, DType::F32);
+            }
+            DType::F64 => {
+                assert_eq!((*trial_space).dtype, DType::F64);
+                assert_eq!((*test_space).dtype, DType::F64);
+            }
+            DType::C32 => {
+                assert_eq!((*trial_space).dtype, DType::F32);
+                assert_eq!((*test_space).dtype, DType::F32);
+            }
+            DType::C64 => {
+                assert_eq!((*trial_space).dtype, DType::F64);
+                assert_eq!((*test_space).dtype, DType::F64);
+            }
+        }
+        assert_eq!((*test_space).stype, (*trial_space).stype);
+        assert_eq!((*test_space).gtype, (*trial_space).gtype);
+        match (*assembler).ktype {
+            KernelType::Laplace => match (*assembler).itype {
+                BoundaryOperator::SingleLayer => match (*assembler).dtype {
+                    DType::F32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_internal_real::<
+                                    f32,
+                                    SingleLayerBoundaryIntegrand<f32>,
+                                    KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                                    SerialFunctionSpace<
+                                        f32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(assembler, trial_space, test_space)
+                            }
+                        },
+                    },
+                    DType::F64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_internal_real::<
+                                    f64,
+                                    SingleLayerBoundaryIntegrand<f64>,
+                                    KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                                    SerialFunctionSpace<
+                                        f64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(assembler, trial_space, test_space)
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::F32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_internal_real::<
+                                    f32,
+                                    DoubleLayerBoundaryIntegrand<f32>,
+                                    KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                                    SerialFunctionSpace<
+                                        f32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(assembler, trial_space, test_space)
+                            }
+                        },
+                    },
+                    DType::F64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_internal_real::<
+                                    f64,
+                                    DoubleLayerBoundaryIntegrand<f64>,
+                                    KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                                    SerialFunctionSpace<
+                                        f64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(assembler, trial_space, test_space)
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::AdjointDoubleLayer => match (*assembler).dtype {
+                    DType::F32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_internal_real::<
+                                    f32,
+                                    AdjointDoubleLayerBoundaryIntegrand<f32>,
+                                    KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                                    SerialFunctionSpace<
+                                        f32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(assembler, trial_space, test_space)
+                            }
+                        },
+                    },
+                    DType::F64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_internal_real::<
+                                    f64,
+                                    AdjointDoubleLayerBoundaryIntegrand<f64>,
+                                    KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                                    SerialFunctionSpace<
+                                        f64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(assembler, trial_space, test_space)
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::Hypersingular => match (*assembler).dtype {
+                    DType::F32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_internal_real::<
+                                    f32,
+                                    LaplaceHypersingularBoundaryIntegrand<f32>,
+                                    KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                                    SerialFunctionSpace<
+                                        f32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(assembler, trial_space, test_space)
+                            }
+                        },
+                    },
+                    DType::F64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_internal_real::<
+                                    f64,
+                                    LaplaceHypersingularBoundaryIntegrand<f64>,
+                                    KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                                    SerialFunctionSpace<
+                                        f64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(assembler, trial_space, test_space)
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+            KernelType::Helmholtz => match (*assembler).itype {
+                BoundaryOperator::SingleLayer => match (*assembler).dtype {
+                    DType::C32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_internal_complex::<
+                                    c32,
+                                    SingleLayerBoundaryIntegrand<c32>,
+                                    KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                                    SerialFunctionSpace<
+                                        c32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(assembler, trial_space, test_space)
+                            }
+                        },
+                    },
+                    DType::C64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_internal_complex::<
+                                    c64,
+                                    SingleLayerBoundaryIntegrand<c64>,
+                                    KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                                    SerialFunctionSpace<
+                                        c64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(assembler, trial_space, test_space)
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::C32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_internal_complex::<
+                                    c32,
+                                    DoubleLayerBoundaryIntegrand<c32>,
+                                    KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                                    SerialFunctionSpace<
+                                        c32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(assembler, trial_space, test_space)
+                            }
+                        },
+                    },
+                    DType::C64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_internal_complex::<
+                                    c64,
+                                    DoubleLayerBoundaryIntegrand<c64>,
+                                    KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                                    SerialFunctionSpace<
+                                        c64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(assembler, trial_space, test_space)
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::AdjointDoubleLayer => match (*assembler).dtype {
+                    DType::C32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_internal_complex::<
+                                    c32,
+                                    AdjointDoubleLayerBoundaryIntegrand<c32>,
+                                    KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                                    SerialFunctionSpace<
+                                        c32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(assembler, trial_space, test_space)
+                            }
+                        },
+                    },
+                    DType::C64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_internal_complex::<
+                                    c64,
+                                    AdjointDoubleLayerBoundaryIntegrand<c64>,
+                                    KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                                    SerialFunctionSpace<
+                                        c64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(assembler, trial_space, test_space)
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::Hypersingular => match (*assembler).dtype {
+                    DType::C32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_internal_complex::<
+                                    c32,
+                                    HelmholtzHypersingularBoundaryIntegrand<c32>,
+                                    KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                                    SerialFunctionSpace<
+                                        c32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(assembler, trial_space, test_space)
+                            }
+                        },
+                    },
+                    DType::C64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_internal_complex::<
+                                    c64,
+                                    HelmholtzHypersingularBoundaryIntegrand<c64>,
+                                    KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                                    SerialFunctionSpace<
+                                        c64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(assembler, trial_space, test_space)
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+        }
+    }
+
+    unsafe fn boundary_assembler_assemble_singular_correction_internal_real<
+        T: RlstScalar<Real = T> + MatrixInverse,
+        Integrand: BoundaryIntegrand<T = T>,
+        Kernel: KernelEvaluatorTrait<T = T>,
+        Space: FunctionSpace<T = T> + Sync,
+    >(
+        assembler: *mut BoundaryAssemblerWrapper,
+        trial_space: *const FunctionSpaceWrapper,
+        test_space: *const FunctionSpaceWrapper,
+    ) -> *const SparseMatrixWrapper {
+        let trial_dim = (*extract_space::<Space>(trial_space)).global_size();
+        let test_dim = (*extract_space::<Space>(test_space)).global_size();
+        Box::into_raw(Box::new(SparseMatrixWrapper {
+            matrix: Box::into_raw(Box::new(
+                (*extract_boundary_assembler::<T, Integrand, Kernel>(assembler))
+                    .assemble_singular_correction(
+                        [test_dim, trial_dim],
+                        &*extract_space::<Space>(trial_space),
+                        &*extract_space::<Space>(test_space),
+                    ),
+            )) as *const c_void,
+            dtype: (*assembler).dtype,
+        }))
+    }
+    unsafe fn boundary_assembler_assemble_singular_correction_internal_complex<
+        T: RlstScalar<Complex = T> + MatrixInverse,
+        Integrand: BoundaryIntegrand<T = T>,
+        Kernel: KernelEvaluatorTrait<T = T>,
+        Space: FunctionSpace<T = T> + Sync,
+    >(
+        assembler: *mut BoundaryAssemblerWrapper,
+        trial_space: *const FunctionSpaceWrapper,
+        test_space: *const FunctionSpaceWrapper,
+    ) -> *const SparseMatrixWrapper {
+        let trial_dim = (*extract_space::<Space>(trial_space)).global_size();
+        let test_dim = (*extract_space::<Space>(test_space)).global_size();
+        Box::into_raw(Box::new(SparseMatrixWrapper {
+            matrix: Box::into_raw(Box::new(
+                (*extract_boundary_assembler::<T, Integrand, Kernel>(assembler))
+                    .assemble_singular_correction(
+                        [test_dim, trial_dim],
+                        &*extract_space::<Space>(trial_space),
+                        &*extract_space::<Space>(test_space),
+                    ),
+            )) as *const c_void,
+            dtype: (*assembler).dtype,
+        }))
+    }
+    #[no_mangle]
+    pub unsafe extern "C" fn boundary_assembler_assemble_singular_correction(
+        assembler: *mut BoundaryAssemblerWrapper,
+        trial_space: *const FunctionSpaceWrapper,
+        test_space: *const FunctionSpaceWrapper,
+    ) -> *const SparseMatrixWrapper {
+        match (*assembler).dtype {
+            DType::F32 => {
+                assert_eq!((*trial_space).dtype, DType::F32);
+                assert_eq!((*test_space).dtype, DType::F32);
+            }
+            DType::F64 => {
+                assert_eq!((*trial_space).dtype, DType::F64);
+                assert_eq!((*test_space).dtype, DType::F64);
+            }
+            DType::C32 => {
+                assert_eq!((*trial_space).dtype, DType::F32);
+                assert_eq!((*test_space).dtype, DType::F32);
+            }
+            DType::C64 => {
+                assert_eq!((*trial_space).dtype, DType::F64);
+                assert_eq!((*test_space).dtype, DType::F64);
+            }
+        }
+        assert_eq!((*test_space).stype, (*trial_space).stype);
+        assert_eq!((*test_space).gtype, (*trial_space).gtype);
+        match (*assembler).ktype {
+            KernelType::Laplace => match (*assembler).itype {
+                BoundaryOperator::SingleLayer => match (*assembler).dtype {
+                    DType::F32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_correction_internal_real::<
+                                    f32,
+                                    SingleLayerBoundaryIntegrand<f32>,
+                                    KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                                    SerialFunctionSpace<
+                                        f32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(assembler, trial_space, test_space)
+                            }
+                        },
+                    },
+                    DType::F64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_correction_internal_real::<
+                                    f64,
+                                    SingleLayerBoundaryIntegrand<f64>,
+                                    KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                                    SerialFunctionSpace<
+                                        f64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(assembler, trial_space, test_space)
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::F32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_correction_internal_real::<
+                                    f32,
+                                    DoubleLayerBoundaryIntegrand<f32>,
+                                    KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                                    SerialFunctionSpace<
+                                        f32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(assembler, trial_space, test_space)
+                            }
+                        },
+                    },
+                    DType::F64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_correction_internal_real::<
+                                    f64,
+                                    DoubleLayerBoundaryIntegrand<f64>,
+                                    KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                                    SerialFunctionSpace<
+                                        f64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(assembler, trial_space, test_space)
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::AdjointDoubleLayer => match (*assembler).dtype {
+                    DType::F32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_correction_internal_real::<
+                                    f32,
+                                    AdjointDoubleLayerBoundaryIntegrand<f32>,
+                                    KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                                    SerialFunctionSpace<
+                                        f32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(assembler, trial_space, test_space)
+                            }
+                        },
+                    },
+                    DType::F64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_correction_internal_real::<
+                                    f64,
+                                    AdjointDoubleLayerBoundaryIntegrand<f64>,
+                                    KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                                    SerialFunctionSpace<
+                                        f64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(assembler, trial_space, test_space)
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::Hypersingular => match (*assembler).dtype {
+                    DType::F32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_correction_internal_real::<
+                                    f32,
+                                    LaplaceHypersingularBoundaryIntegrand<f32>,
+                                    KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                                    SerialFunctionSpace<
+                                        f32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(assembler, trial_space, test_space)
+                            }
+                        },
+                    },
+                    DType::F64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_correction_internal_real::<
+                                    f64,
+                                    LaplaceHypersingularBoundaryIntegrand<f64>,
+                                    KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                                    SerialFunctionSpace<
+                                        f64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(assembler, trial_space, test_space)
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+            KernelType::Helmholtz => match (*assembler).itype {
+                BoundaryOperator::SingleLayer => match (*assembler).dtype {
+                    DType::C32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_correction_internal_complex::<
+                                    c32,
+                                    SingleLayerBoundaryIntegrand<c32>,
+                                    KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                                    SerialFunctionSpace<
+                                        c32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(assembler, trial_space, test_space)
+                            }
+                        },
+                    },
+                    DType::C64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_correction_internal_complex::<
+                                    c64,
+                                    SingleLayerBoundaryIntegrand<c64>,
+                                    KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                                    SerialFunctionSpace<
+                                        c64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(assembler, trial_space, test_space)
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::C32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_correction_internal_complex::<
+                                    c32,
+                                    DoubleLayerBoundaryIntegrand<c32>,
+                                    KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                                    SerialFunctionSpace<
+                                        c32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(assembler, trial_space, test_space)
+                            }
+                        },
+                    },
+                    DType::C64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_correction_internal_complex::<
+                                    c64,
+                                    DoubleLayerBoundaryIntegrand<c64>,
+                                    KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                                    SerialFunctionSpace<
+                                        c64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(assembler, trial_space, test_space)
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::AdjointDoubleLayer => match (*assembler).dtype {
+                    DType::C32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_correction_internal_complex::<
+                                    c32,
+                                    AdjointDoubleLayerBoundaryIntegrand<c32>,
+                                    KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                                    SerialFunctionSpace<
+                                        c32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(assembler, trial_space, test_space)
+                            }
+                        },
+                    },
+                    DType::C64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_correction_internal_complex::<
+                                    c64,
+                                    AdjointDoubleLayerBoundaryIntegrand<c64>,
+                                    KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                                    SerialFunctionSpace<
+                                        c64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(assembler, trial_space, test_space)
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                BoundaryOperator::Hypersingular => match (*assembler).dtype {
+                    DType::C32 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_correction_internal_complex::<
+                                    c32,
+                                    HelmholtzHypersingularBoundaryIntegrand<c32>,
+                                    KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                                    SerialFunctionSpace<
+                                        c32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(assembler, trial_space, test_space)
+                            }
+                        },
+                    },
+                    DType::C64 => match (*trial_space).stype {
+                        SpaceType::SerialFunctionSpace => match (*trial_space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                boundary_assembler_assemble_singular_correction_internal_complex::<
+                                    c64,
+                                    HelmholtzHypersingularBoundaryIntegrand<c64>,
+                                    KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                                    SerialFunctionSpace<
+                                        c64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(assembler, trial_space, test_space)
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+        }
+    }
+
+    #[no_mangle]
+    pub unsafe extern "C" fn boundary_assembler_dtype(
+        assembler: *mut BoundaryAssemblerWrapper,
+    ) -> u8 {
+        (*assembler).dtype as u8
+    }
+
+    unsafe fn laplace_boundary_assembler_new_internal<T: RlstScalar + MatrixInverse>(
+        operator: BoundaryOperator,
+        dtype: DType,
+    ) -> *const BoundaryAssemblerWrapper {
+        Box::into_raw(Box::new(BoundaryAssemblerWrapper {
+            assembler: match operator {
+                BoundaryOperator::SingleLayer => Box::into_raw(Box::new(
+                    BoundaryAssembler::<T, _, _>::new_laplace_single_layer(),
+                )) as *const c_void,
+                BoundaryOperator::DoubleLayer => Box::into_raw(Box::new(
+                    BoundaryAssembler::<T, _, _>::new_laplace_double_layer(),
+                )) as *const c_void,
+                BoundaryOperator::AdjointDoubleLayer => Box::into_raw(Box::new(
+                    BoundaryAssembler::<T, _, _>::new_laplace_adjoint_double_layer(),
+                )) as *const c_void,
+                BoundaryOperator::Hypersingular => Box::into_raw(Box::new(
+                    BoundaryAssembler::<T, _, _>::new_laplace_hypersingular(),
+                )) as *const c_void,
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+            itype: operator,
+            ktype: KernelType::Laplace,
+            dtype,
+        }))
+    }
+
+    #[no_mangle]
+    pub unsafe extern "C" fn laplace_boundary_assembler_new(
+        operator: u8,
+        dtype: u8,
+    ) -> *const BoundaryAssemblerWrapper {
+        let operator = BoundaryOperator::from(operator).unwrap();
+        let dtype = DType::from(dtype).unwrap();
+        match dtype {
+            DType::F32 => laplace_boundary_assembler_new_internal::<f32>(operator, dtype),
+            DType::F64 => laplace_boundary_assembler_new_internal::<f64>(operator, dtype),
+            _ => {
+                panic!("Invalid data type");
+            }
+        }
+    }
+
+    unsafe fn helmholtz_boundary_assembler_new_internal<
+        T: RlstScalar<Complex = T> + MatrixInverse,
+    >(
+        wavenumber: T::Real,
+        operator: BoundaryOperator,
+        dtype: DType,
+    ) -> *const BoundaryAssemblerWrapper {
+        Box::into_raw(Box::new(BoundaryAssemblerWrapper {
+            assembler: match operator {
+                BoundaryOperator::SingleLayer => Box::into_raw(Box::new(
+                    BoundaryAssembler::<T, _, _>::new_helmholtz_single_layer(wavenumber),
+                )) as *const c_void,
+                BoundaryOperator::DoubleLayer => Box::into_raw(Box::new(
+                    BoundaryAssembler::<T, _, _>::new_helmholtz_double_layer(wavenumber),
+                )) as *const c_void,
+                BoundaryOperator::AdjointDoubleLayer => Box::into_raw(Box::new(
+                    BoundaryAssembler::<T, _, _>::new_helmholtz_adjoint_double_layer(wavenumber),
+                )) as *const c_void,
+                BoundaryOperator::Hypersingular => Box::into_raw(Box::new(
+                    BoundaryAssembler::<T, _, _>::new_helmholtz_hypersingular(wavenumber),
+                )) as *const c_void,
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+            itype: operator,
+            ktype: KernelType::Helmholtz,
+            dtype,
+        }))
+    }
+
+    #[no_mangle]
+    pub unsafe extern "C" fn helmholtz_boundary_assembler_new(
+        wavenumber: *const c_void,
+        operator: u8,
+        dtype: u8,
+    ) -> *const BoundaryAssemblerWrapper {
+        let operator = BoundaryOperator::from(operator).unwrap();
+        let dtype = DType::from(dtype).unwrap();
+        match dtype {
+            DType::C32 => helmholtz_boundary_assembler_new_internal::<c32>(
+                *(wavenumber as *const f32),
+                operator,
+                dtype,
+            ),
+            DType::C64 => helmholtz_boundary_assembler_new_internal::<c64>(
+                *(wavenumber as *const f64),
+                operator,
+                dtype,
+            ),
+            _ => {
+                panic!("Invalid data type");
+            }
+        }
+    }
+}
diff --git a/src/traits/assembly/boundary.rs b/src/traits/assembly/boundary.rs
index 9d235f37..35eb02e5 100644
--- a/src/traits/assembly/boundary.rs
+++ b/src/traits/assembly/boundary.rs
@@ -10,7 +10,7 @@ use crate::traits::ParallelFunctionSpace;
 #[cfg(feature = "mpi")]
 use mpi::traits::Communicator;
 use ndelement::types::ReferenceCellType;
-use rlst::{CsrMatrix, RlstScalar};
+use rlst::{CsrMatrix, RandomAccessMut, RawAccessMut, RlstScalar, Shape};
 use std::collections::HashMap;
 
 pub trait CellPairAssembler {
@@ -32,9 +32,12 @@ pub trait BoundaryAssembly {
     type T: RlstScalar;
 
     /// Assemble the singular contributions into a dense matrix
-    fn assemble_singular_into_dense<Space: FunctionSpace<T = Self::T> + Sync>(
+    fn assemble_singular_into_dense<
+        Space: FunctionSpace<T = Self::T> + Sync,
+        Array2: RandomAccessMut<2, Item = Self::T> + Shape<2> + RawAccessMut<Item = Self::T>,
+    >(
         &self,
-        output: &mut RlstArray<Self::T, 2>,
+        output: &mut Array2,
         trial_space: &Space,
         test_space: &Space,
     );
@@ -49,9 +52,12 @@ pub trait BoundaryAssembly {
     /// Assemble the singular correction into a dense matrix
     ///
     /// 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_into_dense<Space: FunctionSpace<T = Self::T> + Sync>(
+    fn assemble_singular_correction_into_dense<
+        Space: FunctionSpace<T = Self::T> + Sync,
+        Array2: RandomAccessMut<2, Item = Self::T> + Shape<2> + RawAccessMut<Item = Self::T>,
+    >(
         &self,
-        output: &mut RlstArray<Self::T, 2>,
+        output: &mut Array2,
         trial_space: &Space,
         test_space: &Space,
     );
@@ -66,17 +72,23 @@ pub trait BoundaryAssembly {
     ) -> CsrMatrix<Self::T>;
 
     /// Assemble into a dense matrix
-    fn assemble_into_dense<Space: FunctionSpace<T = Self::T> + Sync>(
+    fn assemble_into_dense<
+        Space: FunctionSpace<T = Self::T> + Sync,
+        Array2: RandomAccessMut<2, Item = Self::T> + Shape<2> + RawAccessMut<Item = Self::T>,
+    >(
         &self,
-        output: &mut RlstArray<Self::T, 2>,
+        output: &mut Array2,
         trial_space: &Space,
         test_space: &Space,
     );
 
     /// Assemble the non-singular contributions into a dense matrix
-    fn assemble_nonsingular_into_dense<Space: FunctionSpace<T = Self::T> + Sync>(
+    fn assemble_nonsingular_into_dense<
+        Space: FunctionSpace<T = Self::T> + Sync,
+        Array2: RandomAccessMut<2, Item = Self::T> + Shape<2> + RawAccessMut<Item = Self::T>,
+    >(
         &self,
-        output: &mut RlstArray<Self::T, 2>,
+        output: &mut Array2,
         trial_space: &Space,
         test_space: &Space,
         trial_colouring: &HashMap<ReferenceCellType, Vec<Vec<usize>>>,


Test patch:
diff --git a/.github/workflows/run-tests.yml b/.github/workflows/run-tests.yml
index 13769c6e..35f87e8b 100644
--- a/.github/workflows/run-tests.yml
+++ b/.github/workflows/run-tests.yml
@@ -54,7 +54,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: ["3.10", "3.11", "3.12"]
+        python-version: ["3.8", "3.9", "3.10", "3.11", "3.12"]
     steps:
       - name: Set up Python
         uses: actions/setup-python@v4
@@ -70,7 +70,7 @@ jobs:
       - name: Install python package
         run: |
           source .venv/bin/activate
-          maturin develop
+          maturin develop --release
       - name: Run Python tests
         run: |
           source .venv/bin/activate
diff --git a/.github/workflows/run-weekly-tests.yml b/.github/workflows/run-weekly-tests.yml
index f9a26445..0d910fd8 100644
--- a/.github/workflows/run-weekly-tests.yml
+++ b/.github/workflows/run-weekly-tests.yml
@@ -57,3 +57,31 @@ jobs:
 
       - name: Build docs
         run: cargo doc --features "mpi,strict" --no-deps
+
+  run-tests-python:
+    name: Run Python tests
+    runs-on: ubuntu-latest
+    strategy:
+      matrix:
+        python-version: ["3.8", "3.9", "3.10", "3.11", "3.12"]
+    steps:
+      - name: Set up Python
+        uses: actions/setup-python@v4
+        with:
+          python-version: ${{ matrix.python-version }}
+      - uses: actions/checkout@v3
+      - name: Install uv
+        run: pip install uv "maturin>=1.7"
+      - name: Make virtual environment
+        run: |
+          uv venv .venv
+          uv pip install pip pytest
+      - name: Install python package
+        run: |
+          source .venv/bin/activate
+          maturin develop --release
+      - name: Run Python tests
+        run: |
+          source .venv/bin/activate
+          python -m pytest python/test
+
diff --git a/python/test/test_boundary_assembler.py b/python/test/test_boundary_assembler.py
new file mode 100644
index 00000000..5e75bede
--- /dev/null
+++ b/python/test/test_boundary_assembler.py
@@ -0,0 +1,272 @@
+import pytest
+import numpy as np
+from ndelement.reference_cell import ReferenceCellType
+from bempp.assembly.boundary import OperatorType, create_laplace_assembler
+from bempp.function_space import function_space
+from ndgrid.shapes import regular_sphere
+from ndelement.ciarlet import create_family, Family, Continuity
+
+
+@pytest.mark.parametrize(
+    "otype",
+    [
+        OperatorType.SingleLayer,
+        OperatorType.DoubleLayer,
+        OperatorType.AdjointDoubleLayer,
+        OperatorType.Hypersingular,
+    ],
+)
+def test_create_assembler(otype):
+    a = create_laplace_assembler(otype)
+
+    assert a.quadrature_degree(ReferenceCellType.Triangle) != 3
+    a.set_quadrature_degree(ReferenceCellType.Triangle, 3)
+    assert a.quadrature_degree(ReferenceCellType.Triangle) == 3
+    with pytest.raises(ValueError):
+        a.set_quadrature_degree(ReferenceCellType.Interval, 3)
+    with pytest.raises(ValueError):
+        a.quadrature_degree(ReferenceCellType.Interval)
+
+    assert a.singular_quadrature_degree(ReferenceCellType.Triangle, ReferenceCellType.Triangle) != 3
+    a.set_singular_quadrature_degree(ReferenceCellType.Triangle, ReferenceCellType.Triangle, 3)
+    assert a.singular_quadrature_degree(ReferenceCellType.Triangle, ReferenceCellType.Triangle) == 3
+    with pytest.raises(ValueError):
+        a.set_singular_quadrature_degree(ReferenceCellType.Interval, ReferenceCellType.Interval, 3)
+
+    assert a.batch_size != 4
+    a.set_batch_size(4)
+    assert a.batch_size == 4
+
+    assert a.dtype == np.float64
+
+
+def test_single_layer_sphere0_dp0():
+    grid = regular_sphere(0)
+    element = create_family(Family.Lagrange, 0, Continuity.Discontinuous)
+    space = function_space(grid, element)
+
+    a = create_laplace_assembler(OperatorType.SingleLayer)
+
+    mat = a.assemble_into_dense(space, space)
+
+    from_cl = np.array(
+        [
+            [
+                0.1854538822982487,
+                0.08755414595678074,
+                0.05963897421514472,
+                0.08755414595678074,
+                0.08755414595678074,
+                0.05963897421514473,
+                0.04670742127454548,
+                0.05963897421514472,
+            ],
+            [
+                0.08755414595678074,
+                0.1854538822982487,
+                0.08755414595678074,
+                0.05963897421514472,
+                0.05963897421514472,
+                0.08755414595678074,
+                0.05963897421514473,
+                0.04670742127454548,
+            ],
+            [
+                0.05963897421514472,
+                0.08755414595678074,
+                0.1854538822982487,
+                0.08755414595678074,
+                0.04670742127454548,
+                0.05963897421514472,
+                0.08755414595678074,
+                0.05963897421514473,
+            ],
+            [
+                0.08755414595678074,
+                0.05963897421514472,
+                0.08755414595678074,
+                0.1854538822982487,
+                0.05963897421514473,
+                0.04670742127454548,
+                0.05963897421514472,
+                0.08755414595678074,
+            ],
+            [
+                0.08755414595678074,
+                0.05963897421514472,
+                0.046707421274545476,
+                0.05963897421514473,
+                0.1854538822982487,
+                0.08755414595678074,
+                0.05963897421514472,
+                0.08755414595678074,
+            ],
+            [
+                0.05963897421514473,
+                0.08755414595678074,
+                0.05963897421514472,
+                0.046707421274545476,
+                0.08755414595678074,
+                0.1854538822982487,
+                0.08755414595678074,
+                0.05963897421514472,
+            ],
+            [
+                0.046707421274545476,
+                0.05963897421514473,
+                0.08755414595678074,
+                0.05963897421514472,
+                0.05963897421514472,
+                0.08755414595678074,
+                0.1854538822982487,
+                0.08755414595678074,
+            ],
+            [
+                0.05963897421514472,
+                0.046707421274545476,
+                0.05963897421514473,
+                0.08755414595678074,
+                0.08755414595678074,
+                0.05963897421514472,
+                0.08755414595678074,
+                0.1854538822982487,
+            ],
+        ]
+    )
+
+    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)