Merge branch 'main' into mscroggs/unsafe-trait

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),
         )

src/assembly/boundary/integrands.rs

@@ -108,3 +108,75 @@ unsafe impl<T: RlstScalar, I0: BoundaryIntegrand<T = T>, I1: BoundaryIntegrand<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 }
+    }
+}
+
+unsafe impl<T: RlstScalar, I: BoundaryIntegrand<T = T>> BoundaryIntegrand
+    for BoundaryIntegrandScalarProduct<T, I>
+{
+    type T = T;
+
+    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 {
+        unsafe {
+            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,
+                )
+        }
+    }
+
+    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 {
+        unsafe {
+            self.scalar
+                * self.integrand.evaluate_singular(
+                    test_table,
+                    trial_table,
+                    point_index,
+                    test_basis_index,
+                    trial_basis_index,
+                    k,
+                    test_geometry,
+                    trial_geometry,
+                )
+        }
+    }
+}

src/assembly/boundary/integrands/hypersingular.rs

@@ -135,19 +135,23 @@ unsafe impl<T: RlstScalar> BoundaryIntegrand for HypersingularCurlCurlBoundaryIn
 }
 
 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()
+    }
+}
+
 unsafe impl<T: RlstScalar> BoundaryIntegrand for HypersingularNormalNormalBoundaryIntegrand<T> {
     type T = T;