CurrentLoop: use a better justified measurement algorithm

'course the best way to justify it is with tests: hopefully those will
come shortly.

crates/coremem/src/meas.rs

@@ -3,6 +3,7 @@ use crate::real::{Real as _, ToFloat as _};
 use crate::cross::vec::{Vec3, Vec3u};
 use crate::sim::AbstractSim;
 use serde::{Serialize, Deserialize};
+use std::ops::AddAssign;
 
 // TODO: do we really need both Send and Sync?
 pub trait AbstractMeasurement<S>: Send + Sync {
@@ -243,6 +244,24 @@ impl Current {
     }
 }
 
+// TODO: clean up these FieldSample types
+#[derive(Default)]
+struct TupleSum<T>(T);
+
+impl<T0: Default + AddAssign, T1: Default + AddAssign> std::iter::Sum for TupleSum<(T0, T1)> {
+    fn sum<I>(iter: I) -> Self
+        where I: Iterator<Item = Self>
+    {
+        let mut s = Self::default();
+        for TupleSum((a0, a1)) in iter {
+            s.0.0 += a0;
+            s.0.1 += a1;
+        }
+        s
+    }
+}
+
+
 #[derive(Default)]
 struct FieldSample(u32, f64, Vec3<f64>);
 
@@ -323,33 +342,30 @@ impl<R> CurrentLoop<R> {
 }
 impl<R: Region + HasCrossSection> CurrentLoop<R> {
     fn data<S: AbstractSim>(&self, state: &S) -> f32 {
-        // i use a statistical lens for this:
-        // 1. current is the rate of flow of charge into a surface.
-        // 2. in any context where it makes sense to think of current, the current through each
-        //    cross-sectional **is the same**.
-        // 3. each point in our 3d region belongs to exactly one cross-sectional surface.
-        // 4. so, given a point: what's the expected current through the cross section it belongs to?
-        //    - answer: that point's current density times the cross section's area.
-        // 5. average the above over the whole volume, and you get an "average current".
-        //
-        // we're sampling uniformly over the cell space -- not the set of cross sections.
-        // - however, if all cross sections have equal area, this is equivalent.
-        // sampling all points (instead of just a single point):
-        //   1) removes bias from step #4: current *within* a cross section is not uniform, but if
-        //      we sample every point within the cross section and weight them equally, then the
-        //      average is the truth.
-        //   2) probably combats grid quantization / artifacting.
-        let FieldSample(num_samples, sum_cross_sectional_current, _current_vec) = state.map_sum_over_enumerated(&self.region, |coord: Meters, _cell| {
+        // - current exists as a property of a 2d surface.
+        // - the user has provided us a 3d volume which behaves as though it's an extruded surface:
+        //   for any point in the volume we can query the normal vector of the cross section
+        //   containing that point.
+        // - we choose that measuring the "current" on such a volume means to measure the average
+        //   current through all its cross sections (for most boring materials, each
+        //   cross section has nearly identical current).
+        // - therefore, enumerate the entire volume and compute the "net" current (the sum over
+        //   each cell of whatever current in that cell is along the cross-section normal).
+        //   then divide by the number of complete cross sections we measured, to average.
+        let feature_area = state.feature_size() * state.feature_size();
+        let TupleSum((net_current, cross_sections)) = state.map_sum_over_enumerated(&self.region, move |coord: Meters, _cell| {
             // `normal` represents both the size of the cross section (m^2) this cell belongs to,
             // and the normal direction of the cross section.
             let normal = self.region.cross_section_normal(coord);  // [m^2]
+            // calculate the amount of normal current through this specific cell
             let current_density = state.current_density(coord);    // [A/m^2]
-            // now we have an estimation of the entire current flowing through the cross section
-            // this cell belongs to.
-            let cross_sectional_current = current_density.dot(normal.cast()); // [A]
-            FieldSample(1, cross_sectional_current.cast(), current_density.cast())
+            let cross_sectional_current = feature_area * current_density.dot(normal.norm()); // [A]
+            // keep track of how many cross sections we enumerate, since each additional cross
+            // sections represents a double-count of the current.
+            let num_cross_sections_filled = feature_area / normal.mag();
+            TupleSum((cross_sectional_current, num_cross_sections_filled))
         });
-        let mean_cross_sectional_current = sum_cross_sectional_current.cast::<f32>() / f32::from_primitive(num_samples);
+        let mean_cross_sectional_current = net_current.cast::<f32>() / cross_sections;
         mean_cross_sectional_current
     }
 }

crates/cross/src/real.rs

@@ -107,6 +107,10 @@ pub trait Real:
         self == Self::zero()
     }
 
+    fn inv(self) -> Self {
+        Self::one() / self
+    }
+
     fn zero() -> Self;
     fn one() -> Self;
     fn two() -> Self;