From 135a39e75eb3ed2e99ed1c96db2694c717d7bca7 Mon Sep 17 00:00:00 2001
From: Timo Betcke <timo.betcke@gmail.com>
Date: Thu, 31 Oct 2024 21:37:45 +0000
Subject: [PATCH] WIP: Some cleanup

---
 src/assembly/boundary/assemblers.rs           | 53 ++++++++-----------
 src/assembly/boundary/cell_pair_assemblers.rs | 16 +++---
 2 files changed, 29 insertions(+), 40 deletions(-)

diff --git a/src/assembly/boundary/assemblers.rs b/src/assembly/boundary/assemblers.rs
index a48090ff..90ac18fe 100644
--- a/src/assembly/boundary/assemblers.rs
+++ b/src/assembly/boundary/assemblers.rs
@@ -20,7 +20,7 @@ use bempp_quadrature::duffy::{
 };
 use bempp_quadrature::simplex_rules::simplex_rule;
 use bempp_quadrature::types::{CellToCellConnectivity, TestTrialNumericalQuadratureDefinition};
-use itertools::izip;
+use itertools::{izip, Itertools};
 #[cfg(feature = "mpi")]
 use mpi::traits::Communicator;
 use ndelement::reference_cell;
@@ -853,47 +853,36 @@ impl<
 
         for test_cell_type in test_space.grid().entity_types(2) {
             let npts_test = self.options.quadrature_degrees[test_cell_type];
+            let qrule_test = simplex_rule(*test_cell_type, npts_test).unwrap();
+            let mut qpoints_test = rlst_dynamic_array2!(<T as RlstScalar>::Real, [2, npts_test]);
+            for (target, source) in izip!(qpoints_test.data_mut(), qrule_test.points.iter()) {
+                *target = num::cast::<f64, <T as RlstScalar>::Real>(*source).unwrap();
+            }
+            let qweights_test = qrule_test
+                .weights
+                .iter()
+                .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!(
+                T,
+                test_element.tabulate_array_shape(self.table_derivs, npts_test)
+            );
+            test_element.tabulate(&qpoints_test, self.table_derivs, &mut test_table);
+
             for trial_cell_type in trial_space.grid().entity_types(2) {
                 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!(<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, <T as RlstScalar>::Real>(qrule_test.points[2 * i + j])
-                                .unwrap();
-                    }
-                }
-                let qweights_test = qrule_test
-                    .weights
-                    .iter()
-                    .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!(<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, <T as RlstScalar>::Real>(
-                                qrule_trial.points[2 * i + j],
-                            )
-                            .unwrap();
-                    }
+                for (target, source) in izip!(qpoints_trial.data_mut(), qrule_trial.points.iter()) {
+                    *target = num::cast::<f64, <T as RlstScalar>::Real>(*source).unwrap();
                 }
                 let qweights_trial = qrule_trial
                     .weights
                     .iter()
                     .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!(
-                    T,
-                    test_element.tabulate_array_shape(self.table_derivs, npts_test)
-                );
-                test_element.tabulate(&qpoints_test, self.table_derivs, &mut test_table);
+                    .collect_vec();
 
                 let trial_element = trial_space.element(*trial_cell_type);
                 let mut trial_table = rlst_dynamic_array4!(
diff --git a/src/assembly/boundary/cell_pair_assemblers.rs b/src/assembly/boundary/cell_pair_assemblers.rs
index 5391dcee..734e598b 100644
--- a/src/assembly/boundary/cell_pair_assemblers.rs
+++ b/src/assembly/boundary/cell_pair_assemblers.rs
@@ -136,7 +136,13 @@ impl<
         for (trial_i, mut col) in local_mat.col_iter_mut().enumerate() {
             for (test_i, entry) in col.iter_mut().enumerate() {
                 *entry = T::zero();
-                for (index, wt) in self.weights.iter().enumerate() {
+                for (index, (&wt, &test_det, &trial_det)) in izip!(
+                    self.weights.iter(),
+                    self.test_jdet.iter(),
+                    self.trial_jdet.iter()
+                )
+                .enumerate()
+                {
                     *entry += self.integrand.evaluate_singular(
                         self.test_table,
                         self.trial_table,
@@ -146,13 +152,7 @@ impl<
                         &self.k,
                         &test_geometry,
                         &trial_geometry,
-                    ) * num::cast::<T::Real, T>(
-                        *wt * unsafe {
-                            *self.test_jdet.get_unchecked(index)
-                                * *self.trial_jdet.get_unchecked(index)
-                        },
-                    )
-                    .unwrap();
+                    ) * num::cast::<T::Real, T>(wt * test_det * trial_det).unwrap();
                 }
             }
         }