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

Issue:
Refactor boundary assembly into better traits
New trait structure (defined in src/traits/assembly.rs):
- `KernelEvaluator` evaluates Green's kernels with the correct EvalType
- `BoundaryIntegrand` evaluates the integrand for a given test and trial point on a given cell
- `CellPairAssembler` assembles the contribution to the matrix due to a pair of cells
- `BoundaryAssembler` and `ParallelBoundaryAssembler` contain user-facing full assembly functions


Discussion / hints:

Fix patch:
diff --git a/benches/assembly_benchmark.rs b/benches/assembly_benchmark.rs
index c7adac70..26ea9d2c 100644
--- a/benches/assembly_benchmark.rs
+++ b/benches/assembly_benchmark.rs
@@ -1,4 +1,4 @@
-use bempp::assembly::{boundary, boundary::BoundaryAssembler};
+use bempp::assembly::boundary::BoundaryAssembler;
 use bempp::function::SerialFunctionSpace;
 use bempp::traits::{BoundaryAssembly, FunctionSpace};
 use criterion::{criterion_group, criterion_main, Criterion};
@@ -19,7 +19,7 @@ pub fn assembly_parts_benchmark(c: &mut Criterion) {
         let mut matrix = rlst_dynamic_array2!(f64, [space.global_size(), space.global_size()]);
 
         let colouring = space.cell_colouring();
-        let mut a = boundary::SingleLayerAssembler::<f64, _>::new_laplace();
+        let mut a = BoundaryAssembler::<f64, _, _>::new_laplace_single_layer();
         a.quadrature_degree(ReferenceCellType::Triangle, 16);
         a.singular_quadrature_degree(
             (ReferenceCellType::Triangle, ReferenceCellType::Triangle),
diff --git a/examples/assembly.rs b/examples/assembly.rs
index 570a4088..586eaa09 100644
--- a/examples/assembly.rs
+++ b/examples/assembly.rs
@@ -1,4 +1,4 @@
-use bempp::assembly::{boundary, boundary::BoundaryAssembler};
+use bempp::assembly::boundary::BoundaryAssembler;
 use bempp::function::SerialFunctionSpace;
 use bempp::traits::{BoundaryAssembly, FunctionSpace};
 use ndelement::ciarlet::LagrangeElementFamily;
@@ -17,7 +17,7 @@ fn main() {
     let mut matrix = rlst_dynamic_array2!(f64, [ndofs, ndofs]);
 
     // Create an assembler for the Laplace single layer operator
-    let mut a = boundary::SingleLayerAssembler::<f64, _>::new_laplace();
+    let mut a = BoundaryAssembler::<f64, _, _>::new_laplace_single_layer();
 
     // Adjust the quadrature degree for non-singular integrals on a triangle.
     // This makes the integrals use a quadrature rule with 16 points
diff --git a/src/assembly.rs b/src/assembly.rs
index 7570c7ad..9e39be3a 100644
--- a/src/assembly.rs
+++ b/src/assembly.rs
@@ -133,12 +133,12 @@ mod test {
         };
         (Laplace, $operator:ident, $dtype:ident) => {
             paste! {
-                boundary::[<$operator Assembler>]::<[<$dtype>], _>::new_laplace()
+                boundary::BoundaryAssembler::<[<$dtype>], _, _>::[<new_laplace_ $operator>]()
             }
         };
         (Helmholtz, $operator:ident, $dtype:ident) => {
             paste! {
-                boundary::[<$operator Assembler>]::<[<$dtype>], _>::new_helmholtz(3.0)
+                boundary::BoundaryAssembler::<[<$dtype>], _, _>::[<new_helmholtz_ $operator>](3.0)
             }
         };
     }
@@ -169,76 +169,76 @@ mod test {
     }
 
     test_assembly!(
-        (f64, Laplace, SingleLayer, Triangle),
-        (f32, Laplace, SingleLayer, Triangle),
-        //(c64, Laplace, SingleLayer, Triangle),
-        //(c32, Laplace, SingleLayer, Triangle),
-        (f64, Laplace, DoubleLayer, Triangle),
-        (f32, Laplace, DoubleLayer, Triangle),
-        //(c64, Laplace, DoubleLayer, Triangle),
-        //(c32, Laplace, DoubleLayer, Triangle),
-        (f64, Laplace, AdjointDoubleLayer, Triangle),
-        (f32, Laplace, AdjointDoubleLayer, Triangle),
-        //(c64, Laplace, AdjointDoubleLayer, Triangle),
-        //(c32, Laplace, AdjointDoubleLayer, Triangle),
-        (f64, Laplace, Hypersingular, Triangle),
-        (f32, Laplace, Hypersingular, Triangle),
-        //(c64, Laplace, Hypersingular, Triangle),
-        //(c32, Laplace, Hypersingular, Triangle),
-        (c64, Helmholtz, SingleLayer, Triangle),
-        (c32, Helmholtz, SingleLayer, Triangle),
-        (c64, Helmholtz, DoubleLayer, Triangle),
-        (c32, Helmholtz, DoubleLayer, Triangle),
-        (c64, Helmholtz, AdjointDoubleLayer, Triangle),
-        (c32, Helmholtz, AdjointDoubleLayer, Triangle),
-        (c64, Helmholtz, Hypersingular, Triangle),
-        (c32, Helmholtz, Hypersingular, Triangle),
-        (f64, Laplace, SingleLayer, Quadrilateral),
-        (f32, Laplace, SingleLayer, Quadrilateral),
-        //(c64, Laplace, SingleLayer, Quadrilateral),
-        //(c32, Laplace, SingleLayer, Quadrilateral),
-        (f64, Laplace, DoubleLayer, Quadrilateral),
-        (f32, Laplace, DoubleLayer, Quadrilateral),
-        //(c64, Laplace, DoubleLayer, Quadrilateral),
-        //(c32, Laplace, DoubleLayer, Quadrilateral),
-        (f64, Laplace, AdjointDoubleLayer, Quadrilateral),
-        (f32, Laplace, AdjointDoubleLayer, Quadrilateral),
-        //(c64, Laplace, AdjointDoubleLayer, Quadrilateral),
-        //(c32, Laplace, AdjointDoubleLayer, Quadrilateral),
-        (f64, Laplace, Hypersingular, Quadrilateral),
-        (f32, Laplace, Hypersingular, Quadrilateral),
-        //(c64, Laplace, Hypersingular, Quadrilateral),
-        //(c32, Laplace, Hypersingular, Quadrilateral),
-        (c64, Helmholtz, SingleLayer, Quadrilateral),
-        (c32, Helmholtz, SingleLayer, Quadrilateral),
-        (c64, Helmholtz, DoubleLayer, Quadrilateral),
-        (c32, Helmholtz, DoubleLayer, Quadrilateral),
-        (c64, Helmholtz, AdjointDoubleLayer, Quadrilateral),
-        (c32, Helmholtz, AdjointDoubleLayer, Quadrilateral),
-        (c64, Helmholtz, Hypersingular, Quadrilateral),
-        (c32, Helmholtz, Hypersingular, Quadrilateral) //(f64, Laplace, SingleLayer, Mixed),
-                                                       //(f32, Laplace, SingleLayer, Mixed),
-                                                       //(c64, Laplace, SingleLayer, Mixed),
-                                                       //(c32, Laplace, SingleLayer, Mixed),
-                                                       //(f64, Laplace, DoubleLayer, Mixed),
-                                                       //(f32, Laplace, DoubleLayer, Mixed),
-                                                       //(c64, Laplace, DoubleLayer, Mixed),
-                                                       //(c32, Laplace, DoubleLayer, Mixed),
-                                                       //(f64, Laplace, AdjointDoubleLayer, Mixed),
-                                                       //(f32, Laplace, AdjointDoubleLayer, Mixed),
-                                                       //(c64, Laplace, AdjointDoubleLayer, Mixed),
-                                                       //(c32, Laplace, AdjointDoubleLayer, Mixed),
-                                                       //(f64, Laplace, Hypersingular, Mixed),
-                                                       //(f32, Laplace, Hypersingular, Mixed),
-                                                       //(c64, Laplace, Hypersingular, Mixed),
-                                                       //(c32, Laplace, Hypersingular, Mixed),
-                                                       //(c64, Helmholtz, SingleLayer, Mixed),
-                                                       //(c32, Helmholtz, SingleLayer, Mixed),
-                                                       //(c64, Helmholtz, DoubleLayer, Mixed),
-                                                       //(c32, Helmholtz, DoubleLayer, Mixed),
-                                                       //(c64, Helmholtz, AdjointDoubleLayer, Mixed),
-                                                       //(c32, Helmholtz, AdjointDoubleLayer, Mixed),
-                                                       //(c64, Helmholtz, Hypersingular, Mixed),
-                                                       //(c32, Helmholtz, Hypersingular, Mixed)
+        (f64, Laplace, single_layer, Triangle),
+        (f32, Laplace, single_layer, Triangle),
+        //(c64, Laplace, single_layer, Triangle),
+        //(c32, Laplace, single_layer, Triangle),
+        (f64, Laplace, double_layer, Triangle),
+        (f32, Laplace, double_layer, Triangle),
+        //(c64, Laplace, double_layer, Triangle),
+        //(c32, Laplace, double_layer, Triangle),
+        (f64, Laplace, adjoint_double_layer, Triangle),
+        (f32, Laplace, adjoint_double_layer, Triangle),
+        //(c64, Laplace, adjoint_double_layer, Triangle),
+        //(c32, Laplace, adjoint_double_layer, Triangle),
+        (f64, Laplace, hypersingular, Triangle),
+        (f32, Laplace, hypersingular, Triangle),
+        //(c64, Laplace, hypersingular, Triangle),
+        //(c32, Laplace, hypersingular, Triangle),
+        (c64, Helmholtz, single_layer, Triangle),
+        (c32, Helmholtz, single_layer, Triangle),
+        (c64, Helmholtz, double_layer, Triangle),
+        (c32, Helmholtz, double_layer, Triangle),
+        (c64, Helmholtz, adjoint_double_layer, Triangle),
+        (c32, Helmholtz, adjoint_double_layer, Triangle),
+        (c64, Helmholtz, hypersingular, Triangle),
+        (c32, Helmholtz, hypersingular, Triangle),
+        (f64, Laplace, single_layer, Quadrilateral),
+        (f32, Laplace, single_layer, Quadrilateral),
+        //(c64, Laplace, single_layer, Quadrilateral),
+        //(c32, Laplace, single_layer, Quadrilateral),
+        (f64, Laplace, double_layer, Quadrilateral),
+        (f32, Laplace, double_layer, Quadrilateral),
+        //(c64, Laplace, double_layer, Quadrilateral),
+        //(c32, Laplace, double_layer, Quadrilateral),
+        (f64, Laplace, adjoint_double_layer, Quadrilateral),
+        (f32, Laplace, adjoint_double_layer, Quadrilateral),
+        //(c64, Laplace, adjoint_double_layer, Quadrilateral),
+        //(c32, Laplace, adjoint_double_layer, Quadrilateral),
+        (f64, Laplace, hypersingular, Quadrilateral),
+        (f32, Laplace, hypersingular, Quadrilateral),
+        //(c64, Laplace, hypersingular, Quadrilateral),
+        //(c32, Laplace, hypersingular, Quadrilateral),
+        (c64, Helmholtz, single_layer, Quadrilateral),
+        (c32, Helmholtz, single_layer, Quadrilateral),
+        (c64, Helmholtz, double_layer, Quadrilateral),
+        (c32, Helmholtz, double_layer, Quadrilateral),
+        (c64, Helmholtz, adjoint_double_layer, Quadrilateral),
+        (c32, Helmholtz, adjoint_double_layer, Quadrilateral),
+        (c64, Helmholtz, hypersingular, Quadrilateral),
+        (c32, Helmholtz, hypersingular, Quadrilateral) //(f64, Laplace, single_layer, Mixed),
+                                                       //(f32, Laplace, single_layer, Mixed),
+                                                       //(c64, Laplace, single_layer, Mixed),
+                                                       //(c32, Laplace, single_layer, Mixed),
+                                                       //(f64, Laplace, double_layer, Mixed),
+                                                       //(f32, Laplace, double_layer, Mixed),
+                                                       //(c64, Laplace, double_layer, Mixed),
+                                                       //(c32, Laplace, double_layer, Mixed),
+                                                       //(f64, Laplace, adjoint_double_layer, Mixed),
+                                                       //(f32, Laplace, adjoint_double_layer, Mixed),
+                                                       //(c64, Laplace, adjoint_double_layer, Mixed),
+                                                       //(c32, Laplace, adjoint_double_layer, Mixed),
+                                                       //(f64, Laplace, hypersingular, Mixed),
+                                                       //(f32, Laplace, hypersingular, Mixed),
+                                                       //(c64, Laplace, hypersingular, Mixed),
+                                                       //(c32, Laplace, hypersingular, Mixed),
+                                                       //(c64, Helmholtz, single_layer, Mixed),
+                                                       //(c32, Helmholtz, single_layer, Mixed),
+                                                       //(c64, Helmholtz, double_layer, Mixed),
+                                                       //(c32, Helmholtz, double_layer, Mixed),
+                                                       //(c64, Helmholtz, adjoint_double_layer, Mixed),
+                                                       //(c32, Helmholtz, adjoint_double_layer, Mixed),
+                                                       //(c64, Helmholtz, hypersingular, Mixed),
+                                                       //(c32, Helmholtz, hypersingular, Mixed)
     );
 }
diff --git a/src/assembly/boundary.rs b/src/assembly/boundary.rs
index ba7668d3..6162d01b 100644
--- a/src/assembly/boundary.rs
+++ b/src/assembly/boundary.rs
@@ -4,9 +4,6 @@ pub(crate) mod cell_pair_assemblers;
 mod integrands;
 
 pub use assemblers::BoundaryAssembler;
-pub use assemblers::{
-    AdjointDoubleLayerAssembler, DoubleLayerAssembler, HypersingularAssembler, SingleLayerAssembler,
-};
 
 #[cfg(test)]
 mod test {
@@ -29,7 +26,7 @@ mod test {
         let ndofs = space.global_size();
 
         let mut matrix = rlst_dynamic_array2!(f64, [ndofs, ndofs]);
-        let assembler = SingleLayerAssembler::<f64, _>::new_laplace();
+        let assembler = BoundaryAssembler::<f64, _, _>::new_laplace_single_layer();
         assembler.assemble_singular_into_dense(&mut matrix, &space, &space);
         let csr = assembler.assemble_singular_into_csr(&space, &space);
 
@@ -55,7 +52,7 @@ mod test {
         let ndofs = space.global_size();
 
         let mut matrix = rlst_dynamic_array2!(f64, [ndofs, ndofs]);
-        let assembler = SingleLayerAssembler::<f64, _>::new_laplace();
+        let assembler = BoundaryAssembler::<f64, _, _>::new_laplace_single_layer();
         assembler.assemble_singular_into_dense(&mut matrix, &space, &space);
         let csr = assembler.assemble_singular_into_csr(&space, &space);
 
@@ -84,7 +81,7 @@ mod test {
         let ndofs1 = space1.global_size();
 
         let mut matrix = rlst_dynamic_array2!(f64, [ndofs1, ndofs0]);
-        let assembler = SingleLayerAssembler::<f64, _>::new_laplace();
+        let assembler = BoundaryAssembler::<f64, _, _>::new_laplace_single_layer();
         assembler.assemble_singular_into_dense(&mut matrix, &space0, &space1);
         let csr = assembler.assemble_singular_into_csr(&space0, &space1);
         let indptr = csr.indptr();
diff --git a/src/assembly/boundary/assemblers.rs b/src/assembly/boundary/assemblers.rs
index f9ffe15d..99ada30f 100644
--- a/src/assembly/boundary/assemblers.rs
+++ b/src/assembly/boundary/assemblers.rs
@@ -3,10 +3,6 @@ pub mod adjoint_double_layer;
 pub mod double_layer;
 pub mod hypersingular;
 pub mod single_layer;
-pub use adjoint_double_layer::AdjointDoubleLayerAssembler;
-pub use double_layer::DoubleLayerAssembler;
-pub use hypersingular::HypersingularAssembler;
-pub use single_layer::SingleLayerAssembler;
 
 use super::cell_pair_assemblers::{NonsingularCellPairAssembler, SingularCellPairAssembler};
 use crate::assembly::common::{equal_grids, RawData2D, RlstArray, SparseMatrixData};
@@ -151,8 +147,13 @@ where
 
 /// Assemble the contribution to the terms of a matrix for a batch of pairs of adjacent cells
 #[allow(clippy::too_many_arguments)]
-fn assemble_batch_singular<T: RlstScalar + MatrixInverse, Space: FunctionSpace<T = T>>(
-    assembler: &impl BoundaryAssembler<T = T>,
+fn assemble_batch_singular<
+    T: RlstScalar + MatrixInverse,
+    Space: FunctionSpace<T = T>,
+    Integrand: BoundaryIntegrand<T = T>,
+    Kernel: KernelEvaluator<T = T>,
+>(
+    assembler: &BoundaryAssembler<T, Integrand, Kernel>,
     deriv_size: usize,
     shape: [usize; 2],
     trial_cell_type: ReferenceCellType,
@@ -187,8 +188,8 @@ fn assemble_batch_singular<T: RlstScalar + MatrixInverse, Space: FunctionSpace<T
     let mut a = SingularCellPairAssembler::new(
         npts,
         deriv_size,
-        assembler.integrand(),
-        assembler.kernel(),
+        &assembler.integrand,
+        &assembler.kernel,
         test_evaluator,
         trial_evaluator,
         test_table,
@@ -225,8 +226,13 @@ fn assemble_batch_singular<T: RlstScalar + MatrixInverse, Space: FunctionSpace<T
 
 /// Assemble the contribution to the terms of a matrix for a batch of non-adjacent cells
 #[allow(clippy::too_many_arguments)]
-fn assemble_batch_nonadjacent<T: RlstScalar + MatrixInverse, Space: FunctionSpace<T = T>>(
-    assembler: &impl BoundaryAssembler<T = T>,
+fn assemble_batch_nonadjacent<
+    T: RlstScalar + MatrixInverse,
+    Space: FunctionSpace<T = T>,
+    Integrand: BoundaryIntegrand<T = T>,
+    Kernel: KernelEvaluator<T = T>,
+>(
+    assembler: &BoundaryAssembler<T, Integrand, Kernel>,
     deriv_size: usize,
     output: &RawData2D<T>,
     trial_cell_type: ReferenceCellType,
@@ -262,8 +268,8 @@ fn assemble_batch_nonadjacent<T: RlstScalar + MatrixInverse, Space: FunctionSpac
         npts_test,
         npts_trial,
         deriv_size,
-        assembler.integrand(),
-        assembler.kernel(),
+        &assembler.integrand,
+        &assembler.kernel,
         test_evaluator,
         trial_evaluator,
         test_table,
@@ -310,8 +316,10 @@ fn assemble_batch_nonadjacent<T: RlstScalar + MatrixInverse, Space: FunctionSpac
 fn assemble_batch_singular_correction<
     T: RlstScalar + MatrixInverse,
     Space: FunctionSpace<T = T>,
+    Integrand: BoundaryIntegrand<T = T>,
+    Kernel: KernelEvaluator<T = T>,
 >(
-    assembler: &impl BoundaryAssembler<T = T>,
+    assembler: &BoundaryAssembler<T, Integrand, Kernel>,
     deriv_size: usize,
     shape: [usize; 2],
     trial_cell_type: ReferenceCellType,
@@ -348,8 +356,8 @@ fn assemble_batch_singular_correction<
         npts_test,
         npts_trial,
         deriv_size,
-        assembler.integrand(),
-        assembler.kernel(),
+        &assembler.integrand,
+        &assembler.kernel,
         test_evaluator,
         trial_evaluator,
         test_table,
@@ -431,71 +439,76 @@ impl Default for BoundaryAssemblerOptions {
     }
 }
 
-// TODO: make this a struct, with HasBoundaryAssemblerOptions trait
-pub trait BoundaryAssembler: Sync + Sized {
-    //! Boundary assembler
-    //!
-    //! Assemble operators by processing batches of cells in parallel
-
-    /// Scalar type
-    type T: RlstScalar + MatrixInverse;
-    /// Integrand type
-    type Integrand: BoundaryIntegrand<T = Self::T>;
-    /// Kernel type
-    type Kernel: KernelEvaluator<T = Self::T>;
-    /// Number of derivatives
-    const DERIV_SIZE: usize;
-    /// Number of derivatives needed in basis function tables
-    const TABLE_DERIVS: usize;
-
-    /// Get integrand
-    fn integrand(&self) -> &Self::Integrand;
-
-    /// Get integrand
-    fn kernel(&self) -> &Self::Kernel;
+/// Boundary assembler
+///
+/// Assembles operators by processing batches of cells in parallel
+pub struct BoundaryAssembler<
+    T: RlstScalar + MatrixInverse,
+    Integrand: BoundaryIntegrand<T = T>,
+    Kernel: KernelEvaluator<T = T>,
+> {
+    pub(crate) integrand: Integrand,
+    pub(crate) kernel: Kernel,
+    pub(crate) options: BoundaryAssemblerOptions,
+    pub(crate) deriv_size: usize,
+    pub(crate) table_derivs: usize,
+}
 
-    /// Get assembler options
-    fn options(&self) -> &BoundaryAssemblerOptions;
+unsafe impl<
+        T: RlstScalar + MatrixInverse,
+        Integrand: BoundaryIntegrand<T = T>,
+        Kernel: KernelEvaluator<T = T>,
+    > Sync for BoundaryAssembler<T, Integrand, Kernel>
+{
+}
 
-    /// Get mutable assembler options
-    fn options_mut(&mut self) -> &mut BoundaryAssemblerOptions;
+impl<
+        T: RlstScalar + MatrixInverse,
+        Integrand: BoundaryIntegrand<T = T>,
+        Kernel: KernelEvaluator<T = T>,
+    > BoundaryAssembler<T, Integrand, Kernel>
+{
+    /// Create new
+    fn new(integrand: Integrand, kernel: Kernel, deriv_size: usize, table_derivs: usize) -> Self {
+        Self {
+            integrand,
+            kernel,
+            options: BoundaryAssemblerOptions::default(),
+            deriv_size,
+            table_derivs,
+        }
+    }
 
     /// Set (non-singular) quadrature degree for a cell type
-    fn quadrature_degree(&mut self, cell: ReferenceCellType, degree: usize) {
-        *self
-            .options_mut()
-            .quadrature_degrees
-            .get_mut(&cell)
-            .unwrap() = degree;
+    pub fn quadrature_degree(&mut self, cell: ReferenceCellType, degree: usize) {
+        *self.options.quadrature_degrees.get_mut(&cell).unwrap() = degree;
     }
 
     /// Set singular quadrature degree for a pair of cell types
-    fn singular_quadrature_degree(
+    pub fn singular_quadrature_degree(
         &mut self,
         cells: (ReferenceCellType, ReferenceCellType),
         degree: usize,
     ) {
         *self
-            .options_mut()
+            .options
             .singular_quadrature_degrees
             .get_mut(&cells)
             .unwrap() = degree;
     }
 
     /// Set the maximum size of a batch of cells to send to an assembly function
-    fn batch_size(&mut self, size: usize) {
-        self.options_mut().batch_size = size;
+    pub fn batch_size(&mut self, size: usize) {
+        self.options.batch_size = size;
     }
-}
 
-trait InternalAssemblyFunctions: BoundaryAssembler {
     /// Assemble the singular contributions
-    fn assemble_singular<Space: FunctionSpace<T = Self::T> + Sync>(
+    fn assemble_singular<Space: FunctionSpace<T = T> + Sync>(
         &self,
         shape: [usize; 2],
         trial_space: &Space,
         test_space: &Space,
-    ) -> SparseMatrixData<Self::T> {
+    ) -> SparseMatrixData<T> {
         if !equal_grids(test_space.grid(), trial_space.grid()) {
             // If the test and trial grids are different, there are no neighbouring triangles
             return SparseMatrixData::new(shape);
@@ -522,8 +535,8 @@ trait InternalAssemblyFunctions: BoundaryAssembler {
 
         for test_cell_type in grid.entity_types(2) {
             for trial_cell_type in grid.entity_types(2) {
-                let qdegree = self.options().singular_quadrature_degrees
-                    [&(*test_cell_type, *trial_cell_type)];
+                let qdegree =
+                    self.options.singular_quadrature_degrees[&(*test_cell_type, *trial_cell_type)];
                 let offset = qweights.len();
 
                 let mut possible_pairs = vec![];
@@ -567,11 +580,11 @@ trait InternalAssemblyFunctions: BoundaryAssembler {
                     );
                     let npts = qrule.weights.len();
 
-                    let mut points = rlst_dynamic_array2!(<Self::T as RlstScalar>::Real, [2, npts]);
+                    let mut points = rlst_dynamic_array2!(<T as RlstScalar>::Real, [2, npts]);
                     for i in 0..npts {
                         for j in 0..2 {
                             *points.get_mut([j, i]).unwrap() =
-                                num::cast::<f64, <Self::T as RlstScalar>::Real>(
+                                num::cast::<f64, <T as RlstScalar>::Real>(
                                     qrule.trial_points[2 * i + j],
                                 )
                                 .unwrap();
@@ -579,18 +592,18 @@ trait InternalAssemblyFunctions: BoundaryAssembler {
                     }
                     let trial_element = trial_space.element(*trial_cell_type);
                     let mut table = rlst_dynamic_array4!(
-                        Self::T,
-                        trial_element.tabulate_array_shape(Self::TABLE_DERIVS, points.shape()[1])
+                        T,
+                        trial_element.tabulate_array_shape(self.table_derivs, points.shape()[1])
                     );
-                    trial_element.tabulate(&points, Self::TABLE_DERIVS, &mut table);
+                    trial_element.tabulate(&points, self.table_derivs, &mut table);
                     trial_points.push(points);
                     trial_tables.push(table);
 
-                    let mut points = rlst_dynamic_array2!(<Self::T as RlstScalar>::Real, [2, npts]);
+                    let mut points = rlst_dynamic_array2!(<T as RlstScalar>::Real, [2, npts]);
                     for i in 0..npts {
                         for j in 0..2 {
                             *points.get_mut([j, i]).unwrap() =
-                                num::cast::<f64, <Self::T as RlstScalar>::Real>(
+                                num::cast::<f64, <T as RlstScalar>::Real>(
                                     qrule.test_points[2 * i + j],
                                 )
                                 .unwrap();
@@ -598,17 +611,17 @@ trait InternalAssemblyFunctions: BoundaryAssembler {
                     }
                     let test_element = test_space.element(*test_cell_type);
                     let mut table = rlst_dynamic_array4!(
-                        Self::T,
-                        test_element.tabulate_array_shape(Self::TABLE_DERIVS, points.shape()[1])
+                        T,
+                        test_element.tabulate_array_shape(self.table_derivs, points.shape()[1])
                     );
-                    test_element.tabulate(&points, Self::TABLE_DERIVS, &mut table);
+                    test_element.tabulate(&points, self.table_derivs, &mut table);
                     test_points.push(points);
                     test_tables.push(table);
                     qweights.push(
                         qrule
                             .weights
                             .iter()
-                            .map(|w| num::cast::<f64, <Self::T as RlstScalar>::Real>(*w).unwrap())
+                            .map(|w| num::cast::<f64, <T as RlstScalar>::Real>(*w).unwrap())
                             .collect::<Vec<_>>(),
                     );
                 }
@@ -620,7 +633,7 @@ trait InternalAssemblyFunctions: BoundaryAssembler {
             },
             pair_indices.len(),
             grid,
-            self.options().batch_size,
+            self.options.batch_size,
         );
 
         cell_blocks
@@ -628,7 +641,7 @@ trait InternalAssemblyFunctions: BoundaryAssembler {
             .map(|(i, cell_block)| {
                 assemble_batch_singular(
                     self,
-                    Self::DERIV_SIZE,
+                    self.deriv_size,
                     shape,
                     trial_cell_types[i],
                     test_cell_types[i],
@@ -643,7 +656,7 @@ trait InternalAssemblyFunctions: BoundaryAssembler {
                 )
             })
             .reduce(
-                || SparseMatrixData::<Self::T>::new(shape),
+                || SparseMatrixData::<T>::new(shape),
                 |mut a, b| {
                     a.add(b);
                     a
@@ -654,12 +667,12 @@ trait InternalAssemblyFunctions: BoundaryAssembler {
     /// 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 = Self::T> + Sync>(
+    fn assemble_singular_correction<Space: FunctionSpace<T = T> + Sync>(
         &self,
         shape: [usize; 2],
         trial_space: &Space,
         test_space: &Space,
-    ) -> SparseMatrixData<Self::T> {
+    ) -> SparseMatrixData<T> {
         if !equal_grids(test_space.grid(), trial_space.grid()) {
             // If the test and trial grids are different, there are no neighbouring triangles
             return SparseMatrixData::new(shape);
@@ -686,48 +699,44 @@ trait InternalAssemblyFunctions: BoundaryAssembler {
         let mut cell_type_indices = HashMap::new();
 
         for test_cell_type in grid.entity_types(2) {
-            let npts_test = self.options().quadrature_degrees[test_cell_type];
+            let npts_test = self.options.quadrature_degrees[test_cell_type];
             for trial_cell_type in grid.entity_types(2) {
-                let npts_trial = self.options().quadrature_degrees[trial_cell_type];
+                let npts_trial = self.options.quadrature_degrees[trial_cell_type];
                 test_cell_types.push(*test_cell_type);
                 trial_cell_types.push(*trial_cell_type);
                 cell_type_indices.insert((*test_cell_type, *trial_cell_type), qweights_test.len());
 
                 let qrule_test = simplex_rule(*test_cell_type, npts_test).unwrap();
-                let mut test_pts =
-                    rlst_dynamic_array2!(<Self::T as RlstScalar>::Real, [2, npts_test]);
+                let mut test_pts = rlst_dynamic_array2!(<T as RlstScalar>::Real, [2, npts_test]);
                 for i in 0..npts_test {
                     for j in 0..2 {
                         *test_pts.get_mut([j, i]).unwrap() =
-                            num::cast::<f64, <Self::T as RlstScalar>::Real>(
-                                qrule_test.points[2 * i + j],
-                            )
-                            .unwrap();
+                            num::cast::<f64, <T as RlstScalar>::Real>(qrule_test.points[2 * i + j])
+                                .unwrap();
                     }
                 }
                 qweights_test.push(
                     qrule_test
                         .weights
                         .iter()
-                        .map(|w| num::cast::<f64, <Self::T as RlstScalar>::Real>(*w).unwrap())
+                        .map(|w| num::cast::<f64, <T as RlstScalar>::Real>(*w).unwrap())
                         .collect::<Vec<_>>(),
                 );
                 let test_element = test_space.element(*test_cell_type);
                 let mut test_table = rlst_dynamic_array4!(
-                    Self::T,
-                    test_element.tabulate_array_shape(Self::TABLE_DERIVS, npts_test)
+                    T,
+                    test_element.tabulate_array_shape(self.table_derivs, npts_test)
                 );
-                test_element.tabulate(&test_pts, Self::TABLE_DERIVS, &mut test_table);
+                test_element.tabulate(&test_pts, self.table_derivs, &mut test_table);
                 test_tables.push(test_table);
                 qpoints_test.push(test_pts);
 
                 let qrule_trial = simplex_rule(*trial_cell_type, npts_trial).unwrap();
-                let mut trial_pts =
-                    rlst_dynamic_array2!(<Self::T as RlstScalar>::Real, [2, npts_trial]);
+                let mut trial_pts = rlst_dynamic_array2!(<T as RlstScalar>::Real, [2, npts_trial]);
                 for i in 0..npts_trial {
                     for j in 0..2 {
                         *trial_pts.get_mut([j, i]).unwrap() =
-                            num::cast::<f64, <Self::T as RlstScalar>::Real>(
+                            num::cast::<f64, <T as RlstScalar>::Real>(
                                 qrule_trial.points[2 * i + j],
                             )
                             .unwrap();
@@ -737,15 +746,15 @@ trait InternalAssemblyFunctions: BoundaryAssembler {
                     qrule_trial
                         .weights
                         .iter()
-                        .map(|w| num::cast::<f64, <Self::T as RlstScalar>::Real>(*w).unwrap())
+                        .map(|w| num::cast::<f64, <T as RlstScalar>::Real>(*w).unwrap())
                         .collect::<Vec<_>>(),
                 );
                 let trial_element = trial_space.element(*trial_cell_type);
                 let mut trial_table = rlst_dynamic_array4!(
-                    Self::T,
-                    trial_element.tabulate_array_shape(Self::TABLE_DERIVS, npts_trial)
+                    T,
+                    trial_element.tabulate_array_shape(self.table_derivs, npts_trial)
                 );
-                trial_element.tabulate(&trial_pts, Self::TABLE_DERIVS, &mut trial_table);
+                trial_element.tabulate(&trial_pts, self.table_derivs, &mut trial_table);
                 trial_tables.push(trial_table);
                 qpoints_trial.push(trial_pts);
             }
@@ -757,7 +766,7 @@ trait InternalAssemblyFunctions: BoundaryAssembler {
             },
             qweights_test.len(),
             grid,
-            self.options().batch_size,
+            self.options.batch_size,
         );
 
         cell_blocks
@@ -765,7 +774,7 @@ trait InternalAssemblyFunctions: BoundaryAssembler {
             .map(|(i, cell_block)| {
                 assemble_batch_singular_correction(
                     self,
-                    Self::DERIV_SIZE,
+                    self.deriv_size,
                     shape,
                     trial_cell_types[i],
                     test_cell_types[i],
@@ -781,7 +790,7 @@ trait InternalAssemblyFunctions: BoundaryAssembler {
                 )
             })
             .reduce(
-                || SparseMatrixData::<Self::T>::new(shape),
+                || SparseMatrixData::<T>::new(shape),
                 |mut a, b| {
                     a.add(b);
                     a
@@ -789,9 +798,9 @@ trait InternalAssemblyFunctions: BoundaryAssembler {
             )
     }
     /// Assemble the non-singular contributions into a dense matrix
-    fn assemble_nonsingular<Space: FunctionSpace<T = Self::T> + Sync>(
+    fn assemble_nonsingular<Space: FunctionSpace<T = T> + Sync>(
         &self,
-        output: &RawData2D<Self::T>,
+        output: &RawData2D<T>,
         trial_space: &Space,
         test_space: &Space,
         trial_colouring: &HashMap<ReferenceCellType, Vec<Vec<usize>>>,
@@ -806,36 +815,34 @@ trait InternalAssemblyFunctions: BoundaryAssembler {
             panic!("Matrix has wrong shape");
         }
 
-        let batch_size = self.options().batch_size;
+        let batch_size = self.options.batch_size;
 
         for test_cell_type in test_space.grid().entity_types(2) {
-            let npts_test = self.options().quadrature_degrees[test_cell_type];
+            let npts_test = self.options.quadrature_degrees[test_cell_type];
             for trial_cell_type in trial_space.grid().entity_types(2) {
-                let npts_trial = self.options().quadrature_degrees[trial_cell_type];
+                let npts_trial = self.options.quadrature_degrees[trial_cell_type];
                 let qrule_test = simplex_rule(*test_cell_type, npts_test).unwrap();
                 let mut qpoints_test =
-                    rlst_dynamic_array2!(<Self::T as RlstScalar>::Real, [2, npts_test]);
+                    rlst_dynamic_array2!(<T as RlstScalar>::Real, [2, npts_test]);
                 for i in 0..npts_test {
                     for j in 0..2 {
                         *qpoints_test.get_mut([j, i]).unwrap() =
-                            num::cast::<f64, <Self::T as RlstScalar>::Real>(
-                                qrule_test.points[2 * i + j],
-                            )
-                            .unwrap();
+                            num::cast::<f64, <T as RlstScalar>::Real>(qrule_test.points[2 * i + j])
+                                .unwrap();
                     }
                 }
                 let qweights_test = qrule_test
                     .weights
                     .iter()
-                    .map(|w| num::cast::<f64, <Self::T as RlstScalar>::Real>(*w).unwrap())
+                    .map(|w| num::cast::<f64, <T as RlstScalar>::Real>(*w).unwrap())
                     .collect::<Vec<_>>();
                 let qrule_trial = simplex_rule(*trial_cell_type, npts_trial).unwrap();
                 let mut qpoints_trial =
-                    rlst_dynamic_array2!(<Self::T as RlstScalar>::Real, [2, npts_trial]);
+                    rlst_dynamic_array2!(<T as RlstScalar>::Real, [2, npts_trial]);
                 for i in 0..npts_trial {
                     for j in 0..2 {
                         *qpoints_trial.get_mut([j, i]).unwrap() =
-                            num::cast::<f64, <Self::T as RlstScalar>::Real>(
+                            num::cast::<f64, <T as RlstScalar>::Real>(
                                 qrule_trial.points[2 * i + j],
                             )
                             .unwrap();
@@ -844,22 +851,22 @@ trait InternalAssemblyFunctions: BoundaryAssembler {
                 let qweights_trial = qrule_trial
                     .weights
                     .iter()
-                    .map(|w| num::cast::<f64, <Self::T as RlstScalar>::Real>(*w).unwrap())
+                    .map(|w| num::cast::<f64, <T as RlstScalar>::Real>(*w).unwrap())
                     .collect::<Vec<_>>();
 
                 let test_element = test_space.element(*test_cell_type);
                 let mut test_table = rlst_dynamic_array4!(
-                    Self::T,
-                    test_element.tabulate_array_shape(Self::TABLE_DERIVS, npts_test)
+                    T,
+                    test_element.tabulate_array_shape(self.table_derivs, npts_test)
                 );
-                test_element.tabulate(&qpoints_test, Self::TABLE_DERIVS, &mut test_table);
+                test_element.tabulate(&qpoints_test, self.table_derivs, &mut test_table);
 
                 let trial_element = trial_space.element(*trial_cell_type);
                 let mut trial_table = rlst_dynamic_array4!(
-                    Self::T,
-                    trial_element.tabulate_array_shape(Self::TABLE_DERIVS, npts_trial)
+                    T,
+                    trial_element.tabulate_array_shape(self.table_derivs, npts_trial)
                 );
-                trial_element.tabulate(&qpoints_test, Self::TABLE_DERIVS, &mut trial_table);
+                trial_element.tabulate(&qpoints_test, self.table_derivs, &mut trial_table);
 
                 for test_c in &test_colouring[test_cell_type] {
                     for trial_c in &trial_colouring[trial_cell_type] {
@@ -894,7 +901,7 @@ trait InternalAssemblyFunctions: BoundaryAssembler {
                             .map(&|t| {
                                 assemble_batch_nonadjacent(
                                     self,
-                                    Self::DERIV_SIZE,
+                                    self.deriv_size,
                                     output,
                                     *test_cell_type,
                                     *trial_cell_type,
@@ -919,12 +926,16 @@ trait InternalAssemblyFunctions: BoundaryAssembler {
     }
 }
 
-impl<A: BoundaryAssembler> InternalAssemblyFunctions for A {}
-impl<A: BoundaryAssembler + InternalAssemblyFunctions> BoundaryAssembly for A {
-    type T = A::T;
-    fn assemble_singular_into_dense<Space: FunctionSpace<T = A::T> + Sync>(
+impl<
+        T: RlstScalar + MatrixInverse,
+        Integrand: BoundaryIntegrand<T = T>,
+        Kernel: KernelEvaluator<T = T>,
+    > BoundaryAssembly for BoundaryAssembler<T, Integrand, Kernel>
+{
+    type T = T;
+    fn assemble_singular_into_dense<Space: FunctionSpace<T = T> + Sync>(
         &self,
-        output: &mut RlstArray<Self::T, 2>,
+        output: &mut RlstArray<T, 2>,
         trial_space: &Space,
         test_space: &Space,
     ) {
@@ -937,15 +948,15 @@ impl<A: BoundaryAssembler + InternalAssemblyFunctions> BoundaryAssembly for A {
         }
     }
 
-    fn assemble_singular_into_csr<Space: FunctionSpace<T = A::T> + Sync>(
+    fn assemble_singular_into_csr<Space: FunctionSpace<T = T> + Sync>(
         &self,
         trial_space: &Space,
         test_space: &Space,
-    ) -> CsrMatrix<Self::T> {
+    ) -> CsrMatrix<T> {
         let shape = [test_space.global_size(), trial_space.global_size()];
         let sparse_matrix = self.assemble_singular(shape, trial_space, test_space);
 
-        CsrMatrix::<Self::T>::from_aij(
+        CsrMatrix::<T>::from_aij(
             sparse_matrix.shape,
             &sparse_matrix.rows,
             &sparse_matrix.cols,
@@ -954,9 +965,9 @@ impl<A: BoundaryAssembler + InternalAssemblyFunctions> BoundaryAssembly for A {
         .unwrap()
     }
 
-    fn assemble_singular_correction_into_dense<Space: FunctionSpace<T = A::T> + Sync>(
+    fn assemble_singular_correction_into_dense<Space: FunctionSpace<T = T> + Sync>(
         &self,
-        output: &mut RlstArray<Self::T, 2>,
+        output: &mut RlstArray<T, 2>,
         trial_space: &Space,
         test_space: &Space,
     ) {
@@ -970,15 +981,15 @@ impl<A: BoundaryAssembler + InternalAssemblyFunctions> BoundaryAssembly for A {
         }
     }
 
-    fn assemble_singular_correction_into_csr<Space: FunctionSpace<T = A::T> + Sync>(
+    fn assemble_singular_correction_into_csr<Space: FunctionSpace<T = T> + Sync>(
         &self,
         trial_space: &Space,
         test_space: &Space,
-    ) -> CsrMatrix<Self::T> {
+    ) -> CsrMatrix<T> {
         let shape = [test_space.global_size(), trial_space.global_size()];
         let sparse_matrix = self.assemble_singular_correction(shape, trial_space, test_space);
 
-        CsrMatrix::<Self::T>::from_aij(
+        CsrMatrix::<T>::from_aij(
             sparse_matrix.shape,
             &sparse_matrix.rows,
             &sparse_matrix.cols,
@@ -987,9 +998,9 @@ impl<A: BoundaryAssembler + InternalAssemblyFunctions> BoundaryAssembly for A {
         .unwrap()
     }
 
-    fn assemble_into_dense<Space: FunctionSpace<T = A::T> + Sync>(
+    fn assemble_into_dense<Space: FunctionSpace<T = T> + Sync>(
         &self,
-        output: &mut RlstArray<Self::T, 2>,
+        output: &mut RlstArray<T, 2>,
         trial_space: &Space,
         test_space: &Space,
     ) {
@@ -1006,9 +1017,9 @@ impl<A: BoundaryAssembler + InternalAssemblyFunctions> BoundaryAssembly for A {
         self.assemble_singular_into_dense(output, trial_space, test_space);
     }
 
-    fn assemble_nonsingular_into_dense<Space: FunctionSpace<T = A::T> + Sync>(
+    fn assemble_nonsingular_into_dense<Space: FunctionSpace<T = T> + Sync>(
         &self,
-        output: &mut RlstArray<Self::T, 2>,
+        output: &mut RlstArray<T, 2>,
         trial_space: &Space,
         test_space: &Space,
         trial_colouring: &HashMap<ReferenceCellType, Vec<Vec<usize>>>,
@@ -1038,26 +1049,31 @@ impl<A: BoundaryAssembler + InternalAssemblyFunctions> BoundaryAssembly for A {
     }
 }
 #[cfg(feature = "mpi")]
-impl<A: BoundaryAssembler + InternalAssemblyFunctions> ParallelBoundaryAssembly for A {
+impl<
+        T: RlstScalar + MatrixInverse,
+        Integrand: BoundaryIntegrand<T = T>,
+        Kernel: KernelEvaluator<T = T>,
+    > ParallelBoundaryAssembly for BoundaryAssembler<T, Integrand, Kernel>
+{
     fn parallel_assemble_singular_into_csr<
         C: Communicator,
-        Space: ParallelFunctionSpace<C, T = A::T>,
+        Space: ParallelFunctionSpace<C, T = T>,
     >(
         &self,
         trial_space: &Space,
         test_space: &Space,
-    ) -> CsrMatrix<Self::T> {
+    ) -> CsrMatrix<T> {
         self.assemble_singular_into_csr(trial_space.local_space(), test_space.local_space())
     }
 
     fn parallel_assemble_singular_correction_into_csr<
         C: Communicator,
-        Space: ParallelFunctionSpace<C, T = A::T>,
+        Space: ParallelFunctionSpace<C, T = T>,
     >(
         &self,
         trial_space: &Space,
         test_space: &Space,
-    ) -> CsrMatrix<Self::T> {
+    ) -> CsrMatrix<T> {
         self.assemble_singular_correction_into_csr(
             trial_space.local_space(),
             test_space.local_space(),
diff --git a/src/assembly/boundary/assemblers/adjoint_double_layer.rs b/src/assembly/boundary/assemblers/adjoint_double_layer.rs
index ee6ed46b..f1e01018 100644
--- a/src/assembly/boundary/assemblers/adjoint_double_layer.rs
+++ b/src/assembly/boundary/assemblers/adjoint_double_layer.rs
@@ -1,5 +1,5 @@
 //! Adjoint double layer assemblers
-use super::{BoundaryAssembler, BoundaryAssemblerOptions};
+use super::BoundaryAssembler;
 use crate::assembly::{
     boundary::integrands::AdjointDoubleLayerBoundaryIntegrand, common::GreenKernelEvalType,
     kernels::KernelEvaluator,
@@ -7,60 +7,40 @@ use crate::assembly::{
 use green_kernels::{helmholtz_3d::Helmholtz3dKernel, laplace_3d::Laplace3dKernel, traits::Kernel};
 use rlst::{MatrixInverse, RlstScalar};
 
-/// Assembler for a adjoint double layer operator
-pub struct AdjointDoubleLayerAssembler<T: RlstScalar + MatrixInverse, K: Kernel<T = T>> {
-    integrand: AdjointDoubleLayerBoundaryIntegrand<T>,
-    kernel: KernelEvaluator<T, K>,
-    options: BoundaryAssemblerOptions,
-}
-impl<T: RlstScalar + MatrixInverse, K: Kernel<T = T>> AdjointDoubleLayerAssembler<T, K> {
+impl<T: RlstScalar + MatrixInverse, K: Kernel<T = T>>
+    BoundaryAssembler<T, AdjointDoubleLayerBoundaryIntegrand<T>, KernelEvaluator<T, K>>
+{
     /// Create a new adjoint double layer assembler
-    pub fn new(kernel: KernelEvaluator<T, K>) -> Self {
-        Self {
-            integrand: AdjointDoubleLayerBoundaryIntegrand::new(),
-            kernel,
-            options: BoundaryAssemblerOptions::default(),
-        }
+    pub fn new_adjoint_double_layer(kernel: KernelEvaluator<T, K>) -> Self {
+        Self::new(AdjointDoubleLayerBoundaryIntegrand::new(), kernel, 4, 0)
     }
 }
-impl<T: RlstScalar + MatrixInverse> AdjointDoubleLayerAssembler<T, Laplace3dKernel<T>> {
+impl<T: RlstScalar + MatrixInverse>
+    BoundaryAssembler<
+        T,
+        AdjointDoubleLayerBoundaryIntegrand<T>,
+        KernelEvaluator<T, Laplace3dKernel<T>>,
+    >
+{
     /// Create a new Laplace adjoint double layer assembler
-    pub fn new_laplace() -> Self {
-        Self::new(KernelEvaluator::new_laplace(
+    pub fn new_laplace_adjoint_double_layer() -> Self {
+        Self::new_adjoint_double_layer(KernelEvaluator::new_laplace(
             GreenKernelEvalType::ValueDeriv,
         ))
     }
 }
 impl<T: RlstScalar<Complex = T> + MatrixInverse>
-    AdjointDoubleLayerAssembler<T, Helmholtz3dKernel<T>>
+    BoundaryAssembler<
+        T,
+        AdjointDoubleLayerBoundaryIntegrand<T>,
+        KernelEvaluator<T, Helmholtz3dKernel<T>>,
+    >
 {
     /// Create a new Helmholtz adjoint double layer assembler
-    pub fn new_helmholtz(wavenumber: T::Real) -> Self {
-        Self::new(KernelEvaluator::new_helmholtz(
+    pub fn new_helmholtz_adjoint_double_layer(wavenumber: T::Real) -> Self {
+        Self::new_adjoint_double_layer(KernelEvaluator::new_helmholtz(
             wavenumber,
             GreenKernelEvalType::ValueDeriv,
         ))
     }
 }
-
-impl<T: RlstScalar + MatrixInverse, K: Kernel<T = T>> BoundaryAssembler
-    for AdjointDoubleLayerAssembler<T, K>
-{
-    const DERIV_SIZE: usize = 4;
-    const TABLE_DERIVS: usize = 0;
-    type T = T;
-    type Integrand = AdjointDoubleLayerBoundaryIntegrand<T>;
-    type Kernel = KernelEvaluator<T, K>;
-    fn integrand(&self) -> &AdjointDoubleLayerBoundaryIntegrand<T> {
-        &self.integrand
-    }
-    fn kernel(&self) -> &KernelEvaluator<T, K> {
-        &self.kernel
-    }
-    fn options(&self) -> &BoundaryAssemblerOptions {
-        &self.options
-    }
-    fn options_mut(&mut self) -> &mut BoundaryAssemblerOptions {
-        &mut self.options
-    }
-}
diff --git a/src/assembly/boundary/assemblers/double_layer.rs b/src/assembly/boundary/assemblers/double_layer.rs
index 8e441cd2..f837ea5e 100644
--- a/src/assembly/boundary/assemblers/double_layer.rs
+++ b/src/assembly/boundary/assemblers/double_layer.rs
@@ -1,5 +1,5 @@
 //! Double layer assemblers
-use super::{BoundaryAssembler, BoundaryAssemblerOptions};
+use super::BoundaryAssembler;
 use crate::assembly::{
     boundary::integrands::DoubleLayerBoundaryIntegrand, common::GreenKernelEvalType,
     kernels::KernelEvaluator,
@@ -7,58 +7,32 @@ use crate::assembly::{
 use green_kernels::{helmholtz_3d::Helmholtz3dKernel, laplace_3d::Laplace3dKernel, traits::Kernel};
 use rlst::{MatrixInverse, RlstScalar};
 
-/// Assembler for a double layer operator
-pub struct DoubleLayerAssembler<T: RlstScalar + MatrixInverse, K: Kernel<T = T>> {
-    integrand: DoubleLayerBoundaryIntegrand<T>,
-    kernel: KernelEvaluator<T, K>,
-    options: BoundaryAssemblerOptions,
-}
-impl<T: RlstScalar + MatrixInverse, K: Kernel<T = T>> DoubleLayerAssembler<T, K> {
+impl<T: RlstScalar + MatrixInverse, K: Kernel<T = T>>
+    BoundaryAssembler<T, DoubleLayerBoundaryIntegrand<T>, KernelEvaluator<T, K>>
+{
     /// Create a new double layer assembler
-    pub fn new(kernel: KernelEvaluator<T, K>) -> Self {
-        Self {
-            integrand: DoubleLayerBoundaryIntegrand::new(),
-            kernel,
-            options: BoundaryAssemblerOptions::default(),
-        }
+    pub fn new_double_layer(kernel: KernelEvaluator<T, K>) -> Self {
+        Self::new(DoubleLayerBoundaryIntegrand::new(), kernel, 4, 0)
     }
 }
-impl<T: RlstScalar + MatrixInverse> DoubleLayerAssembler<T, Laplace3dKernel<T>> {
+impl<T: RlstScalar + MatrixInverse>
+    BoundaryAssembler<T, DoubleLayerBoundaryIntegrand<T>, KernelEvaluator<T, Laplace3dKernel<T>>>
+{
     /// Create a new Laplace double layer assembler
-    pub fn new_laplace() -> Self {
-        Self::new(KernelEvaluator::new_laplace(
+    pub fn new_laplace_double_layer() -> Self {
+        Self::new_double_layer(KernelEvaluator::new_laplace(
             GreenKernelEvalType::ValueDeriv,
         ))
     }
 }
-impl<T: RlstScalar<Complex = T> + MatrixInverse> DoubleLayerAssembler<T, Helmholtz3dKernel<T>> {
+impl<T: RlstScalar<Complex = T> + MatrixInverse>
+    BoundaryAssembler<T, DoubleLayerBoundaryIntegrand<T>, KernelEvaluator<T, Helmholtz3dKernel<T>>>
+{
     /// Create a new Helmholtz double layer assembler
-    pub fn new_helmholtz(wavenumber: T::Real) -> Self {
-        Self::new(KernelEvaluator::new_helmholtz(
+    pub fn new_helmholtz_double_layer(wavenumber: T::Real) -> Self {
+        Self::new_double_layer(KernelEvaluator::new_helmholtz(
             wavenumber,
             GreenKernelEvalType::ValueDeriv,
         ))
     }
 }
-
-impl<T: RlstScalar + MatrixInverse, K: Kernel<T = T>> BoundaryAssembler
-    for DoubleLayerAssembler<T, K>
-{
-    const DERIV_SIZE: usize = 4;
-    const TABLE_DERIVS: usize = 0;
-    type T = T;
-    type Integrand = DoubleLayerBoundaryIntegrand<T>;
-    type Kernel = KernelEvaluator<T, K>;
-    fn integrand(&self) -> &DoubleLayerBoundaryIntegrand<T> {
-        &self.integrand
-    }
-    fn kernel(&self) -> &KernelEvaluator<T, K> {
-        &self.kernel
-    }
-    fn options(&self) -> &BoundaryAssemblerOptions {
-        &self.options
-    }
-    fn options_mut(&mut self) -> &mut BoundaryAssemblerOptions {
-        &mut self.options
-    }
-}
diff --git a/src/assembly/boundary/assemblers/hypersingular.rs b/src/assembly/boundary/assemblers/hypersingular.rs
index e0fb5a1a..6a6b0ea8 100644
--- a/src/assembly/boundary/assemblers/hypersingular.rs
+++ b/src/assembly/boundary/assemblers/hypersingular.rs
@@ -1,5 +1,5 @@
 //! Hypersingular assemblers
-use super::{BoundaryAssembler, BoundaryAssemblerOptions};
+use super::BoundaryAssembler;
 use crate::assembly::{
     boundary::integrands::{
         BoundaryIntegrandSum, HypersingularCurlCurlBoundaryIntegrand,
@@ -18,73 +18,40 @@ type HelmholtzIntegrand<T> = BoundaryIntegrandSum<
     HypersingularNormalNormalBoundaryIntegrand<T>,
 >;
 
-/// Assembler for a hypersingular operator
-pub struct HypersingularAssembler<
-    T: RlstScalar + MatrixInverse,
-    K: Kernel<T = T>,
-    I: BoundaryIntegrand<T = T> + Sync,
-> {
-    kernel: KernelEvaluator<T, K>,
-    integrand: I,
-    options: BoundaryAssemblerOptions,
-}
-
-impl<T: RlstScalar + MatrixInverse, K: Kernel<T = T>, I: BoundaryIntegrand<T = T> + Sync>
-    HypersingularAssembler<T, K, I>
+impl<T: RlstScalar + MatrixInverse, K: Kernel<T = T>, I: BoundaryIntegrand<T = T>>
+    BoundaryAssembler<T, I, KernelEvaluator<T, K>>
 {
-    /// Create a new hypersingular assembler
-    pub fn new(kernel: KernelEvaluator<T, K>, integrand: I) -> Self {
-        Self {
-            kernel,
-            integrand,
-            options: BoundaryAssemblerOptions::default(),
-        }
+    /// Create a new adjoint double layer assembler
+    pub fn new_hypersingular(integrand: I, kernel: KernelEvaluator<T, K>) -> Self {
+        Self::new(integrand, kernel, 4, 1)
     }
 }
 impl<T: RlstScalar + MatrixInverse>
-    HypersingularAssembler<T, Laplace3dKernel<T>, HypersingularCurlCurlBoundaryIntegrand<T>>
+    BoundaryAssembler<
+        T,
+        HypersingularCurlCurlBoundaryIntegrand<T>,
+        KernelEvaluator<T, Laplace3dKernel<T>>,
+    >
 {
-    /// Create a new Laplace hypersingular assembler
-    pub fn new_laplace() -> Self {
-        Self::new(
-            KernelEvaluator::new_laplace(GreenKernelEvalType::ValueDeriv),
+    /// Create a new Laplace adjoint double layer assembler
+    pub fn new_laplace_hypersingular() -> Self {
+        Self::new_hypersingular(
             HypersingularCurlCurlBoundaryIntegrand::new(),
+            KernelEvaluator::new_laplace(GreenKernelEvalType::ValueDeriv),
         )
     }
 }
 impl<T: RlstScalar<Complex = T> + MatrixInverse>
-    HypersingularAssembler<T, Helmholtz3dKernel<T>, HelmholtzIntegrand<T>>
+    BoundaryAssembler<T, HelmholtzIntegrand<T>, KernelEvaluator<T, Helmholtz3dKernel<T>>>
 {
-    /// Create a new Helmholtz hypersingular assembler
-    pub fn new_helmholtz(wavenumber: T::Real) -> Self {
-        Self::new(
-            KernelEvaluator::new_helmholtz(wavenumber, GreenKernelEvalType::ValueDeriv),
+    /// Create a new Helmholtz adjoint double layer assembler
+    pub fn new_helmholtz_hypersingular(wavenumber: T::Real) -> Self {
+        Self::new_hypersingular(
             BoundaryIntegrandSum::new(
                 HypersingularCurlCurlBoundaryIntegrand::new(),
                 HypersingularNormalNormalBoundaryIntegrand::new(wavenumber),
             ),
+            KernelEvaluator::new_helmholtz(wavenumber, GreenKernelEvalType::ValueDeriv),
         )
     }
 }
-
-impl<T: RlstScalar + MatrixInverse, K: Kernel<T = T>, I: BoundaryIntegrand<T = T> + Sync>
-    BoundaryAssembler for HypersingularAssembler<T, K, I>
-{
-    const DERIV_SIZE: usize = 4;
-    const TABLE_DERIVS: usize = 1;
-    type T = T;
-    type Integrand = I;
-    type Kernel = KernelEvaluator<T, K>;
-    fn integrand(&self) -> &I {
-        &self.integrand
-    }
-    fn kernel(&self) -> &KernelEvaluator<T, K> {
-        &self.kernel
-    }
-    fn options(&self) -> &BoundaryAssemblerOptions {
-        &self.options
-    }
-    fn options_mut(&mut self) -> &mut BoundaryAssemblerOptions {
-        &mut self.options
-    }
-}
diff --git a/src/assembly/boundary/assemblers/single_layer.rs b/src/assembly/boundary/assemblers/single_layer.rs
index 4bf9f81a..90045f6d 100644
--- a/src/assembly/boundary/assemblers/single_layer.rs
+++ b/src/assembly/boundary/assemblers/single_layer.rs
@@ -1,5 +1,5 @@
 //! Single layer assemblers
-use super::{BoundaryAssembler, BoundaryAssemblerOptions};
+use super::BoundaryAssembler;
 use crate::assembly::{
     boundary::integrands::SingleLayerBoundaryIntegrand, common::GreenKernelEvalType,
     kernels::KernelEvaluator,
@@ -7,56 +7,30 @@ use crate::assembly::{
 use green_kernels::{helmholtz_3d::Helmholtz3dKernel, laplace_3d::Laplace3dKernel, traits::Kernel};
 use rlst::{MatrixInverse, RlstScalar};
 
-/// Assembler for a single layer operator
-pub struct SingleLayerAssembler<T: RlstScalar + MatrixInverse, K: Kernel<T = T>> {
-    integrand: SingleLayerBoundaryIntegrand<T>,
-    kernel: KernelEvaluator<T, K>,
-    options: BoundaryAssemblerOptions,
-}
-impl<T: RlstScalar + MatrixInverse, K: Kernel<T = T>> SingleLayerAssembler<T, K> {
+impl<T: RlstScalar + MatrixInverse, K: Kernel<T = T>>
+    BoundaryAssembler<T, SingleLayerBoundaryIntegrand<T>, KernelEvaluator<T, K>>
+{
     /// Create a new single layer assembler
-    pub fn new(kernel: KernelEvaluator<T, K>) -> Self {
-        Self {
-            integrand: SingleLayerBoundaryIntegrand::new(),
-            kernel,
-            options: BoundaryAssemblerOptions::default(),
-        }
+    pub fn new_single_layer(kernel: KernelEvaluator<T, K>) -> Self {
+        Self::new(SingleLayerBoundaryIntegrand::new(), kernel, 1, 0)
     }
 }
-impl<T: RlstScalar + MatrixInverse> SingleLayerAssembler<T, Laplace3dKernel<T>> {
+impl<T: RlstScalar + MatrixInverse>
+    BoundaryAssembler<T, SingleLayerBoundaryIntegrand<T>, KernelEvaluator<T, Laplace3dKernel<T>>>
+{
     /// Create a new Laplace single layer assembler
-    pub fn new_laplace() -> Self {
-        Self::new(KernelEvaluator::new_laplace(GreenKernelEvalType::Value))
+    pub fn new_laplace_single_layer() -> Self {
+        Self::new_single_layer(KernelEvaluator::new_laplace(GreenKernelEvalType::Value))
     }
 }
-impl<T: RlstScalar<Complex = T> + MatrixInverse> SingleLayerAssembler<T, Helmholtz3dKernel<T>> {
+impl<T: RlstScalar<Complex = T> + MatrixInverse>
+    BoundaryAssembler<T, SingleLayerBoundaryIntegrand<T>, KernelEvaluator<T, Helmholtz3dKernel<T>>>
+{
     /// Create a new Helmholtz single layer assembler
-    pub fn new_helmholtz(wavenumber: T::Real) -> Self {
-        Self::new(KernelEvaluator::new_helmholtz(
+    pub fn new_helmholtz_single_layer(wavenumber: T::Real) -> Self {
+        Self::new_single_layer(KernelEvaluator::new_helmholtz(
             wavenumber,
             GreenKernelEvalType::Value,
         ))
     }
 }
-
-impl<T: RlstScalar + MatrixInverse, K: Kernel<T = T>> BoundaryAssembler
-    for SingleLayerAssembler<T, K>
-{
-    const DERIV_SIZE: usize = 1;
-    const TABLE_DERIVS: usize = 0;
-    type T = T;
-    type Integrand = SingleLayerBoundaryIntegrand<T>;
-    type Kernel = KernelEvaluator<T, K>;
-    fn integrand(&self) -> &SingleLayerBoundaryIntegrand<T> {
-        &self.integrand
-    }
-    fn kernel(&self) -> &KernelEvaluator<T, K> {
-        &self.kernel
-    }
-    fn options(&self) -> &BoundaryAssemblerOptions {
-        &self.options
-    }
-    fn options_mut(&mut self) -> &mut BoundaryAssemblerOptions {
-        &mut self.options
-    }
-}


Test patch:
diff --git a/examples/test_parallel_assembly.rs b/examples/test_parallel_assembly.rs
index f92357b2..dad3ed94 100644
--- a/examples/test_parallel_assembly.rs
+++ b/examples/test_parallel_assembly.rs
@@ -4,7 +4,7 @@
 use approx::assert_relative_eq;
 #[cfg(feature = "mpi")]
 use bempp::{
-    assembly::boundary,
+    assembly::boundary::BoundaryAssembler,
     function::{ParallelFunctionSpace, SerialFunctionSpace},
     traits::{BoundaryAssembly, FunctionSpace, ParallelBoundaryAssembly},
 };
@@ -91,7 +91,7 @@ fn test_parallel_assembly_single_element_grid<C: Communicator>(
     let element = LagrangeElementFamily::<f64>::new(degree, cont);
     let space = ParallelFunctionSpace::new(&grid, &element);
 
-    let a = boundary::SingleLayerAssembler::<f64, _>::new_laplace();
+    let a = BoundaryAssembler::<f64, _, _>::new_laplace_single_layer();
 
     let matrix = a.parallel_assemble_singular_into_csr(&space, &space);
 
diff --git a/tests/compare_to_bempp_cl.rs b/tests/compare_to_bempp_cl.rs
index 83f6b791..d5e3d032 100644
--- a/tests/compare_to_bempp_cl.rs
+++ b/tests/compare_to_bempp_cl.rs
@@ -1,5 +1,5 @@
 use approx::*;
-use bempp::assembly::{boundary, potential, potential::PotentialAssembler};
+use bempp::assembly::{boundary::BoundaryAssembler, potential, potential::PotentialAssembler};
 use bempp::function::SerialFunctionSpace;
 use bempp::traits::{BoundaryAssembly, FunctionSpace};
 use cauchy::c64;
@@ -18,7 +18,7 @@ fn test_laplace_single_layer_dp0_dp0() {
 
     let mut matrix = rlst_dynamic_array2!(f64, [ndofs, ndofs]);
 
-    let a = boundary::SingleLayerAssembler::<f64, _>::new_laplace();
+    let a = BoundaryAssembler::<f64, _, _>::new_laplace_single_layer();
     a.assemble_into_dense(&mut matrix, &space, &space);
 
     // Compare to result from bempp-cl
@@ -41,7 +41,7 @@ fn test_laplace_double_layer_dp0_dp0() {
     let ndofs = space.global_size();
 
     let mut matrix = rlst_dynamic_array2!(f64, [ndofs, ndofs]);
-    let a = boundary::DoubleLayerAssembler::<f64, _>::new_laplace();
+    let a = BoundaryAssembler::<f64, _, _>::new_laplace_double_layer();
     a.assemble_into_dense(&mut matrix, &space, &space);
 
     // Compare to result from bempp-cl
@@ -63,7 +63,7 @@ fn test_laplace_adjoint_double_layer_dp0_dp0() {
     let ndofs = space.global_size();
 
     let mut matrix = rlst_dynamic_array2!(f64, [ndofs, ndofs]);
-    let a = boundary::AdjointDoubleLayerAssembler::<f64, _>::new_laplace();
+    let a = BoundaryAssembler::<f64, _, _>::new_laplace_adjoint_double_layer();
     a.assemble_into_dense(&mut matrix, &space, &space);
 
     // Compare to result from bempp-cl
@@ -86,7 +86,7 @@ fn test_laplace_hypersingular_dp0_dp0() {
     let ndofs = space.global_size();
 
     let mut matrix = rlst_dynamic_array2!(f64, [ndofs, ndofs]);
-    let a = boundary::HypersingularAssembler::<f64, _, _>::new_laplace();
+    let a = BoundaryAssembler::<f64, _, _>::new_laplace_hypersingular();
     a.assemble_into_dense(&mut matrix, &space, &space);
 
     for i in 0..ndofs {
@@ -105,7 +105,7 @@ fn test_laplace_hypersingular_p1_p1() {
     let ndofs = space.global_size();
 
     let mut matrix = rlst_dynamic_array2!(f64, [ndofs, ndofs]);
-    let a = boundary::HypersingularAssembler::<f64, _, _>::new_laplace();
+    let a = BoundaryAssembler::<f64, _, _>::new_laplace_hypersingular();
     a.assemble_into_dense(&mut matrix, &space, &space);
 
     // Compare to result from bempp-cl
@@ -134,7 +134,7 @@ fn test_helmholtz_single_layer_dp0_dp0() {
     let ndofs = space.global_size();
     let mut matrix = rlst_dynamic_array2!(c64, [ndofs, ndofs]);
 
-    let a = boundary::SingleLayerAssembler::<c64, _>::new_helmholtz(3.0);
+    let a = BoundaryAssembler::<c64, _, _>::new_helmholtz_single_layer(3.0);
     a.assemble_into_dense(&mut matrix, &space, &space);
 
     // Compare to result from bempp-cl
@@ -157,7 +157,7 @@ fn test_helmholtz_double_layer_dp0_dp0() {
     let ndofs = space.global_size();
     let mut matrix = rlst_dynamic_array2!(c64, [ndofs, ndofs]);
 
-    let a = boundary::DoubleLayerAssembler::<c64, _>::new_helmholtz(3.0);
+    let a = BoundaryAssembler::<c64, _, _>::new_helmholtz_double_layer(3.0);
     a.assemble_into_dense(&mut matrix, &space, &space);
 
     // Compare to result from bempp-cl
@@ -179,7 +179,7 @@ fn test_helmholtz_adjoint_double_layer_dp0_dp0() {
     let ndofs = space.global_size();
     let mut matrix = rlst_dynamic_array2!(c64, [ndofs, ndofs]);
 
-    let a = boundary::AdjointDoubleLayerAssembler::<c64, _>::new_helmholtz(3.0);
+    let a = BoundaryAssembler::<c64, _, _>::new_helmholtz_adjoint_double_layer(3.0);
     a.assemble_into_dense(&mut matrix, &space, &space);
 
     // Compare to result from bempp-cl
@@ -202,7 +202,7 @@ fn test_helmholtz_hypersingular_p1_p1() {
     let ndofs = space.global_size();
     let mut matrix = rlst_dynamic_array2!(c64, [ndofs, ndofs]);
 
-    let a = boundary::HypersingularAssembler::<c64, _, _>::new_helmholtz(3.0);
+    let a = BoundaryAssembler::<c64, _, _>::new_helmholtz_hypersingular(3.0);
     a.assemble_into_dense(&mut matrix, &space, &space);
 
     // Compare to result from bempp-cl