Move wavenumber out of HypersingularNormalNormalBoundaryIntegrand (#289)

* BoundaryIntegrandScalarProduct * remove extra minus sign * clippy
bempp · Aug 28, 2024 · 5f88369 · 5f88369
1 parent f2017b7
commit 5f88369
Show file tree

Hide file tree

Showing 3 changed files with 101 additions and 28 deletions.
diff --git a/src/assembly/boundary/assemblers/hypersingular.rs b/src/assembly/boundary/assemblers/hypersingular.rs
@@ -2,8 +2,8 @@
 use super::BoundaryAssembler;
 use crate::assembly::{
     boundary::integrands::{
-        BoundaryIntegrandSum, HypersingularCurlCurlBoundaryIntegrand,
-        HypersingularNormalNormalBoundaryIntegrand,
+        BoundaryIntegrandScalarProduct, BoundaryIntegrandSum,
+        HypersingularCurlCurlBoundaryIntegrand, HypersingularNormalNormalBoundaryIntegrand,
     },
     common::GreenKernelEvalType,
     kernels::KernelEvaluator,
@@ -15,7 +15,7 @@ use rlst::{MatrixInverse, RlstScalar};
 type HelmholtzIntegrand<T> = BoundaryIntegrandSum<
     T,
     HypersingularCurlCurlBoundaryIntegrand<T>,
-    HypersingularNormalNormalBoundaryIntegrand<T>,
+    BoundaryIntegrandScalarProduct<T, HypersingularNormalNormalBoundaryIntegrand<T>>,
 >;
 
 impl<T: RlstScalar + MatrixInverse, K: Kernel<T = T>, I: BoundaryIntegrand<T = T>>
@@ -49,7 +49,10 @@ impl<T: RlstScalar<Complex = T> + MatrixInverse>
         Self::new_hypersingular(
             BoundaryIntegrandSum::new(
                 HypersingularCurlCurlBoundaryIntegrand::new(),
-                HypersingularNormalNormalBoundaryIntegrand::new(wavenumber),
+                BoundaryIntegrandScalarProduct::new(
+                    num::cast::<T::Real, T>(-wavenumber.powi(2)).unwrap(),
+                    HypersingularNormalNormalBoundaryIntegrand::new(),
+                ),
             ),
             KernelEvaluator::new_helmholtz(wavenumber, GreenKernelEvalType::ValueDeriv),
         )

diff --git a/src/assembly/boundary/integrands.rs b/src/assembly/boundary/integrands.rs
@@ -108,3 +108,71 @@ impl<T: RlstScalar, I0: BoundaryIntegrand<T = T>, I1: BoundaryIntegrand<T = T>>
         )
     }
 }
+
+/// An integrand multiplied by a scalar
+pub struct BoundaryIntegrandScalarProduct<T: RlstScalar, I: BoundaryIntegrand<T = T>> {
+    scalar: T,
+    integrand: I,
+}
+
+impl<T: RlstScalar, I: BoundaryIntegrand<T = T>> BoundaryIntegrandScalarProduct<T, I> {
+    pub fn new(scalar: T, integrand: I) -> Self {
+        Self { scalar, integrand }
+    }
+}
+
+impl<T: RlstScalar, I: BoundaryIntegrand<T = T>> BoundaryIntegrand
+    for BoundaryIntegrandScalarProduct<T, I>
+{
+    type T = T;
+
+    unsafe fn evaluate_nonsingular(
+        &self,
+        test_table: &RlstArray<T, 4>,
+        trial_table: &RlstArray<T, 4>,
+        test_point_index: usize,
+        trial_point_index: usize,
+        test_basis_index: usize,
+        trial_basis_index: usize,
+        k: &RlstArray<T, 3>,
+        test_geometry: &impl CellGeometry<T = T::Real>,
+        trial_geometry: &impl CellGeometry<T = T::Real>,
+    ) -> T {
+        self.scalar
+            * self.integrand.evaluate_nonsingular(
+                test_table,
+                trial_table,
+                test_point_index,
+                trial_point_index,
+                test_basis_index,
+                trial_basis_index,
+                k,
+                test_geometry,
+                trial_geometry,
+            )
+    }
+
+    unsafe fn evaluate_singular(
+        &self,
+        test_table: &RlstArray<T, 4>,
+        trial_table: &RlstArray<T, 4>,
+        point_index: usize,
+        test_basis_index: usize,
+        trial_basis_index: usize,
+        k: &RlstArray<Self::T, 2>,
+        test_geometry: &impl CellGeometry<T = T::Real>,
+        trial_geometry: &impl CellGeometry<T = T::Real>,
+    ) -> T {
+        self.scalar
+            * self.integrand.evaluate_singular(
+                test_table,
+                trial_table,
+                point_index,
+                test_basis_index,
+                trial_basis_index,
+                k,
+                test_geometry,
+                trial_geometry,
+            )
+    }
+}
diff --git a/src/assembly/boundary/integrands/hypersingular.rs b/src/assembly/boundary/integrands/hypersingular.rs
@@ -131,19 +131,23 @@ impl<T: RlstScalar> BoundaryIntegrand for HypersingularCurlCurlBoundaryIntegrand
 }
 
 pub struct HypersingularNormalNormalBoundaryIntegrand<T: RlstScalar> {
-    wavenumber: T::Real,
     _t: std::marker::PhantomData<T>,
 }
 
 impl<T: RlstScalar> HypersingularNormalNormalBoundaryIntegrand<T> {
-    pub fn new(wavenumber: T::Real) -> Self {
+    pub fn new() -> Self {
         Self {
-            wavenumber,
             _t: std::marker::PhantomData,
         }
     }
 }
 
+impl<T: RlstScalar> Default for HypersingularNormalNormalBoundaryIntegrand<T> {
+    fn default() -> Self {
+        Self::new()
+    }
+}
+
 impl<T: RlstScalar> BoundaryIntegrand for HypersingularNormalNormalBoundaryIntegrand<T> {
     type T = T;
 
@@ -159,21 +163,20 @@ impl<T: RlstScalar> BoundaryIntegrand for HypersingularNormalNormalBoundaryInteg
         test_geometry: &impl CellGeometry<T = T::Real>,
         trial_geometry: &impl CellGeometry<T = T::Real>,
     ) -> T {
-        -*k.get_unchecked([0, test_point_index, trial_point_index])
+        *k.get_unchecked([0, test_point_index, trial_point_index])
             * num::cast::<T::Real, T>(
-                self.wavenumber.powi(2)
-                    * (*trial_geometry
+                *trial_geometry
+                    .normals()
+                    .get_unchecked([0, trial_point_index])
+                    * *test_geometry.normals().get_unchecked([0, test_point_index])
+                    + *trial_geometry
+                        .normals()
+                        .get_unchecked([1, trial_point_index])
+                        * *test_geometry.normals().get_unchecked([1, test_point_index])
+                    + *trial_geometry
                         .normals()
-                        .get_unchecked([0, trial_point_index])
-                        * *test_geometry.normals().get_unchecked([0, test_point_index])
-                        + *trial_geometry
-                            .normals()
-                            .get_unchecked([1, trial_point_index])
-                            * *test_geometry.normals().get_unchecked([1, test_point_index])
-                        + *trial_geometry
-                            .normals()
-                            .get_unchecked([2, trial_point_index])
-                            * *test_geometry.normals().get_unchecked([2, test_point_index])),
+                        .get_unchecked([2, trial_point_index])
+                        * *test_geometry.normals().get_unchecked([2, test_point_index]),
             )
             .unwrap()
             * *test_table.get_unchecked([0, test_point_index, test_basis_index, 0])
@@ -191,15 +194,14 @@ impl<T: RlstScalar> BoundaryIntegrand for HypersingularNormalNormalBoundaryInteg
         test_geometry: &impl CellGeometry<T = T::Real>,
         trial_geometry: &impl CellGeometry<T = T::Real>,
     ) -> T {
-        -*k.get_unchecked([0, point_index])
+        *k.get_unchecked([0, point_index])
             * num::cast::<T::Real, T>(
-                self.wavenumber.powi(2)
-                    * (*trial_geometry.normals().get_unchecked([0, point_index])
-                        * *test_geometry.normals().get_unchecked([0, point_index])
-                        + *trial_geometry.normals().get_unchecked([1, point_index])
-                            * *test_geometry.normals().get_unchecked([1, point_index])
-                        + *trial_geometry.normals().get_unchecked([2, point_index])
-                            * *test_geometry.normals().get_unchecked([2, point_index])),
+                *trial_geometry.normals().get_unchecked([0, point_index])
+                    * *test_geometry.normals().get_unchecked([0, point_index])
+                    + *trial_geometry.normals().get_unchecked([1, point_index])
+                        * *test_geometry.normals().get_unchecked([1, point_index])
+                    + *trial_geometry.normals().get_unchecked([2, point_index])
+                        * *test_geometry.normals().get_unchecked([2, point_index]),
             )
             .unwrap()
             * *test_table.get_unchecked([0, point_index, test_basis_index, 0])