Non-uniform current stimuli can be specified easier

examples/wrapped_torus.rs

@@ -1,6 +1,6 @@
 use coremem::{Driver, Flt, mat, meas};
 use coremem::geom::{Coord, Index, Meters, Sphere, Torus, Vec2, Vec3};
-use coremem::stim::{Stimulus, Sinusoid};
+use coremem::stim::{CurlStimulus, Sinusoid};
 use log::trace;
 
 fn main() {
@@ -29,7 +29,7 @@ fn main() {
     let size_px = Index((width_px, width_px, depth_px).into());
     let mut driver = Driver::new(size_px, feat_size);
     driver.set_steps_per_frame(120);
-    let base = "wrapped_torus-3";
+    let base = "wrapped_torus-4";
     let ferro_region = Torus::new_xy(Meters((half_width, half_width, half_depth).into()), ferro_major, ferro_minor);
     let drive_region = Torus::new_xz(Meters((half_width - ferro_major, half_width, half_depth).into()), wire_major, wire_minor);
     let sense_region = Torus::new_xz(Meters((half_width + ferro_major, half_width, half_depth).into()), wire_major, wire_minor);
@@ -43,19 +43,15 @@ fn main() {
     println!("boundary: {}um; {}um", m_to_um(boundary_xy), m_to_um(boundary_z));
     driver.add_upml_boundary(Meters((boundary_xy, boundary_xy, boundary_z).into()));
 
-    driver.add_stimulus(Stimulus::new(
+    driver.add_stimulus(CurlStimulus::new(
-        Sphere::new(
+        drive_region.clone(),
-            Meters((half_width - ferro_major + wire_major, half_width, half_depth).into()),
+        Sinusoid::from_wavelength(Vec3::new(0.0, 0.0, peak_current), current_duration * 2.0),
-            wire_minor),
+        drive_region.center(),
-        Sinusoid::from_wavelength(Vec3::new(0.0, 0.0, peak_current), current_duration * 2.0)
+        drive_region.axis(),
     ));
+
     driver.add_measurement(meas::Current::new("sense", sense_region.clone()));
     driver.add_measurement(meas::Current::new("drive", drive_region.clone()));
-    // TODO: simulate current in the drive element
-    //driver.add_stimulus(Stimulus::new(
-    //    conductor_region.clone(),
-    //    Sinusoid::from_wavelength(Vec3::new(0.0, 0.0, peak_current), current_duration * 2.0
-    //)));
 
     let prefix = format!("{}/{}-flt{}-{}-feat{}um-{}mA-{}ps--radii{}um-{}um-{}um-{}um",
         base,

src/geom/mod.rs

@@ -244,6 +244,12 @@ impl Torus {
     pub fn new_xz(center: Meters, major_rad: Flt, minor_rad: Flt) -> Self {
         Self::new(center, Meters(Vec3::new(0.0, 1.0, 0.0)), major_rad, minor_rad)
     }
+    pub fn center(&self) -> Meters {
+        self.center
+    }
+    pub fn axis(&self) -> Meters {
+        self.normal
+    }
 }
 
 #[typetag::serde]
@@ -285,6 +291,17 @@ impl Region for Sphere {
     }
 }
 
+/// Region describing the entire simulation space
+#[derive(Copy, Clone, Serialize, Deserialize)]
+pub struct WorldRegion;
+
+#[typetag::serde]
+impl Region for WorldRegion {
+    fn contains(&self, _: Meters) -> bool {
+        true
+    }
+}
+
 #[cfg(test)]
 mod test {
     use super::*;

src/geom/vec.rs

@@ -162,6 +162,16 @@ impl Vec3 {
         self.elem_mul(other).component_sum()
     }
 
+    pub fn cross(&self, other: Self) -> Self {
+        // det(   i    j    k  |
+        //    |  s.x  s.y  s.z |
+        //    |  o.x  o.y  o.z )
+        let x = self.y * other.z - other.y * self.z;
+        let y = other.z * self.z - self.x * other.z;
+        let z = self.x * other.y - other.x * self.y;
+        Self { x, y, z }
+    }
+
     /// Perform element-wise multiplication with `other`.
     pub fn elem_mul(&self, other: Self) -> Self {
         Self {

src/sim.rs

@@ -1,5 +1,5 @@
 use crate::{flt::{Flt, Real}, consts};
-use crate::geom::{Coord, Index, Meters, Vec3, Vec3u};
+use crate::geom::{Coord, Index, Meters, Region, Vec3, Vec3u};
 use crate::mat::{self, GenericMaterial, Material};
 use dyn_clone::{self, DynClone};
 use log::trace;
@@ -54,11 +54,19 @@ impl<'a> dyn GenericSim + 'a {
         self.sample(at.to_meters(self.feature_size()))
     }
     /// Apply `F` to each Cell, and sum the results.
-    pub fn map_sum<F: Fn(&Cell<mat::Static>) -> R + , R: Sum<R>>(&self, f: F) -> R {
+    pub fn map_sum<F, Ret>(&self, f: F) -> Ret
+    where
+          F: Fn(&Cell<mat::Static>) -> Ret,
+          Ret: Sum<Ret>,
+    {
         self.map_sum_enumerated(|_at: Index, cell| f(cell))
     }
 
-    pub fn map_sum_enumerated<C: Coord, F: Fn(C, &Cell<mat::Static>) -> R, R: Sum<R>>(&self, f: F) -> R {
+    pub fn map_sum_enumerated<C, F, Ret>(&self, f: F) -> Ret
+    where C: Coord,
+          F: Fn(C, &Cell<mat::Static>) -> Ret,
+          Ret: Sum<Ret>
+    {
         let (w, h, d) = (self.width(), self.height(), self.depth());
         (0..d).map(|z| (0..h).map(|y| (0..w).map(|x| {
             let at = Index(Vec3u::new(x, y, z));
@@ -66,6 +74,38 @@ impl<'a> dyn GenericSim + 'a {
         }).sum()).sum()).sum()
     }
 
+    pub fn volume_of_region<R: Region>(&self, region: &R) -> u32 {
+        self.map_sum_over(region, |_| 1)
+    }
+
+    /// Apply `F` to each Cell, and sum the results.
+    pub fn map_sum_over<F, Ret, Reg>(&self, region: &Reg, f: F) -> Ret
+    where
+          F: Fn(&Cell<mat::Static>) -> Ret,
+          Ret: Sum<Ret> + Default,
+          Reg: Region
+    {
+        self.map_sum_over_enumerated(region, |_at: Index, cell| f(cell))
+    }
+
+    pub fn map_sum_over_enumerated<C, F, Ret, Reg>(&self, region: &Reg, f: F) -> Ret
+    where C: Coord,
+          F: Fn(C, &Cell<mat::Static>) -> Ret,
+          Ret: Sum<Ret> + Default,
+          Reg: Region,
+    {
+        let (w, h, d) = (self.width(), self.height(), self.depth());
+        (0..d).map(|z| (0..h).map(|y| (0..w).map(|x| {
+            let at = Index(Vec3u::new(x, y, z));
+            let meters = at.to_meters(self.feature_size());
+            if region.contains(meters) {
+                f(C::from_index(at, self.feature_size()), &self.get(at))
+            } else {
+                Default::default()
+            }
+        }).sum()).sum()).sum()
+    }
+
     pub fn current<C: Coord>(&self, c: C) -> Vec3 {
         self.get(c).current_density() * (self.feature_size() * self.feature_size())
     }

src/stim.rs

@@ -1,6 +1,6 @@
 use crate::consts;
 use crate::flt::Flt;
-use crate::geom::{Coord as _, Index, Region, Vec3, Vec3u};
+use crate::geom::{Coord as _, Index, Meters, Region, Vec3, Vec3u};
 use crate::sim::GenericSim;
 
 use log::debug;
@@ -14,23 +14,23 @@ pub trait TimeVarying {
     fn at(&mut self, t_sec: Flt) -> Vec3;
 }
 
+/// Apply a time-varying stimulus uniformly across some region
 pub struct Stimulus<R, T> {
     region: R,
-    time: T,
+    stim: T,
 }
 
 impl<R, T> Stimulus<R, T> {
-    pub fn new(region: R, time: T) -> Self {
+    pub fn new(region: R, stim: T) -> Self {
         Self {
-            region, time
+            region, stim
         }
     }
 }
 
 impl<R: Region, T: TimeVarying> AbstractStimulus for Stimulus<R, T> {
     fn apply(&mut self, sim: &mut dyn GenericSim) {
-        // TODO: it would be nice to divide this amount across the # of cells matched
+        let amount = self.stim.at(sim.time());
-        let amount = self.time.at(sim.time());
         debug!("stim: {:?}", amount);
         for z in 0..sim.depth() {
             for y in 0..sim.height() {
@@ -45,6 +45,45 @@ impl<R: Region, T: TimeVarying> AbstractStimulus for Stimulus<R, T> {
     }
 }
 
+/// Apply a time-varying stimulus across some region.
+/// The stimulus seen at each point is based on its angle about the specified ray.
+pub struct CurlStimulus<R, T> {
+    region: R,
+    stim: T,
+    center: Meters,
+    axis: Meters,
+}
+
+impl<R, T> CurlStimulus<R, T> {
+    pub fn new(region: R, stim: T, center: Meters, axis: Meters) -> Self {
+        Self { region, stim, center, axis }
+    }
+}
+
+impl<R: Region, T: TimeVarying> AbstractStimulus for CurlStimulus<R, T> {
+    fn apply(&mut self, sim: &mut dyn GenericSim) {
+        let amount = self.stim.at(sim.time());
+        debug!("stim: {:?}", amount);
+        // TODO: should be possible to lift more of this into the sim to parallelize it 
+        for z in 0..sim.depth() {
+            for y in 0..sim.height() {
+                for x in 0..sim.width() {
+                    let loc = Index(Vec3u::new(x, y, z));
+                    let meters = loc.to_meters(sim.feature_size());
+                    if self.region.contains(meters) {
+                        let from_center_to_point = *meters - *self.center;
+                        let rotational = from_center_to_point.cross(*self.axis);
+                        let impulse = rotational.with_mag(amount.mag());
+                        // TODO: should somehow preserve the *components* of the time-varying
+                        // thing, not just the magnitude.
+                        sim.impulse_e(loc, impulse);
+                    }
+                }
+            }
+        }
+    }
+}
+
 pub struct Sinusoid {
     amp: Vec3,
     omega: Flt,