wrap potential assembly

bempp · Sep 19, 2024 · 11c14c1 · 11c14c1
1 parent 865c43c
commit 11c14c1
Show file tree

Hide file tree

Showing 7 changed files with 497 additions and 119 deletions.
diff --git a/python/bempp/assembly/potential.py b/python/bempp/assembly/potential.py
@@ -8,7 +8,6 @@
 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,
@@ -48,6 +47,23 @@ def __del__(self):
         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):

diff --git 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
@@ -30,3 +30,52 @@ def test_create_assembler(otype):
     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)
diff --git 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],
@@ -154,11 +155,14 @@ impl<
         self.options.batch_size
     }
 
-    fn assemble<Space: FunctionSpace<T = T> + Sync>(
+    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() {
@@ -237,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");