diff --git a/examples/toroid25d.rs b/examples/toroid25d.rs index d5ca1d0..0d8867d 100644 --- a/examples/toroid25d.rs +++ b/examples/toroid25d.rs @@ -1,5 +1,5 @@ -use coremem::{consts, Driver, Flt, mat, meas}; -use coremem::geom::{Coord, CylinderZ, Vec2, Vec3}; +use coremem::{Driver, Flt, mat, meas}; +use coremem::geom::{Coord, CylinderZ, Index, Meters, Vec2, Vec3}; use coremem::stim::{Stimulus, Sinusoid}; fn main() { @@ -22,14 +22,14 @@ fn main() { let width_px = from_m(width); let depth_px = from_m(depth); - let size_px = (width_px, width_px, depth_px).into(); + 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(8); //driver.set_steps_per_frame(40); //driver.set_steps_per_frame(200); driver.add_y4m_renderer(&*format!("toroid25d-flt{}-{}-feat{}um-{}A--radii{}um-{}um-{}um.y4m", std::mem::size_of::() * 8, - size_px, + *size_px, m_to_um(feat_size), peak_current as u32, m_to_um(conductor_outer_rad), @@ -43,31 +43,31 @@ fn main() { driver.add_measurement(meas::Label(format!("Conductivity: {}, Imax: {:.2e}", conductivity, peak_current))); driver.add_measurement(meas::Current(conductor_region.clone())); driver.add_measurement(meas::Magnetization( - (half_width + ferro_inner_rad + 2.0*feat_size, half_width, half_depth).into() + Meters((half_width + ferro_inner_rad + 2.0*feat_size, half_width, half_depth).into()) )); driver.add_measurement(meas::MagneticFlux( - (half_width + ferro_inner_rad + 2.0*feat_size, half_width, half_depth).into() + Meters((half_width + ferro_inner_rad + 2.0*feat_size, half_width, half_depth).into()) )); driver.add_measurement(meas::MagneticStrength( - (half_width + ferro_inner_rad + 2.0*feat_size, half_width, half_depth).into() + Meters((half_width + ferro_inner_rad + 2.0*feat_size, half_width, half_depth).into()) )); driver.add_measurement(meas::Magnetization( - (half_width + ferro_inner_rad + 1.0*feat_size, half_width, half_depth).into() + Meters((half_width + ferro_inner_rad + 1.0*feat_size, half_width, half_depth).into()) )); driver.add_measurement(meas::MagneticFlux( - (half_width + ferro_inner_rad + 1.0*feat_size, half_width, half_depth).into() + Meters((half_width + ferro_inner_rad + 1.0*feat_size, half_width, half_depth).into()) )); driver.add_measurement(meas::MagneticStrength( - (half_width + ferro_inner_rad + 1.0*feat_size, half_width, half_depth).into() + Meters((half_width + ferro_inner_rad + 1.0*feat_size, half_width, half_depth).into()) )); driver.add_measurement(meas::Magnetization( - (half_width + 0.5 * (ferro_inner_rad + ferro_outer_rad), half_width, half_depth).into() + Meters((half_width + 0.5 * (ferro_inner_rad + ferro_outer_rad), half_width, half_depth).into()) )); driver.add_measurement(meas::MagneticFlux( - (half_width + 0.5 * (ferro_inner_rad + ferro_outer_rad), half_width, half_depth).into() + Meters((half_width + 0.5 * (ferro_inner_rad + ferro_outer_rad), half_width, half_depth).into()) )); driver.add_measurement(meas::MagneticStrength( - (half_width + 0.5 * (ferro_inner_rad + ferro_outer_rad), half_width, half_depth).into() + Meters((half_width + 0.5 * (ferro_inner_rad + ferro_outer_rad), half_width, half_depth).into()) )); let center = Vec2::new(half_width as _, half_width as _); @@ -75,7 +75,7 @@ fn main() { for z_px in 0..depth_px { for y_px in 0..width_px { for x_px in 0..width_px { - let loc = Coord::new(x_px, y_px, z_px); + let loc = Index((x_px, y_px, z_px).into()); let d = Vec2::new(to_m(x_px), to_m(y_px)) - center; let r = d.mag(); if (conductor_inner_rad as _..conductor_outer_rad as _).contains(&r) { @@ -90,7 +90,7 @@ fn main() { } let boundary_xy = half_width - ferro_outer_rad - buffer; println!("boundary: {}um", m_to_um(boundary_xy)); - let boundary = Coord::new(from_m(boundary_xy), from_m(boundary_xy), 20); + let boundary = Index((from_m(boundary_xy), from_m(boundary_xy), 20).into()); driver.add_upml_boundary(boundary); driver.add_stimulus(Stimulus::new( diff --git a/src/driver.rs b/src/driver.rs index d3dcc4c..bdd8fff 100644 --- a/src/driver.rs +++ b/src/driver.rs @@ -1,6 +1,6 @@ use crate::{flt::Flt, mat}; use crate::consts; -use crate::geom::{Coord, Vec3}; +use crate::geom::{Coord, Index, Vec3, Vec3u}; use crate::meas::{self, AbstractMeasurement}; use crate::render::{self, MultiRenderer, Renderer}; use crate::sim::{GenericSim as _, SimState}; @@ -23,9 +23,9 @@ pub struct Driver { } impl Driver { - pub fn new(size: Coord, feature_size: Flt) -> Self { + pub fn new(size: C, feature_size: Flt) -> Self { Driver { - state: SimState::new(size, feature_size), + state: SimState::new(size.to_index(feature_size), feature_size), renderer: Default::default(), steps_per_frame: 1, time_spent_stepping: Default::default(), @@ -83,8 +83,9 @@ impl Driver { // } // } - pub fn add_upml_boundary(&mut self, thickness: Coord) { + pub fn add_upml_boundary(&mut self, thickness: C) { // Based on explanation here (slide 63): https://empossible.net/wp-content/uploads/2020/01/Lecture-The-Perfectly-Matched-Layer.pdf + let thickness = thickness.to_index(self.state.feature_size()); let h = self.state.height(); let w = self.state.width(); let d = self.state.depth(); @@ -119,7 +120,7 @@ impl Driver { let cond_y = scale * (depth_y as Flt/thickness.y() as Flt).powf(3.0); let cond_z = scale * (depth_z as Flt/thickness.z() as Flt).powf(3.0); let conductor = mat::Static::anisotropic_conductor(Vec3::new(cond_x, cond_y, cond_z)); - *self.state.get_mut((x, y, z).into()).mat_mut() = conductor.into(); + *self.state.get_mut(Index(Vec3u::new(x, y, z))).mat_mut() = conductor.into(); } } } diff --git a/src/geom/mod.rs b/src/geom/mod.rs index 17d29ef..f5c9a6b 100644 --- a/src/geom/mod.rs +++ b/src/geom/mod.rs @@ -1,8 +1,10 @@ -mod coord; +mod units; mod vec; +mod vecu; -pub use coord::Coord; +pub use units::{Coord, Meters, Index}; pub use vec::{Vec2, Vec3}; +pub use vecu::Vec3u; use crate::flt::{Flt, Real}; use std::fmt::{self, Display}; @@ -180,7 +182,7 @@ impl Polygon2d { } pub trait Region { - fn contains(&self, p: Vec3) -> bool; + fn contains(&self, p: Meters) -> bool; } #[derive(Copy, Clone)] @@ -190,6 +192,7 @@ pub struct CylinderZ { } impl CylinderZ { + // TODO: should be unit-typed pub fn new(center: Vec2, radius: Flt) -> Self { Self { center, @@ -199,7 +202,7 @@ impl CylinderZ { } impl Region for CylinderZ { - fn contains(&self, p: Vec3) -> bool { + fn contains(&self, p: Meters) -> bool { p.xy().distance_sq(self.center) <= (self.radius * self.radius).into() } } diff --git a/src/geom/units.rs b/src/geom/units.rs new file mode 100644 index 0000000..fd31d46 --- /dev/null +++ b/src/geom/units.rs @@ -0,0 +1,60 @@ +use crate::flt::Flt; +use super::{Vec3, Vec3u}; +use std::ops::Deref; + +pub trait Coord: Copy + Clone { + fn to_meters(&self, feature_size: Flt) -> Meters; + fn to_index(&self, feature_size: Flt) -> Index; + fn from_meters(other: Meters, feature_size: Flt) -> Self; + fn from_index(other: Index, feature_size: Flt) -> Self; +} + +#[derive(Copy, Clone)] +pub struct Meters(pub Vec3); + +impl Coord for Meters { + fn to_meters(&self, _feature_size: Flt) -> Meters { + *self + } + fn to_index(&self, feature_size: Flt) -> Index { + Index((self.0 / feature_size).into()) + } + fn from_meters(other: Meters, _feature_size: Flt) -> Self { + other + } + fn from_index(other: Index, feature_size: Flt) -> Self { + other.to_meters(feature_size) + } +} + +impl Deref for Meters { + type Target = Vec3; + fn deref(&self) -> &Vec3 { + &self.0 + } +} + +#[derive(Copy, Clone)] +pub struct Index(pub Vec3u); + +impl Coord for Index { + fn to_meters(&self, feature_size: Flt) -> Meters { + Meters(Vec3::from(self.0) * feature_size) + } + fn to_index(&self, _feature_size: Flt) -> Index { + *self + } + fn from_meters(other: Meters, feature_size: Flt) -> Self { + other.to_index(feature_size) + } + fn from_index(other: Index, _feature_size: Flt) -> Self { + other + } +} + +impl Deref for Index { + type Target = Vec3u; + fn deref(&self) -> &Vec3u { + &self.0 + } +} diff --git a/src/geom/vec.rs b/src/geom/vec.rs index a234229..fab4e83 100644 --- a/src/geom/vec.rs +++ b/src/geom/vec.rs @@ -1,5 +1,5 @@ use crate::flt::{Flt, Real}; -use super::Coord; +use super::Vec3u; use std::convert::From; use std::iter::Sum; use std::ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub}; @@ -197,9 +197,9 @@ impl From<(Real, Real, Real)> for Vec3 { } } -impl From for Vec3 { - fn from(c: Coord) -> Self { - Self::new(c.x().into(), c.y().into(), c.z().into()) +impl From for Vec3 { + fn from(v: Vec3u) -> Self { + Self::new(v.x().into(), v.y().into(), v.z().into()) } } diff --git a/src/geom/coord.rs b/src/geom/vecu.rs similarity index 76% rename from src/geom/coord.rs rename to src/geom/vecu.rs index 6820cd4..c1da694 100644 --- a/src/geom/coord.rs +++ b/src/geom/vecu.rs @@ -1,12 +1,13 @@ +use super::Vec3; use std::fmt::{self, Display}; #[derive(Copy, Clone, Default, Eq, PartialEq)] -pub struct Coord { +pub struct Vec3u { y: u32, x: u32, z: u32, } -impl Coord { +impl Vec3u { pub fn new(x: u32, y: u32, z: u32) -> Self { Self { x, @@ -29,13 +30,19 @@ impl Coord { } } -impl From<(u32, u32, u32)> for Coord { +impl From<(u32, u32, u32)> for Vec3u { fn from((x, y, z): (u32, u32, u32)) -> Self { Self::new(x, y, z) } } -impl Display for Coord { +impl From for Vec3u { + fn from(v: Vec3) -> Self { + Self::new(v.x() as _, v.y() as _, v.z() as _) + } +} + +impl Display for Vec3u { fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { write!(f, "({}, {}, {})", self.x, self.y, self.z) } diff --git a/src/meas.rs b/src/meas.rs index a3c48be..4a3a756 100644 --- a/src/meas.rs +++ b/src/meas.rs @@ -1,5 +1,5 @@ use crate::flt::Flt; -use crate::geom::{Coord, Vec3, Region}; +use crate::geom::{Meters, Region}; use crate::mat::Material as _; use crate::sim::{Cell, GenericSim}; use std::fmt::Display; @@ -39,12 +39,14 @@ impl AbstractMeasurement for Label { } } -fn sum_over_region, R: Region, F: Fn(Coord, &Cell) -> T>(state: &dyn GenericSim, r: &R, f: F) -> T { +fn sum_over_region, R: Region, F: Fn(Meters, &Cell) -> T>(state: &dyn GenericSim, r: &R, f: F) -> T { // TODO: z coordinate? - state.map_sum_enumerated(|coord, cell| if r.contains(Vec3::from(coord)*state.feature_size()) { - f(coord, cell) - } else { - Default::default() + state.map_sum_enumerated(|coord, cell| { + if r.contains(coord) { + f(coord, cell) + } else { + Default::default() + } }) } @@ -57,13 +59,13 @@ impl AbstractMeasurement for Current { } } -fn loc(v: Vec3) -> String { - let (x, y, z): (Flt, Flt, Flt) = (v * 1e6).into(); +fn loc(v: Meters) -> String { + let (x, y, z): (Flt, Flt, Flt) = (*v * 1e6).into(); format!("({}, {}, {}) um", x as i32, y as i32, z as i32) } /// M -pub struct Magnetization(pub Vec3); +pub struct Magnetization(pub Meters); impl AbstractMeasurement for Magnetization { fn eval(&self, state: &dyn GenericSim) -> String { @@ -73,7 +75,7 @@ impl AbstractMeasurement for Magnetization { } /// B -pub struct MagneticFlux(pub Vec3); +pub struct MagneticFlux(pub Meters); impl AbstractMeasurement for MagneticFlux { fn eval(&self, state: &dyn GenericSim) -> String { @@ -83,7 +85,7 @@ impl AbstractMeasurement for MagneticFlux { } /// H -pub struct MagneticStrength(pub Vec3); +pub struct MagneticStrength(pub Meters); impl AbstractMeasurement for MagneticStrength { fn eval(&self, state: &dyn GenericSim) -> String { @@ -92,7 +94,7 @@ impl AbstractMeasurement for MagneticStrength { } } -pub struct ElectricField(pub Vec3); +pub struct ElectricField(pub Meters); impl AbstractMeasurement for ElectricField { fn eval(&self, state: &dyn GenericSim) -> String { diff --git a/src/render.rs b/src/render.rs index 1fac9f0..ee94217 100644 --- a/src/render.rs +++ b/src/render.rs @@ -1,5 +1,5 @@ use ansi_term::Color::RGB; -use crate::geom::{Vec2, Vec3}; +use crate::geom::{Meters, Vec2, Vec3}; use crate::{flt::{Flt, Real}, Material as _}; use crate::mat; use crate::sim::{Cell, GenericSim}; @@ -92,7 +92,7 @@ impl<'a> RenderSteps<'a> { let y_prop = y_px as Flt / self.im.height() as Flt; let y_m = y_prop * (self.sim.height() as Flt * self.sim.feature_size()); let z_m = self.z as Flt * self.sim.feature_size(); - self.sim.sample(Vec3::new(x_m, y_m, z_m)) + self.sim.sample(Meters(Vec3::new(x_m, y_m, z_m))) } ////////////// Ex/Ey/Bz configuration //////////// diff --git a/src/sim.rs b/src/sim.rs index 804d5f4..64d2b9f 100644 --- a/src/sim.rs +++ b/src/sim.rs @@ -1,5 +1,5 @@ use crate::{flt::{Flt, Real}, consts}; -use crate::geom::{Coord, Vec3}; +use crate::geom::{Coord, Index, Meters, Vec3, Vec3u}; use crate::mat::{self, GenericMaterial, Material}; use log::trace; @@ -9,38 +9,61 @@ use std::convert::From; use std::iter::Sum; pub trait GenericSim { - fn sample(&self, pos_meters: Vec3) -> Cell; - /// TODO: DEPRECATED - fn get(&self, at: Coord) -> Cell { - let at_vec: Vec3 = at.into(); - self.sample(at_vec * self.feature_size()) - } - fn width(&self) -> u32; - fn height(&self) -> u32; - fn depth(&self) -> u32; + fn sample(&self, pos: Meters) -> Cell; + fn impulse_e_meters(&mut self, pos: Meters, amount: Vec3); + fn size(&self) -> Index; fn feature_size(&self) -> Flt; - fn time(&self) -> Flt; + fn timestep(&self) -> Flt; fn step_no(&self) -> u64; + + fn width(&self) -> u32 { + self.size().x() + } + fn height(&self) -> u32 { + self.size().y() + } + fn depth(&self) -> u32 { + self.size().z() + } + fn time(&self) -> Flt { + self.timestep() * self.step_no() as Flt + } } impl<'a> dyn GenericSim + 'a { + pub fn get(&self, at: C) -> Cell { + self.sample(at.to_meters(self.feature_size())) + } /// Apply `F` to each Cell, and sum the results. pub fn map_sum) -> R + , R: Sum>(&self, f: F) -> R { - self.map_sum_enumerated(|_at, cell| f(cell)) + self.map_sum_enumerated(|_at: Index, cell| f(cell)) } - pub fn map_sum_enumerated) -> R, R: Sum>(&self, f: F) -> R { + pub fn map_sum_enumerated) -> R, R: Sum>(&self, f: F) -> R { let (w, h, d) = (self.width(), self.height(), self.depth()); (0..d).map(|z| (0..h).map(|y| (0..w).map(|x| { - let at = Coord::new(x, y, z); - f(at, &self.get(at)) + let at = Index(Vec3u::new(x, y, z)); + f(C::from_index(at, self.feature_size()), &self.get(at)) }).sum()).sum()).sum() } - /// TODO: DEPRECATED - pub fn current(&self, c: Coord) -> Vec3 { + pub fn current(&self, c: C) -> Vec3 { self.get(c).current_density() * (self.feature_size() * self.feature_size()) } + + pub fn impulse_e(&mut self, pos: C, amt: Vec3) { + self.impulse_e_meters(pos.to_meters(self.feature_size()), amt) + } + + pub fn impulse_ex(&mut self, c: C, ex: Flt) { + self.impulse_e(c, Vec3::new(ex, 0.0, 0.0)); + } + pub fn impulse_ey(&mut self, c: C, ey: Flt) { + self.impulse_e(c, Vec3::new(0.0, ey, 0.0)); + } + pub fn impulse_ez(&mut self, c: C, ez: Flt) { + self.impulse_e(c, Vec3::new(0.0, 0.0, ez)); + } } #[derive(Default)] @@ -52,7 +75,7 @@ pub struct SimState { } impl SimState { - pub fn new(size: Coord, feature_size: Flt) -> Self { + pub fn new(size: Index, feature_size: Flt) -> Self { Self { cells: Array3::default((size.z() as _, size.y() as _, size.x() as _)), scratch: Array3::default((size.z() as _, size.y() as _, size.x() as _)), @@ -100,23 +123,10 @@ impl SimState { } } -impl SimState { - /// Apply `F` to each Cell, and sum the results. - pub fn map_sum) -> R + Sync + Send, R: Sum + Send>(&self, f: F) -> R { - self.cells.view().into_par_iter().map(f).sum() - } - - pub fn map_sum_enumerated) -> R + Sync + Send, R: Sum + Send>(&self, f: F) -> R { - Zip::from(ndarray::indices_of(&self.cells)).and(&self.cells).into_par_iter().map(|((z, y, x), c)| { - f(Coord::new(x as _, y as _, z as _), c) - }).sum() - } -} - impl GenericSim for SimState { - fn sample(&self, pos_meters: Vec3) -> Cell { + fn sample(&self, pos: Meters) -> Cell { // TODO: smarter sampling than nearest neighbor? - let pos_sim = pos_meters / self.feature_size(); + let pos_sim = pos.to_index(self.feature_size()); let idx = [pos_sim.z() as usize, pos_sim.y() as _, pos_sim.x() as _]; match self.cells.get(idx) { Some(cell) => Cell { @@ -130,20 +140,22 @@ impl GenericSim for SimState { } } - fn width(&self) -> u32 { - self.cells.shape()[2] as _ + fn impulse_e_meters(&mut self, pos: Meters, amount: Vec3) { + self.get_mut(pos).state.e += amount; } - fn height(&self) -> u32 { - self.cells.shape()[1] as _ - } - fn depth(&self) -> u32 { - self.cells.shape()[0] as _ + + fn size(&self) -> Index { + Index(Vec3u::new( + self.cells.shape()[2] as _, + self.cells.shape()[1] as _, + self.cells.shape()[0] as _, + )) } fn feature_size(&self) -> Flt { self.feature_size.into() } - fn time(&self) -> Flt { - self.timestep() * self.step_no as Flt + fn timestep(&self) -> Flt { + (self.feature_size / consts::real::C()).into() } fn step_no(&self) -> u64 { self.step_no @@ -151,31 +163,19 @@ impl GenericSim for SimState { } impl SimState { - pub fn impulse_bz(&mut self, c: Coord, bz: Flt) { + pub fn impulse_bz(&mut self, c: C, bz: Flt) { self.get_mut(c).impulse_bz(bz); } } impl SimState { - pub fn impulse_ex(&mut self, c: Coord, ex: Flt) { - self.get_mut(c).state.e += Vec3::new(ex, 0.0, 0.0); + pub fn get(&self, c: C) -> &Cell { + let at = c.to_index(self.feature_size.into()); + &self.cells[at.row_major_idx()] } - pub fn impulse_ey(&mut self, c: Coord, ey: Flt) { - self.get_mut(c).state.e += Vec3::new(0.0, ey, 0.0); - } - pub fn impulse_ez(&mut self, c: Coord, ez: Flt) { - self.get_mut(c).state.e += Vec3::new(0.0, 0.0, ez); - } - - pub fn get(&self, c: Coord) -> &Cell { - &self.cells[c.row_major_idx()] - } - pub fn get_mut(&mut self, c: Coord) -> &mut Cell { - &mut self.cells[c.row_major_idx()] - } - - pub fn timestep(&self) -> Flt { - (self.feature_size / consts::real::C()).into() + pub fn get_mut(&mut self, c: C) -> &mut Cell { + let at = c.to_index(self.feature_size.into()); + &mut self.cells[at.row_major_idx()] } }