start adding potential assemblers

python/bempp/assembly/potential.py

@@ -0,0 +1,89 @@
+"""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
+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")
+
+
+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 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))
+    )

python/test/test_potential_assembler.py

@@ -0,0 +1,32 @@
+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

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;

src/assembly/potential/assemblers.rs

@@ -135,15 +135,25 @@ 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;
     }
 
+    /// 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>(
         &self,
         output: &RawData2D<T>,

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,

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,