From a2c30dd6b68b8e364f9c880960b5576235ebf5ec Mon Sep 17 00:00:00 2001
From: Matthew Scroggs <matthew.w.scroggs@gmail.com>
Date: Wed, 28 Aug 2024 17:00:34 +0100
Subject: [PATCH] store test_integrand rather than remultiplying every time
 (#292)

---
 Cargo.toml                                    |   3 +-
 src/assembly/boundary/cell_pair_assemblers.rs | 103 +++++++++---------
 2 files changed, 53 insertions(+), 53 deletions(-)

diff --git a/Cargo.toml b/Cargo.toml
index e65e0b20..d8b3f05f 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -25,7 +25,8 @@ mpi = { version = "0.8.*", optional = true }
 num = "0.4"
 lazy_static = "1.4"
 ndelement = "0.1.0"
-ndgrid = "0.1.0"
+# ndgrid = "0.1.0"
+ndgrid = { git = "https://github.com/bempp/ndgrid.git" }
 rayon = "1.9"
 rlst = "0.2"
 green-kernels = "0.2.0"
diff --git a/src/assembly/boundary/cell_pair_assemblers.rs b/src/assembly/boundary/cell_pair_assemblers.rs
index c60d4605..187b29fc 100644
--- a/src/assembly/boundary/cell_pair_assemblers.rs
+++ b/src/assembly/boundary/cell_pair_assemblers.rs
@@ -137,8 +137,8 @@ impl<
             for (test_i, entry) in col.iter_mut().enumerate() {
                 *entry = T::zero();
                 for (index, wt) in self.weights.iter().enumerate() {
-                    unsafe {
-                        *entry += self.integrand.evaluate_singular(
+                    *entry += unsafe {
+                        self.integrand.evaluate_singular(
                             self.test_table,
                             self.trial_table,
                             index,
@@ -147,12 +147,11 @@ impl<
                             &self.k,
                             &test_geometry,
                             &trial_geometry,
-                        ) * num::cast::<T::Real, T>(
-                            *wt * *self.test_jdet.get_unchecked(index)
-                                * *self.trial_jdet.get_unchecked(index),
                         )
-                        .unwrap();
-                    }
+                    } * num::cast::<T::Real, T>(
+                        *wt * self.test_jdet[index] * self.trial_jdet[index],
+                    )
+                    .unwrap();
                 }
             }
         }
@@ -292,26 +291,25 @@ impl<
             for (test_i, entry) in col.iter_mut().enumerate() {
                 *entry = T::zero();
                 for (test_index, test_wt) in self.test_weights.iter().enumerate() {
+                    let test_integrand =
+                        num::cast::<T::Real, T>(*test_wt * self.test_jdet[test_index]).unwrap();
                     for (trial_index, trial_wt) in self.trial_weights.iter().enumerate() {
-                        *entry += unsafe {
-                            self.integrand.evaluate_nonsingular(
-                                self.test_table,
-                                self.trial_table,
-                                test_index,
-                                trial_index,
-                                test_i,
-                                trial_i,
-                                &self.k,
-                                &test_geometry,
-                                &trial_geometry,
-                            ) * num::cast::<T::Real, T>(
-                                *test_wt
-                                    * self.test_jdet[test_index]
-                                    * *trial_wt
-                                    * self.trial_jdet[trial_index],
-                            )
-                            .unwrap()
-                        };
+                        *entry +=
+                            unsafe {
+                                self.integrand.evaluate_nonsingular(
+                                    self.test_table,
+                                    self.trial_table,
+                                    test_index,
+                                    trial_index,
+                                    test_i,
+                                    trial_i,
+                                    &self.k,
+                                    &test_geometry,
+                                    &trial_geometry,
+                                )
+                            } * num::cast::<T::Real, T>(*trial_wt * self.trial_jdet[trial_index])
+                                .unwrap()
+                                * test_integrand;
                     }
                 }
             }
@@ -474,6 +472,8 @@ impl<
             for (test_i, entry) in col.iter_mut().enumerate() {
                 *entry = T::zero();
                 for (test_index, test_wt) in self.test_weights.iter().enumerate() {
+                    let test_integrand =
+                        num::cast::<T::Real, T>(*test_wt * self.test_jdet[test_index]).unwrap();
                     for (trial_index, trial_wt) in self.trial_weights.iter().enumerate() {
                         *entry += unsafe {
                             self.integrand.evaluate_nonsingular(
@@ -486,14 +486,12 @@ impl<
                                 &self.k,
                                 &test_geometry,
                                 &trial_geometry,
-                            ) * num::cast::<T::Real, T>(
-                                *test_wt
-                                    * self.test_jdet[test_index]
-                                    * *trial_wt
-                                    * self.trial_jdet[self.trial_cell][trial_index],
                             )
-                            .unwrap()
-                        };
+                        } * num::cast::<T::Real, T>(
+                            *trial_wt * self.trial_jdet[self.trial_cell][trial_index],
+                        )
+                        .unwrap()
+                            * test_integrand;
                     }
                 }
             }
@@ -654,26 +652,27 @@ impl<
             for (test_i, entry) in col.iter_mut().enumerate() {
                 *entry = T::zero();
                 for (test_index, test_wt) in self.test_weights.iter().enumerate() {
+                    let test_integrand = num::cast::<T::Real, T>(
+                        *test_wt * self.test_jdet[self.test_cell][test_index],
+                    )
+                    .unwrap();
                     for (trial_index, trial_wt) in self.trial_weights.iter().enumerate() {
-                        *entry += unsafe {
-                            self.integrand.evaluate_nonsingular(
-                                self.test_table,
-                                self.trial_table,
-                                test_index,
-                                trial_index,
-                                test_i,
-                                trial_i,
-                                &self.k,
-                                &test_geometry,
-                                &trial_geometry,
-                            ) * num::cast::<T::Real, T>(
-                                *test_wt
-                                    * self.test_jdet[self.test_cell][test_index]
-                                    * *trial_wt
-                                    * self.trial_jdet[trial_index],
-                            )
-                            .unwrap()
-                        };
+                        *entry +=
+                            unsafe {
+                                self.integrand.evaluate_nonsingular(
+                                    self.test_table,
+                                    self.trial_table,
+                                    test_index,
+                                    trial_index,
+                                    test_i,
+                                    trial_i,
+                                    &self.k,
+                                    &test_geometry,
+                                    &trial_geometry,
+                                )
+                            } * num::cast::<T::Real, T>(*trial_wt * self.trial_jdet[trial_index])
+                                .unwrap()
+                                * test_integrand;
                     }
                 }
             }