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/README.md b/README.md
index a436152f..a63e36fd 100644
--- a/README.md
+++ b/README.md
@@ -3,6 +3,22 @@
 Bempp-rs is an open-source boundary element library written in Rust that can be used to assemble all the standard integral
 kernels for Laplace and Helmholtz problems.
 
+## Using bempp-rs
+### Rust
+You can use the latest release of bempp by adding the following to `[dependencies]` section of your Cargo.toml file:
+
+```toml
+bempp = "0.1.0"
+```
+
+### Python
+You can install the latest release of bempp by running:
+
+```bash
+pip3 install bempp-rs
+```
+
+
 ## Documentation
 The latest documentation of the main branch of this repo is available at [bempp.github.io/bempp-rs/](https://bempp.github.io/bempp-rs/).
 
diff --git a/RELEASE.md b/RELEASE.md
index 943ffacb..4fa504f1 100644
--- a/RELEASE.md
+++ b/RELEASE.md
@@ -5,33 +5,37 @@ To make a new release of bempp, follow the following steps:
 0) If you are yet to make a release on your current computer, run `cargo login` and copy an API
    key from https://crates.io/me
 
-1) Checkout the branch `release` and merge the branch `main` into the branch `release`:
+1) Checkout the `main` branch and `git pull`, then checkout a new branch called `release-v[x].[y].[z]`
+   (where `[x]`, `[y]`, and `[z]` are defined in the next step):
    ```bash
-   git checkout release
-   git merge main
+   git checkout main
+   git pull
+   git checkout -b release-v[x].[y].[z]
    ```
 
-2) Update the version number in `Cargo.toml`.
+2) Update the version numbers in `Cargo.toml` and `pyproject.toml`.
    The version numbers have the format `[x].[y].[z]`. If you are releasing a major
    version, you should increment `[x]` and set `[y]` and `[z]` to 0.
    If you are releasing a minor version, you should increment `[y]` and set `[z]`
    to zero. If you are releasing a bugfix, you should increment `[z]`.
 
-4) In `Cargo.toml`, check that the `rlst`, `green-kernels`, and `kifmm`
-   dependencies are at the latest version.
+3) Update the version number in the "Using bempp" section of README.md.
 
-5) Commit your changes and push to GitHub. Wait to see if the CI tests pass.
+4) Commit your changes and push to GitHub. Wait to see if the CI tests pass.
 
-6) [Create a release on GitHub](https://github.com/bempp/kifmm/releases/new) from the `release` branch.
+5) [Create a release on GitHub](https://github.com/bempp/kifmm/releases/new) from the `main` branch.
    The release tag and title should be `v[x].[y].[z]` (where `[x]`, `[y]` and `[z]` are as in step 2).
    In the "Describe this release" box, you should bullet point the main changes since the last
    release.
 
-7) Run `cargo publish --dry-run`, then run `cargo package --list` and
+6) Run `cargo publish --dry-run`, then run `cargo package --list` and
    check that no unwanted extras have been included in the release.
 
-8) If everything is working as expected, run `cargo publish`. This will push the new version to
+7) If everything is working as expected, run `cargo publish`. This will push the new version to
    crates.io. Note: this cannot be undone, but you can use `cargo yank` to mark a version as
    unsuitable for use.
 
-9) Open a pull request to `main` to update the version number in `Cargo.toml` to `[x].[y].[z]-dev`
+8) Open a pull request to `main` to update the version numbers in `Cargo.toml` and `pyproject.toml`
+   to `[x].[y].[z]-dev`
+
+9) Add the release to the next issue of [Scientific Computing in Rust Monthly](https://github.com/rust-scicomp/scientific-computing-in-rust-monthly)
diff --git a/python/bempp/assembly/potential.py b/python/bempp/assembly/potential.py
new file mode 100644
index 00000000..50a4644e
--- /dev/null
+++ b/python/bempp/assembly/potential.py
@@ -0,0 +1,105 @@
+"""Potential 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
+
+_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")
+
+
+class OperatorType(Enum):
+    """Operator type."""
+
+    SingleLayer = 0
+    DoubleLayer = 1
+
+
+class PotentialAssembler(object):
+    """Potential 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_potential_assembler(self._rs_assembler)
+
+    def assemble_into_dense(
+        self,
+        space: FunctionSpace,
+        points: npt.NDArray[np.floating],
+    ) -> npt.NDArray[np.floating]:
+        """Assemble operator into a dense matrix."""
+        assert space.dtype == points.dtype == self.dtype
+        output = np.zeros((points.shape[0], space.global_size), dtype=self.dtype, order="F")
+        _lib.potential_assembler_assemble_into_dense(
+            self._rs_assembler,
+            _ffi.cast("void*", output.ctypes.data),
+            space._rs_space,
+            _ffi.cast("void*", points.ctypes.data),
+            points.shape[0],
+        )
+        return output
+
+    def set_quadrature_degree(self, cell: ReferenceCellType, degree: int):
+        """Set the (non-singular) quadrature degree for a cell type."""
+        if not _lib.potential_assembler_has_quadrature_degree(self._rs_assembler, cell.value):
+            raise ValueError(f"Invalid cell type: {cell}")
+        _lib.potential_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.potential_assembler_has_quadrature_degree(self._rs_assembler, cell.value):
+            raise ValueError(f"Invalid cell type: {cell}")
+        return _lib.potential_assembler_quadrature_degree(self._rs_assembler, cell.value)
+
+    def set_batch_size(self, batch_size: int):
+        """Set the batch size."""
+        _lib.potential_assembler_set_batch_size(self._rs_assembler, batch_size)
+
+    @property
+    def batch_size(self) -> int:
+        """Get the batch size."""
+        return _lib.potential_assembler_batch_size(self._rs_assembler)
+
+    @property
+    def dtype(self):
+        """Data type."""
+        return _dtypes[_lib.potential_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 PotentialAssembler(
+        _lib.laplace_potential_assembler_new(operator_type.value, _dtype_value(dtype))
+    )
diff --git a/src/assembly/potential.rs b/src/assembly/potential.rs
index ab2ff885..db096252 100644
--- a/src/assembly/potential.rs
+++ b/src/assembly/potential.rs
@@ -1,6 +1,6 @@
 //! Assembly of potential operators
 mod assemblers;
 pub(crate) mod cell_assemblers;
-mod integrands;
+pub mod integrands;
 
 pub use assemblers::PotentialAssembler;
diff --git a/src/assembly/potential/assemblers.rs b/src/assembly/potential/assemblers.rs
index db0f4281..9f7373cf 100644
--- a/src/assembly/potential/assemblers.rs
+++ b/src/assembly/potential/assemblers.rs
@@ -14,8 +14,8 @@ use ndelement::types::ReferenceCellType;
 use ndgrid::traits::Grid;
 use rayon::prelude::*;
 use rlst::{
-    rlst_dynamic_array2, rlst_dynamic_array4, DefaultIterator, MatrixInverse, RandomAccessMut,
-    RawAccess, RawAccessMut, RlstScalar, Shape,
+    rlst_dynamic_array2, rlst_dynamic_array4, DefaultIterator, MatrixInverse, RandomAccessByRef,
+    RandomAccessMut, RawAccess, RawAccessMut, RlstScalar, Shape,
 };
 use std::collections::HashMap;
 
@@ -26,13 +26,14 @@ fn assemble_batch<
     Space: FunctionSpace<T = T> + Sync,
     Integrand: PotentialIntegrand<T = T>,
     Kernel: KernelEvaluator<T = T>,
+    Array2: RandomAccessByRef<2, Item = T::Real> + Shape<2> + RawAccess<Item = T::Real> + Sync,
 >(
     assembler: &PotentialAssembler<T, Integrand, Kernel>,
     deriv_size: usize,
     output: &RawData2D<T>,
     cell_type: ReferenceCellType,
     space: &Space,
-    evaluation_points: &RlstArray<T::Real, 2>,
+    evaluation_points: &Array2,
     cells: &[usize],
     points: &RlstArray<T::Real, 2>,
     weights: &[T::Real],
@@ -135,20 +136,33 @@ 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 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;
     }
 
-    fn assemble<Space: FunctionSpace<T = T> + Sync>(
+    /// 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
+    }
+
+    fn assemble<
+        Space: FunctionSpace<T = T> + Sync,
+        Array2: RandomAccessByRef<2, Item = T::Real> + Shape<2> + RawAccess<Item = T::Real> + Sync,
+    >(
         &self,
         output: &RawData2D<T>,
         space: &Space,
-        points: &RlstArray<T::Real, 2>,
+        points: &Array2,
         colouring: &HashMap<ReferenceCellType, Vec<Vec<usize>>>,
     ) {
         if !space.is_serial() {
@@ -227,11 +241,15 @@ impl<
 {
     type T = T;
 
-    fn assemble_into_dense<Space: FunctionSpace<T = T> + Sync>(
+    fn assemble_into_dense<
+        Space: FunctionSpace<T = T> + Sync,
+        Array2Mut: RandomAccessMut<2, Item = T> + Shape<2> + RawAccessMut<Item = T>,
+        Array2: RandomAccessByRef<2, Item = T::Real> + Shape<2> + RawAccess<Item = T::Real> + Sync,
+    >(
         &self,
-        output: &mut RlstArray<T, 2>,
+        output: &mut Array2Mut,
         space: &Space,
-        points: &RlstArray<T::Real, 2>,
+        points: &Array2,
     ) {
         if !space.is_serial() {
             panic!("Dense assembly can only be used for function spaces stored in serial");
diff --git a/src/assembly/potential/cell_assemblers.rs b/src/assembly/potential/cell_assemblers.rs
index 7c89dc28..7da349ec 100644
--- a/src/assembly/potential/cell_assemblers.rs
+++ b/src/assembly/potential/cell_assemblers.rs
@@ -5,8 +5,8 @@ use crate::traits::{CellAssembler, KernelEvaluator, PotentialIntegrand};
 use ndgrid::traits::GeometryMap;
 use num::Zero;
 use rlst::{
-    rlst_dynamic_array2, rlst_dynamic_array3, DefaultIteratorMut, RawAccess, RawAccessMut,
-    RlstScalar,
+    rlst_dynamic_array2, rlst_dynamic_array3, DefaultIteratorMut, RandomAccessByRef, RawAccess,
+    RawAccessMut, RlstScalar, Shape,
 };
 
 pub struct PotentialCellAssembler<
@@ -15,6 +15,7 @@ pub struct PotentialCellAssembler<
     I: PotentialIntegrand<T = T>,
     G: GeometryMap<T = T::Real>,
     K: KernelEvaluator<T = T>,
+    Array2: RandomAccessByRef<2, Item = T::Real> + Shape<2> + RawAccess<Item = T::Real> + Sync,
 > {
     integrand: &'a I,
     kernel: &'a K,
@@ -25,7 +26,7 @@ pub struct PotentialCellAssembler<
     normals: RlstArray<T::Real, 2>,
     jacobians: RlstArray<T::Real, 2>,
     jdet: Vec<T::Real>,
-    evaluation_points: &'a RlstArray<T::Real, 2>,
+    evaluation_points: &'a Array2,
     weights: &'a [T::Real],
     cell: usize,
 }
@@ -36,7 +37,8 @@ impl<
         I: PotentialIntegrand<T = T>,
         G: GeometryMap<T = T::Real>,
         K: KernelEvaluator<T = T>,
-    > PotentialCellAssembler<'a, T, I, G, K>
+        Array2: RandomAccessByRef<2, Item = T::Real> + Shape<2> + RawAccess<Item = T::Real> + Sync,
+    > PotentialCellAssembler<'a, T, I, G, K, Array2>
 {
     #[allow(clippy::too_many_arguments)]
     pub fn new(
@@ -47,7 +49,7 @@ impl<
         kernel: &'a K,
         evaluator: G,
         table: &'a RlstArray<T, 4>,
-        evaluation_points: &'a RlstArray<T::Real, 2>,
+        evaluation_points: &'a Array2,
         weights: &'a [T::Real],
     ) -> Self {
         Self {
@@ -72,7 +74,8 @@ impl<
         I: PotentialIntegrand<T = T>,
         G: GeometryMap<T = T::Real>,
         K: KernelEvaluator<T = T>,
-    > CellAssembler for PotentialCellAssembler<'a, T, I, G, K>
+        Array2: RandomAccessByRef<2, Item = T::Real> + Shape<2> + RawAccess<Item = T::Real> + Sync,
+    > CellAssembler for PotentialCellAssembler<'a, T, I, G, K, Array2>
 {
     type T = T;
     fn set_cell(&mut self, cell: usize) {
diff --git a/src/assembly/potential/integrands/double_layer.rs b/src/assembly/potential/integrands/double_layer.rs
index 4d4adeac..4cd6a048 100644
--- a/src/assembly/potential/integrands/double_layer.rs
+++ b/src/assembly/potential/integrands/double_layer.rs
@@ -3,11 +3,13 @@ use crate::assembly::common::RlstArray;
 use crate::traits::{CellGeometry, PotentialIntegrand};
 use rlst::{RlstScalar, UnsafeRandomAccessByRef};
 
+/// Integrand for a double layer potential operator
 pub struct DoubleLayerPotentialIntegrand<T: RlstScalar> {
     _t: std::marker::PhantomData<T>,
 }
 
 impl<T: RlstScalar> DoubleLayerPotentialIntegrand<T> {
+    /// Create new
     pub fn new() -> Self {
         Self {
             _t: std::marker::PhantomData,
diff --git a/src/assembly/potential/integrands/single_layer.rs b/src/assembly/potential/integrands/single_layer.rs
index 7a1291f0..df6b2084 100644
--- a/src/assembly/potential/integrands/single_layer.rs
+++ b/src/assembly/potential/integrands/single_layer.rs
@@ -3,11 +3,13 @@ use crate::assembly::common::RlstArray;
 use crate::traits::{CellGeometry, PotentialIntegrand};
 use rlst::{RlstScalar, UnsafeRandomAccessByRef};
 
+/// Integrand for a single layer potential operator
 pub struct SingleLayerPotentialIntegrand<T: RlstScalar> {
     _t: std::marker::PhantomData<T>,
 }
 
 impl<T: RlstScalar> SingleLayerPotentialIntegrand<T> {
+    /// Create new
     pub fn new() -> Self {
         Self {
             _t: std::marker::PhantomData,
diff --git a/src/bindings.rs b/src/bindings.rs
index ab48e835..ebac62ef 100644
--- a/src/bindings.rs
+++ b/src/bindings.rs
@@ -669,7 +669,7 @@ pub mod function {
     }
 }
 
-pub mod assembly {
+pub mod boundary_assembly {
     use super::function::{extract_space, FunctionSpaceWrapper, GridType, SpaceType};
     use super::DType;
     use crate::assembly::common::SparseMatrixData;
@@ -4435,3 +4435,993 @@ pub mod assembly {
         }
     }
 }
+
+pub mod potential_assembly {
+    use super::boundary_assembly::KernelType;
+    use super::function::{extract_space, FunctionSpaceWrapper, GridType, SpaceType};
+    use super::DType;
+    use crate::{
+        assembly::kernels::KernelEvaluator,
+        assembly::potential::integrands::{
+            DoubleLayerPotentialIntegrand, SingleLayerPotentialIntegrand,
+        },
+        assembly::potential::PotentialAssembler,
+        function::SerialFunctionSpace,
+        traits::{
+            FunctionSpace, KernelEvaluator as KernelEvaluatorTrait, PotentialAssembly,
+            PotentialIntegrand,
+        },
+    };
+    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_slice2, rlst_array_from_slice_mut2, MatrixInverse, RlstScalar,
+    };
+    use std::ffi::c_void;
+    use std::slice::{from_raw_parts, from_raw_parts_mut};
+
+    #[derive(Debug, PartialEq, Clone, Copy)]
+    #[repr(u8)]
+    pub enum PotentialOperator {
+        SingleLayer = 0,
+        DoubleLayer = 1,
+        ElectricField = 2,
+        MagneticField = 3,
+    }
+
+    impl PotentialOperator {
+        fn from(value: u8) -> Option<Self> {
+            match value {
+                0 => Some(PotentialOperator::SingleLayer),
+                1 => Some(PotentialOperator::DoubleLayer),
+                2 => Some(PotentialOperator::ElectricField),
+                3 => Some(PotentialOperator::MagneticField),
+                _ => None,
+            }
+        }
+    }
+
+    #[repr(C)]
+    pub struct PotentialAssemblerWrapper {
+        pub assembler: *const c_void,
+        pub itype: PotentialOperator,
+        pub ktype: KernelType,
+        pub dtype: DType,
+    }
+    impl Drop for PotentialAssemblerWrapper {
+        fn drop(&mut self) {
+            let Self {
+                assembler,
+                itype,
+                ktype,
+                dtype,
+            } = self;
+            match ktype {
+                KernelType::Laplace => match itype {
+                    PotentialOperator::SingleLayer => match dtype {
+                        DType::F32 => drop(unsafe {
+                            Box::from_raw(
+                                *assembler
+                                    as *mut PotentialAssembler<
+                                        f32,
+                                        SingleLayerPotentialIntegrand<f32>,
+                                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                                    >,
+                            )
+                        }),
+                        DType::F64 => drop(unsafe {
+                            Box::from_raw(
+                                *assembler
+                                    as *mut PotentialAssembler<
+                                        f64,
+                                        SingleLayerPotentialIntegrand<f64>,
+                                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                                    >,
+                            )
+                        }),
+                        _ => {
+                            panic!("Invalid data type");
+                        }
+                    },
+                    PotentialOperator::DoubleLayer => match dtype {
+                        DType::F32 => drop(unsafe {
+                            Box::from_raw(
+                                *assembler
+                                    as *mut PotentialAssembler<
+                                        f32,
+                                        DoubleLayerPotentialIntegrand<f32>,
+                                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                                    >,
+                            )
+                        }),
+                        DType::F64 => drop(unsafe {
+                            Box::from_raw(
+                                *assembler
+                                    as *mut PotentialAssembler<
+                                        f64,
+                                        DoubleLayerPotentialIntegrand<f64>,
+                                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                                    >,
+                            )
+                        }),
+                        _ => {
+                            panic!("Invalid data type");
+                        }
+                    },
+                    _ => {
+                        panic!("Invalid operator");
+                    }
+                },
+                KernelType::Helmholtz => match itype {
+                    PotentialOperator::SingleLayer => match dtype {
+                        DType::C32 => drop(unsafe {
+                            Box::from_raw(
+                                *assembler
+                                    as *mut PotentialAssembler<
+                                        c32,
+                                        SingleLayerPotentialIntegrand<c32>,
+                                        KernelEvaluator<c32, Laplace3dKernel<c32>>,
+                                    >,
+                            )
+                        }),
+                        DType::C64 => drop(unsafe {
+                            Box::from_raw(
+                                *assembler
+                                    as *mut PotentialAssembler<
+                                        c64,
+                                        SingleLayerPotentialIntegrand<c64>,
+                                        KernelEvaluator<c64, Laplace3dKernel<c64>>,
+                                    >,
+                            )
+                        }),
+                        _ => {
+                            panic!("Invalid data type");
+                        }
+                    },
+                    PotentialOperator::DoubleLayer => match dtype {
+                        DType::C32 => drop(unsafe {
+                            Box::from_raw(
+                                *assembler
+                                    as *mut PotentialAssembler<
+                                        c32,
+                                        DoubleLayerPotentialIntegrand<c32>,
+                                        KernelEvaluator<c32, Laplace3dKernel<c32>>,
+                                    >,
+                            )
+                        }),
+                        DType::C64 => drop(unsafe {
+                            Box::from_raw(
+                                *assembler
+                                    as *mut PotentialAssembler<
+                                        c64,
+                                        DoubleLayerPotentialIntegrand<c64>,
+                                        KernelEvaluator<c64, Laplace3dKernel<c64>>,
+                                    >,
+                            )
+                        }),
+                        _ => {
+                            panic!("Invalid data type");
+                        }
+                    },
+                    _ => {
+                        panic!("Invalid operator");
+                    }
+                },
+            }
+        }
+    }
+
+    #[no_mangle]
+    pub unsafe extern "C" fn free_potential_assembler(a: *mut PotentialAssemblerWrapper) {
+        assert!(!a.is_null());
+        unsafe { drop(Box::from_raw(a)) }
+    }
+
+    pub(crate) unsafe fn extract_potential_assembler<
+        T: RlstScalar + MatrixInverse,
+        Integrand: PotentialIntegrand<T = T>,
+        Kernel: KernelEvaluatorTrait<T = T>,
+    >(
+        assembler: *const PotentialAssemblerWrapper,
+    ) -> *const PotentialAssembler<T, Integrand, Kernel> {
+        (*assembler).assembler as *const PotentialAssembler<T, Integrand, Kernel>
+    }
+
+    pub(crate) unsafe fn extract_potential_assembler_mut<
+        T: RlstScalar + MatrixInverse,
+        Integrand: PotentialIntegrand<T = T>,
+        Kernel: KernelEvaluatorTrait<T = T>,
+    >(
+        assembler: *const PotentialAssemblerWrapper,
+    ) -> *mut PotentialAssembler<T, Integrand, Kernel> {
+        (*assembler).assembler as *mut PotentialAssembler<T, Integrand, Kernel>
+    }
+
+    #[no_mangle]
+    pub unsafe extern "C" fn potential_assembler_has_quadrature_degree(
+        assembler: *mut PotentialAssemblerWrapper,
+        cell: u8,
+    ) -> bool {
+        let cell = ReferenceCellType::from(cell).unwrap();
+        match (*assembler).ktype {
+            KernelType::Laplace => match (*assembler).itype {
+                PotentialOperator::SingleLayer => match (*assembler).dtype {
+                    DType::F32 => (*extract_potential_assembler::<
+                        f32,
+                        SingleLayerPotentialIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    DType::F64 => (*extract_potential_assembler::<
+                        f64,
+                        SingleLayerPotentialIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                PotentialOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::F32 => (*extract_potential_assembler::<
+                        f32,
+                        DoubleLayerPotentialIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    DType::F64 => (*extract_potential_assembler::<
+                        f64,
+                        DoubleLayerPotentialIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+            KernelType::Helmholtz => match (*assembler).itype {
+                PotentialOperator::SingleLayer => match (*assembler).dtype {
+                    DType::C32 => (*extract_potential_assembler::<
+                        c32,
+                        SingleLayerPotentialIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    DType::C64 => (*extract_potential_assembler::<
+                        c64,
+                        SingleLayerPotentialIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                PotentialOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::C32 => (*extract_potential_assembler::<
+                        c32,
+                        DoubleLayerPotentialIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    DType::C64 => (*extract_potential_assembler::<
+                        c64,
+                        DoubleLayerPotentialIntegrand<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 potential_assembler_set_quadrature_degree(
+        assembler: *mut PotentialAssemblerWrapper,
+        cell: u8,
+        degree: usize,
+    ) {
+        let cell = ReferenceCellType::from(cell).unwrap();
+        match (*assembler).ktype {
+            KernelType::Laplace => match (*assembler).itype {
+                PotentialOperator::SingleLayer => match (*assembler).dtype {
+                    DType::F32 => (*extract_potential_assembler_mut::<
+                        f32,
+                        SingleLayerPotentialIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .set_quadrature_degree(cell, degree),
+                    DType::F64 => (*extract_potential_assembler_mut::<
+                        f64,
+                        SingleLayerPotentialIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .set_quadrature_degree(cell, degree),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                PotentialOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::F32 => (*extract_potential_assembler_mut::<
+                        f32,
+                        DoubleLayerPotentialIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .set_quadrature_degree(cell, degree),
+                    DType::F64 => (*extract_potential_assembler_mut::<
+                        f64,
+                        DoubleLayerPotentialIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .set_quadrature_degree(cell, degree),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+            KernelType::Helmholtz => match (*assembler).itype {
+                PotentialOperator::SingleLayer => match (*assembler).dtype {
+                    DType::C32 => (*extract_potential_assembler_mut::<
+                        c32,
+                        SingleLayerPotentialIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .set_quadrature_degree(cell, degree),
+                    DType::C64 => (*extract_potential_assembler_mut::<
+                        c64,
+                        SingleLayerPotentialIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .set_quadrature_degree(cell, degree),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                PotentialOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::C32 => (*extract_potential_assembler_mut::<
+                        c32,
+                        DoubleLayerPotentialIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .set_quadrature_degree(cell, degree),
+                    DType::C64 => (*extract_potential_assembler_mut::<
+                        c64,
+                        DoubleLayerPotentialIntegrand<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 potential_assembler_quadrature_degree(
+        assembler: *mut PotentialAssemblerWrapper,
+        cell: u8,
+    ) -> usize {
+        let cell = ReferenceCellType::from(cell).unwrap();
+        match (*assembler).ktype {
+            KernelType::Laplace => match (*assembler).itype {
+                PotentialOperator::SingleLayer => match (*assembler).dtype {
+                    DType::F32 => (*extract_potential_assembler::<
+                        f32,
+                        SingleLayerPotentialIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    DType::F64 => (*extract_potential_assembler::<
+                        f64,
+                        SingleLayerPotentialIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                PotentialOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::F32 => (*extract_potential_assembler::<
+                        f32,
+                        DoubleLayerPotentialIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    DType::F64 => (*extract_potential_assembler::<
+                        f64,
+                        DoubleLayerPotentialIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+            KernelType::Helmholtz => match (*assembler).itype {
+                PotentialOperator::SingleLayer => match (*assembler).dtype {
+                    DType::C32 => (*extract_potential_assembler::<
+                        c32,
+                        SingleLayerPotentialIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    DType::C64 => (*extract_potential_assembler::<
+                        c64,
+                        SingleLayerPotentialIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                PotentialOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::C32 => (*extract_potential_assembler::<
+                        c32,
+                        DoubleLayerPotentialIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    DType::C64 => (*extract_potential_assembler::<
+                        c64,
+                        DoubleLayerPotentialIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .quadrature_degree(cell),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+        }
+        .unwrap()
+    }
+
+    #[no_mangle]
+    pub unsafe extern "C" fn potential_assembler_set_batch_size(
+        assembler: *mut PotentialAssemblerWrapper,
+        batch_size: usize,
+    ) {
+        match (*assembler).ktype {
+            KernelType::Laplace => match (*assembler).itype {
+                PotentialOperator::SingleLayer => match (*assembler).dtype {
+                    DType::F32 => (*extract_potential_assembler_mut::<
+                        f32,
+                        SingleLayerPotentialIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .set_batch_size(batch_size),
+                    DType::F64 => (*extract_potential_assembler_mut::<
+                        f64,
+                        SingleLayerPotentialIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .set_batch_size(batch_size),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                PotentialOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::F32 => (*extract_potential_assembler_mut::<
+                        f32,
+                        DoubleLayerPotentialIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .set_batch_size(batch_size),
+                    DType::F64 => (*extract_potential_assembler_mut::<
+                        f64,
+                        DoubleLayerPotentialIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .set_batch_size(batch_size),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+            KernelType::Helmholtz => match (*assembler).itype {
+                PotentialOperator::SingleLayer => match (*assembler).dtype {
+                    DType::C32 => (*extract_potential_assembler_mut::<
+                        c32,
+                        SingleLayerPotentialIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .set_batch_size(batch_size),
+                    DType::C64 => (*extract_potential_assembler_mut::<
+                        c64,
+                        SingleLayerPotentialIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .set_batch_size(batch_size),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                PotentialOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::C32 => (*extract_potential_assembler_mut::<
+                        c32,
+                        DoubleLayerPotentialIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .set_batch_size(batch_size),
+                    DType::C64 => (*extract_potential_assembler_mut::<
+                        c64,
+                        DoubleLayerPotentialIntegrand<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 potential_assembler_batch_size(
+        assembler: *mut PotentialAssemblerWrapper,
+    ) -> usize {
+        match (*assembler).ktype {
+            KernelType::Laplace => match (*assembler).itype {
+                PotentialOperator::SingleLayer => match (*assembler).dtype {
+                    DType::F32 => (*extract_potential_assembler::<
+                        f32,
+                        SingleLayerPotentialIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .batch_size(),
+                    DType::F64 => (*extract_potential_assembler::<
+                        f64,
+                        SingleLayerPotentialIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .batch_size(),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                PotentialOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::F32 => (*extract_potential_assembler::<
+                        f32,
+                        DoubleLayerPotentialIntegrand<f32>,
+                        KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                    >(assembler))
+                    .batch_size(),
+                    DType::F64 => (*extract_potential_assembler::<
+                        f64,
+                        DoubleLayerPotentialIntegrand<f64>,
+                        KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                    >(assembler))
+                    .batch_size(),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+            KernelType::Helmholtz => match (*assembler).itype {
+                PotentialOperator::SingleLayer => match (*assembler).dtype {
+                    DType::C32 => (*extract_potential_assembler::<
+                        c32,
+                        SingleLayerPotentialIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .batch_size(),
+                    DType::C64 => (*extract_potential_assembler::<
+                        c64,
+                        SingleLayerPotentialIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .batch_size(),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                PotentialOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::C32 => (*extract_potential_assembler::<
+                        c32,
+                        DoubleLayerPotentialIntegrand<c32>,
+                        KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                    >(assembler))
+                    .batch_size(),
+                    DType::C64 => (*extract_potential_assembler::<
+                        c64,
+                        DoubleLayerPotentialIntegrand<c64>,
+                        KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                    >(assembler))
+                    .batch_size(),
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+        }
+    }
+
+    unsafe fn potential_assembler_assemble_into_dense_internal_real<
+        T: RlstScalar<Real = T> + MatrixInverse,
+        Integrand: PotentialIntegrand<T = T>,
+        Kernel: KernelEvaluatorTrait<T = T>,
+        Space: FunctionSpace<T = T> + Sync,
+    >(
+        assembler: *mut PotentialAssemblerWrapper,
+        output: *mut c_void,
+        space: *const FunctionSpaceWrapper,
+        points: *const c_void,
+        npts: usize,
+    ) {
+        let points = rlst_array_from_slice2!(
+            from_raw_parts(points as *const T::Real, 3 * npts),
+            [3, npts]
+        );
+        let dim = (*extract_space::<Space>(space)).global_size();
+        let mut output = rlst_array_from_slice_mut2!(
+            from_raw_parts_mut(output as *mut T, npts * dim),
+            [npts, dim]
+        );
+
+        (*extract_potential_assembler::<T, Integrand, Kernel>(assembler)).assemble_into_dense(
+            &mut output,
+            &*extract_space::<Space>(space),
+            &points,
+        )
+    }
+    unsafe fn potential_assembler_assemble_into_dense_internal_complex<
+        T: RlstScalar<Complex = T> + MatrixInverse,
+        Integrand: PotentialIntegrand<T = T>,
+        Kernel: KernelEvaluatorTrait<T = T>,
+        Space: FunctionSpace<T = T> + Sync,
+    >(
+        assembler: *mut PotentialAssemblerWrapper,
+        output: *mut c_void,
+        space: *const FunctionSpaceWrapper,
+        points: *const c_void,
+        npts: usize,
+    ) {
+        let points = rlst_array_from_slice2!(
+            from_raw_parts(points as *const T::Real, 3 * npts),
+            [3, npts]
+        );
+        let dim = (*extract_space::<Space>(space)).global_size();
+        let mut output = rlst_array_from_slice_mut2!(
+            from_raw_parts_mut(output as *mut T, npts * dim),
+            [npts, dim]
+        );
+
+        (*extract_potential_assembler::<T, Integrand, Kernel>(assembler)).assemble_into_dense(
+            &mut output,
+            &*extract_space::<Space>(space),
+            &points,
+        )
+    }
+    #[no_mangle]
+    pub unsafe extern "C" fn potential_assembler_assemble_into_dense(
+        assembler: *mut PotentialAssemblerWrapper,
+        output: *mut c_void,
+        space: *const FunctionSpaceWrapper,
+        points: *const c_void,
+        npts: usize,
+    ) {
+        match (*assembler).dtype {
+            DType::F32 => {
+                assert_eq!((*space).dtype, DType::F32);
+            }
+            DType::F64 => {
+                assert_eq!((*space).dtype, DType::F64);
+            }
+            DType::C32 => {
+                assert_eq!((*space).dtype, DType::F32);
+            }
+            DType::C64 => {
+                assert_eq!((*space).dtype, DType::F64);
+            }
+        }
+        match (*assembler).ktype {
+            KernelType::Laplace => match (*assembler).itype {
+                PotentialOperator::SingleLayer => match (*assembler).dtype {
+                    DType::F32 => match (*space).stype {
+                        SpaceType::SerialFunctionSpace => match (*space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                potential_assembler_assemble_into_dense_internal_real::<
+                                    f32,
+                                    SingleLayerPotentialIntegrand<f32>,
+                                    KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                                    SerialFunctionSpace<
+                                        f32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(
+                                    assembler, output, space, points, npts
+                                );
+                            }
+                        },
+                    },
+                    DType::F64 => match (*space).stype {
+                        SpaceType::SerialFunctionSpace => match (*space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                potential_assembler_assemble_into_dense_internal_real::<
+                                    f64,
+                                    SingleLayerPotentialIntegrand<f64>,
+                                    KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                                    SerialFunctionSpace<
+                                        f64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(
+                                    assembler, output, space, points, npts
+                                );
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                PotentialOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::F32 => match (*space).stype {
+                        SpaceType::SerialFunctionSpace => match (*space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                potential_assembler_assemble_into_dense_internal_real::<
+                                    f32,
+                                    DoubleLayerPotentialIntegrand<f32>,
+                                    KernelEvaluator<f32, Laplace3dKernel<f32>>,
+                                    SerialFunctionSpace<
+                                        f32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(
+                                    assembler, output, space, points, npts
+                                );
+                            }
+                        },
+                    },
+                    DType::F64 => match (*space).stype {
+                        SpaceType::SerialFunctionSpace => match (*space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                potential_assembler_assemble_into_dense_internal_real::<
+                                    f64,
+                                    DoubleLayerPotentialIntegrand<f64>,
+                                    KernelEvaluator<f64, Laplace3dKernel<f64>>,
+                                    SerialFunctionSpace<
+                                        f64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(
+                                    assembler, output, space, points, npts
+                                );
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+            KernelType::Helmholtz => match (*assembler).itype {
+                PotentialOperator::SingleLayer => match (*assembler).dtype {
+                    DType::C32 => match (*space).stype {
+                        SpaceType::SerialFunctionSpace => match (*space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                potential_assembler_assemble_into_dense_internal_complex::<
+                                    c32,
+                                    SingleLayerPotentialIntegrand<c32>,
+                                    KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                                    SerialFunctionSpace<
+                                        c32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(
+                                    assembler, output, space, points, npts
+                                );
+                            }
+                        },
+                    },
+                    DType::C64 => match (*space).stype {
+                        SpaceType::SerialFunctionSpace => match (*space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                potential_assembler_assemble_into_dense_internal_complex::<
+                                    c64,
+                                    SingleLayerPotentialIntegrand<c64>,
+                                    KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                                    SerialFunctionSpace<
+                                        c64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(
+                                    assembler, output, space, points, npts
+                                );
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                PotentialOperator::DoubleLayer => match (*assembler).dtype {
+                    DType::C32 => match (*space).stype {
+                        SpaceType::SerialFunctionSpace => match (*space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                potential_assembler_assemble_into_dense_internal_complex::<
+                                    c32,
+                                    DoubleLayerPotentialIntegrand<c32>,
+                                    KernelEvaluator<c32, Helmholtz3dKernel<c32>>,
+                                    SerialFunctionSpace<
+                                        c32,
+                                        SingleElementGrid<f32, CiarletElement<f32>>,
+                                    >,
+                                >(
+                                    assembler, output, space, points, npts
+                                );
+                            }
+                        },
+                    },
+                    DType::C64 => match (*space).stype {
+                        SpaceType::SerialFunctionSpace => match (*space).gtype {
+                            GridType::SerialSingleElementGrid => {
+                                potential_assembler_assemble_into_dense_internal_complex::<
+                                    c64,
+                                    DoubleLayerPotentialIntegrand<c64>,
+                                    KernelEvaluator<c64, Helmholtz3dKernel<c64>>,
+                                    SerialFunctionSpace<
+                                        c64,
+                                        SingleElementGrid<f64, CiarletElement<f64>>,
+                                    >,
+                                >(
+                                    assembler, output, space, points, npts
+                                );
+                            }
+                        },
+                    },
+                    _ => {
+                        panic!("Invalid data type");
+                    }
+                },
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+        }
+    }
+
+    #[no_mangle]
+    pub unsafe extern "C" fn potential_assembler_dtype(
+        assembler: *mut PotentialAssemblerWrapper,
+    ) -> u8 {
+        (*assembler).dtype as u8
+    }
+
+    unsafe fn laplace_potential_assembler_new_internal<T: RlstScalar + MatrixInverse>(
+        operator: PotentialOperator,
+        dtype: DType,
+    ) -> *const PotentialAssemblerWrapper {
+        Box::into_raw(Box::new(PotentialAssemblerWrapper {
+            assembler: match operator {
+                PotentialOperator::SingleLayer => {
+                    Box::into_raw(Box::new(
+                        PotentialAssembler::<T, _, _>::new_laplace_single_layer(),
+                    )) as *const c_void
+                }
+                PotentialOperator::DoubleLayer => {
+                    Box::into_raw(Box::new(
+                        PotentialAssembler::<T, _, _>::new_laplace_double_layer(),
+                    )) as *const c_void
+                }
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+            itype: operator,
+            ktype: KernelType::Laplace,
+            dtype,
+        }))
+    }
+
+    #[no_mangle]
+    pub unsafe extern "C" fn laplace_potential_assembler_new(
+        operator: u8,
+        dtype: u8,
+    ) -> *const PotentialAssemblerWrapper {
+        let operator = PotentialOperator::from(operator).unwrap();
+        let dtype = DType::from(dtype).unwrap();
+        match dtype {
+            DType::F32 => laplace_potential_assembler_new_internal::<f32>(operator, dtype),
+            DType::F64 => laplace_potential_assembler_new_internal::<f64>(operator, dtype),
+            _ => {
+                panic!("Invalid data type");
+            }
+        }
+    }
+
+    unsafe fn helmholtz_potential_assembler_new_internal<
+        T: RlstScalar<Complex = T> + MatrixInverse,
+    >(
+        wavenumber: T::Real,
+        operator: PotentialOperator,
+        dtype: DType,
+    ) -> *const PotentialAssemblerWrapper {
+        Box::into_raw(Box::new(PotentialAssemblerWrapper {
+            assembler: match operator {
+                PotentialOperator::SingleLayer => Box::into_raw(Box::new(
+                    PotentialAssembler::<T, _, _>::new_helmholtz_single_layer(wavenumber),
+                )) as *const c_void,
+                PotentialOperator::DoubleLayer => Box::into_raw(Box::new(
+                    PotentialAssembler::<T, _, _>::new_helmholtz_double_layer(wavenumber),
+                )) as *const c_void,
+                _ => {
+                    panic!("Invalid operator");
+                }
+            },
+            itype: operator,
+            ktype: KernelType::Helmholtz,
+            dtype,
+        }))
+    }
+
+    #[no_mangle]
+    pub unsafe extern "C" fn helmholtz_potential_assembler_new(
+        wavenumber: *const c_void,
+        operator: u8,
+        dtype: u8,
+    ) -> *const PotentialAssemblerWrapper {
+        let operator = PotentialOperator::from(operator).unwrap();
+        let dtype = DType::from(dtype).unwrap();
+        match dtype {
+            DType::C32 => helmholtz_potential_assembler_new_internal::<c32>(
+                *(wavenumber as *const f32),
+                operator,
+                dtype,
+            ),
+            DType::C64 => helmholtz_potential_assembler_new_internal::<c64>(
+                *(wavenumber as *const f64),
+                operator,
+                dtype,
+            ),
+            _ => {
+                panic!("Invalid data type");
+            }
+        }
+    }
+}
diff --git a/src/traits/assembly/potential.rs b/src/traits/assembly/potential.rs
index b06b6e22..428a3f8b 100644
--- a/src/traits/assembly/potential.rs
+++ b/src/traits/assembly/potential.rs
@@ -2,7 +2,7 @@
 use super::CellGeometry;
 use crate::assembly::common::RlstArray;
 use crate::traits::FunctionSpace;
-use rlst::RlstScalar;
+use rlst::{RandomAccessByRef, RandomAccessMut, RawAccess, RawAccessMut, RlstScalar, Shape};
 
 pub unsafe trait PotentialIntegrand {
     //! Integrand
@@ -42,10 +42,17 @@ pub trait PotentialAssembly {
     type T: RlstScalar;
 
     /// 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,
+        Array2Mut: RandomAccessMut<2, Item = Self::T> + Shape<2> + RawAccessMut<Item = Self::T>,
+        Array2: RandomAccessByRef<2, Item = <Self::T as RlstScalar>::Real>
+            + Shape<2>
+            + RawAccess<Item = <Self::T as RlstScalar>::Real>
+            + Sync,
+    >(
         &self,
-        output: &mut RlstArray<Self::T, 2>,
+        output: &mut Array2Mut,
         space: &Space,
-        points: &RlstArray<<Self::T as RlstScalar>::Real, 2>,
+        points: &Array2,
     );
 }


Test patch:
diff --git a/python/test/test_boundary_assembler.py b/python/test/test_boundary_assembler.py
index 5e75bede..a2b715da 100644
--- a/python/test/test_boundary_assembler.py
+++ b/python/test/test_boundary_assembler.py
@@ -40,103 +40,6 @@ def test_create_assembler(otype):
     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",
     [
@@ -270,3 +173,100 @@ def test_assemble_singular_and_nonsingular(
 
     mat2 = a.assemble_into_dense(trial_space, test_space)
     assert np.allclose(mat, mat2)
+
+
+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)
diff --git a/python/test/test_potential_assembler.py b/python/test/test_potential_assembler.py
new file mode 100644
index 00000000..1c7de0e7
--- /dev/null
+++ b/python/test/test_potential_assembler.py
@@ -0,0 +1,81 @@
+import pytest
+import numpy as np
+from ndelement.reference_cell import ReferenceCellType
+from bempp.assembly.potential 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,
+    ],
+)
+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.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)
+
+    points = np.array([[2.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 2.0]])
+
+    mat = a.assemble_into_dense(space, points)
+
+    from_cl = np.array(
+        [
+            [
+                0.04038047926587569,
+                0.02879904511649957,
+                0.02879904511649957,
+                0.0403804792658757,
+                0.04038047926587569,
+                0.028799045116499562,
+                0.02879904511649957,
+                0.04038047926587571,
+            ],
+            [
+                0.0403804792658757,
+                0.04038047926587569,
+                0.028799045116499573,
+                0.02879904511649957,
+                0.04038047926587571,
+                0.04038047926587569,
+                0.028799045116499573,
+                0.028799045116499573,
+            ],
+            [
+                0.04038047926587571,
+                0.04038047926587571,
+                0.04038047926587571,
+                0.04038047926587571,
+                0.028799045116499573,
+                0.028799045116499573,
+                0.028799045116499573,
+                0.028799045116499573,
+            ],
+        ]
+    )
+
+    assert np.allclose(mat, from_cl, rtol=1e-4)