diff --git a/src/driver.rs b/src/driver.rs index 5ad85f8..0d0f37a 100644 --- a/src/driver.rs +++ b/src/driver.rs @@ -206,9 +206,7 @@ impl Driver { { trace!("stimuli begin"); let start_time = Instant::now(); - for stim in &mut *self.stimuli { - stim.apply(&mut self.state); - } + self.state.apply_stimulus(&self.stimuli); self.time_spent_on_stimuli += start_time.elapsed(); } diff --git a/src/sim.rs b/src/sim.rs index 16b51c9..7213d16 100644 --- a/src/sim.rs +++ b/src/sim.rs @@ -1,6 +1,7 @@ use crate::{flt::{Flt, Real}, consts}; use crate::geom::{Coord, Index, Meters, Region, Vec3, Vec3u}; use crate::mat::{self, GenericMaterial, Material}; +use crate::stim::AbstractStimulus; use dyn_clone::{self, DynClone}; use log::trace; @@ -191,6 +192,20 @@ impl SimState { self.step_no += 1; } + + pub fn apply_stimulus(&mut self, stim: &S) { + let feature_size = self.feature_size(); + + let t_sec = self.time(); + trace!("apply_stimulus begin"); + Zip::from(ndarray::indices_of(&self.cells)).and(&mut self.cells).par_apply( + |(z, y, x), cell| { + let pos_meters = Index((x as u32, y as u32, z as u32).into()).to_meters(feature_size); + let value = stim.at(t_sec, pos_meters); + cell.state.e += value; + }); + trace!("apply_stimulus end"); + } } impl GenericSim for SimState { diff --git a/src/stim.rs b/src/stim.rs index 7507d77..a72b20d 100644 --- a/src/stim.rs +++ b/src/stim.rs @@ -5,13 +5,27 @@ use crate::sim::GenericSim; use log::debug; -pub trait AbstractStimulus { - fn apply(&mut self, sim: &mut dyn GenericSim); +pub trait AbstractStimulus: Sync { + /// Return the E field which should be added to the provided position/time. + // TODO: this needs to be made independent of the time step. + fn at(&self, t_sec: Flt, pos: Meters) -> Vec3; +} + +impl AbstractStimulus for Vec { + fn at(&self, t_sec: Flt, pos: Meters) -> Vec3 { + self.iter().map(|s| s.at(t_sec, pos)).sum() + } +} + +impl AbstractStimulus for Box { + fn at(&self, t_sec: Flt, pos: Meters) -> Vec3 { + (**self).at(t_sec, pos) + } } pub trait TimeVarying { /// Retrieve the E impulse to apply at the provided time (in seconds). - fn at(&mut self, t_sec: Flt) -> Vec3; + fn at(&self, t_sec: Flt) -> Vec3; } /// Apply a time-varying stimulus uniformly across some region @@ -28,19 +42,12 @@ impl Stimulus { } } -impl AbstractStimulus for Stimulus { - fn apply(&mut self, sim: &mut dyn GenericSim) { - let amount = self.stim.at(sim.time()); - debug!("stim: {:?}", amount); - 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)); - if self.region.contains(loc.to_meters(sim.feature_size())) { - sim.impulse_e(loc, amount); - } - } - } +impl AbstractStimulus for Stimulus { + fn at(&self, t_sec: Flt, pos: Meters) -> Vec3 { + if self.region.contains(pos) { + self.stim.at(t_sec) + } else { + Vec3::zero() } } } @@ -60,26 +67,18 @@ impl CurlStimulus { } } -impl AbstractStimulus for CurlStimulus { - 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); - } - } - } +impl AbstractStimulus for CurlStimulus { + fn at(&self, t_sec: Flt, pos: Meters) -> Vec3 { + if self.region.contains(pos) { + let amount = self.stim.at(t_sec); + let from_center_to_point = *pos - *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. + impulse + } else { + Vec3::zero() } } } @@ -102,7 +101,7 @@ impl Sinusoid { } impl TimeVarying for Sinusoid { - fn at(&mut self, t_sec: Flt) -> Vec3 { + fn at(&self, t_sec: Flt) -> Vec3 { self.amp * (t_sec * self.omega).sin() } }