Implement a BH-curve material (example material is 3R1)

NB: I think there are some errors needing to be worked out, but at least
it doesn't crash!

examples/em_reflection.rs

@@ -29,7 +29,7 @@ fn main() {
     loop {
         step += 1;
         let imp = if step < 50 { 
-            2.5 * ((step as f64)*0.04*std::f64::consts::PI).sin()
+            250000.0 * ((step as f64)*0.04*std::f64::consts::PI).sin()
         } else {
             0.0
         };

examples/ferromagnet.rs

@@ -13,11 +13,23 @@ fn main() {
     }
     for y in 40..60 {
         for x in 62..70 {
-            *state.get_mut(x, y).mat_mut() = mat::ComsolFerromagnet {
+            *state.get_mut(x, y).mat_mut() = mat::PiecewiseLinearFerromagnet::from_bh(&[
-                xi: 150.0,
+                (  35.0, 0.0),
-                ms: 2.0,
+                (  50.0, 0.250),
-                ..mat::ComsolFerromagnet::default()
+                ( 100.0, 0.325),
-            }.into();
+                ( 200.0, 0.350),
+                (1000.0, 0.390),
+                // Falling
+                ( 200.0, 0.360),
+                ( 100.0, 0.345),
+                (  50.0, 0.340),
+                (   0.0, 0.325),
+            ]).into();
+            //*state.get_mut(x, y).mat_mut() = mat::ComsolFerromagnet {
+            //    xi: 150.0,
+            //    ms: 2.0,
+            //    ..mat::ComsolFerromagnet::default()
+            //}.into();
         }
     }
 
@@ -25,7 +37,7 @@ fn main() {
     loop {
         step += 1;
         let imp = if step < 50 { 
-            25.0 * ((step as f64)*0.02*std::f64::consts::PI).sin()
+            250000.0 * ((step as f64)*0.02*std::f64::consts::PI).sin()
         } else {
             0.0
         };

src/geom.rs

@@ -0,0 +1,174 @@
+use decorum::R64;
+use std::ops::{Add, Neg};
+
+#[derive(Copy, Clone, Debug, Default)]
+pub struct Point {
+    pub x: f64,
+    pub y: f64,
+}
+
+impl Add for Point {
+    type Output = Point;
+    fn add(self, other: Point) -> Point {
+        Point {
+            x: self.x + other.x,
+            y: self.y + other.y,
+        }
+    }
+}
+
+impl Neg for Point {
+    type Output = Point;
+    fn neg(self) -> Point {
+        Point {
+            x: -self.x,
+            y: -self.y,
+        }
+    }
+}
+
+impl Point {
+    pub fn new(x: f64, y: f64) -> Self {
+        Self {
+            x,
+            y,
+        }
+    }
+}
+
+#[derive(Copy, Clone, Debug, Default)]
+pub struct Line {
+    from: Point,
+    to: Point,
+}
+
+impl Add for Line {
+    type Output = Line;
+    fn add(self, other: Line) -> Line {
+        Line {
+            from: self.from + Point::new(0.0, other.y(self.from.x)),
+            to: self.to + Point::new(0.0, other.y(self.to.x)),
+        }
+    }
+}
+
+impl Line {
+    pub fn new(from: Point, to: Point) -> Self {
+        Self {
+            from,
+            to,
+        }
+    }
+
+    pub fn x(&self, y: f64) -> f64 {
+        self.shifted(Point::new(0.0, -y)).x_intercept()
+    }
+    pub fn y(&self, x: f64) -> f64 {
+        self.from.y + (x - self.from.x) * self.slope()
+    }
+    pub fn at_x(&self, x: f64) -> Point {
+        Point::new(x, self.y(x))
+    }
+
+    pub fn x_intercept(&self) -> f64 {
+        let t = (0.0 - self.from.y) / (self.to.y - self.from.y);
+        self.from.x + t * (self.to.x - self.from.x)
+    }
+
+    pub fn is_ascending(&self) -> bool {
+        self.to.y > self.from.y
+    }
+
+    pub fn is_ascending_x(&self) -> bool {
+        self.to.x > self.from.x
+    }
+
+    pub fn contains_x(&self, x: f64) -> bool {
+        if self.is_ascending_x() {
+            x >= self.from.x && x < self.to.x
+        } else {
+            x >= self.to.x && x < self.from.x
+        }
+    }
+
+    pub fn contains_y(&self, y: f64) -> bool {
+        if self.is_ascending_x() {
+            y >= self.from.y && y < self.to.y
+        } else {
+            y >= self.to.y && y < self.from.y
+        }
+    }
+
+    pub fn clamp_x(&self, x: f64) -> f64 {
+        match x {
+            v if v < self.from.x => self.from.x,
+            v if v > self.to.x => self.to.x,
+            v => v,
+        }
+    }
+
+    pub fn clamp_by_x(&self, x: f64) -> Point {
+        self.at_x(self.clamp_x(x))
+    }
+
+    pub fn slope(&self) -> f64 {
+        (self.to.y - self.from.y) / (self.to.x - self.from.x)
+    }
+
+    /// Move the coordinate along the line toward some desired `y(x)` value.
+    /// The returned value could be OOB: check `self.contains_x` to find out.
+    pub fn move_toward_y_unclamped(&self, start_x: f64, target_y: f64) -> f64 {
+        let start_y = self.y(start_x);
+        let delta_x = (target_y - start_y) / self.slope();
+        start_x + delta_x
+    }
+
+    pub fn move_toward_y(&self, start_x: f64, target_y: f64) -> f64 {
+        let x = self.move_toward_y(start_x, target_y);
+        self.clamp_x(x)
+    }
+
+    pub fn shifted(&self, p: Point) -> Line {
+        Line {
+            from: self.from + p,
+            to: self.to + p
+        }
+    }
+}
+
+#[derive(Clone, Debug)]
+pub struct Polygon {
+    points: Vec<Point>,
+}
+
+impl Polygon {
+    pub fn new(points: Vec<Point>) -> Self {
+        Self {
+            points
+        }
+    }
+
+    pub fn segments<'a>(&'a self) -> impl Iterator<Item=Line> + 'a {
+        (0..self.points.len()).into_iter().map(move |i| {
+            let from = self.points[i];
+            let to = *self.points.get(i+1).unwrap_or_else(|| &self.points[0]);
+            Line { from, to }
+        })
+    }
+
+    pub fn min_x(&self) -> f64 {
+        self.points.iter().map(|p| R64::from_inner(p.x)).min().unwrap().into()
+    }
+
+    pub fn max_x(&self) -> f64 {
+        self.points.iter().map(|p| R64::from_inner(p.x)).max().unwrap().into()
+    }
+
+    pub fn min_y(&self) -> f64 {
+        self.points.iter().map(|p| R64::from_inner(p.y)).min().unwrap().into()
+    }
+
+    pub fn max_y(&self) -> f64 {
+        self.points.iter().map(|p| R64::from_inner(p.y)).max().unwrap().into()
+    }
+}

src/lib.rs

@@ -7,6 +7,7 @@
 
 use decorum::R64;
 
+pub mod geom;
 pub mod mat;
 pub mod render;
 pub mod sim;

src/mat.rs

@@ -1,6 +1,6 @@
 use crate::{CellState, consts};
+use crate::geom::{Line, Point, Polygon};
 use enum_dispatch::enum_dispatch;
-use piecewise_linear::PiecewiseLinearFunction;
 
 #[enum_dispatch]
 pub trait Material {
@@ -80,67 +80,108 @@ impl Material for ComsolFerromagnet {
     }
 }
 
-pub struct BHCurve {
+/// M as a function of H
-    /// describes the curve as B rises from 0.
+#[derive(Clone)]
-    rising: PiecewiseLinearFunction<f64>,
+struct MHCurve {
-    /// describes the curve as B falls towards 0.
+    geom: Polygon,
-    falling: PiecewiseLinearFunction<f64>,
 }
 
-impl BHCurve {
+impl MHCurve {
-    fn b_bounds_for_h(&self, h: f64) -> (f64, f64) {
+    /// Construct a M(H) curve from a sweep from M = 0 to Ms and back down to M = 0.
-        let (rise_start, rise_end) = self.rising.domain();
+    /// The curve below M = 0 is derived by symmetry.
-        let bottom = match h {
+    pub fn new(points: &[Point]) -> Self {
-            x if (..-rise_end).contains(&x) => {
+        let full_pts: Vec<Point> = points.iter().cloned().chain(points.iter().cloned().map(|p| -p)).collect();
-                // negative saturation
+
-                -self.falling.y_at_x(-rise_end).unwrap() + consts::MU0 * (h + rise_end)
+        Self {
-            },
+            geom: Polygon::new(full_pts)
-            x if (-rise_end..rise_start).contains(&x) => {
+        }
-                // negative
+    }
-                -self.falling.y_at_x(-h).unwrap()
+    /// Moves (h, m) towards some location in the MH curve where H + M = target_hm.
+    /// Returns `Ok((h, m))` if complete; `Err((h, m))` if there's more work to be done (call it
+    /// again).
+    fn step_toward(&self, h: f64, m: f64, target_hm: f64) -> Result<Point, Point> {
+        let is_ascending = target_hm > h + m;
+        if (is_ascending && m == self.geom.max_y()) || (!is_ascending && m == -self.geom.min_y()) {
+            // Fully saturated. m is fixed, while h moves freely
+            return Ok(Point::new(target_hm - m, m));
+        }
+        // Locate the segment which would contain the current point
+        let mut segments = self.geom.segments();
+        let active_segment = loop {
+            let line = segments.next().unwrap_or_else(|| {
+                panic!("failed to find segment for h:{}, m:{}, {:?}", h, m, self.geom.segments().collect::<Vec<_>>());
+            });
+            if line.contains_y(m) && line.is_ascending() == is_ascending {
+                if line.contains_x(h) {
+                    break line;
+                } else {
+                    // need to move the point toward this line
+                    let h_intercept = line.x(m);
+                    break Line::new(Point::new(h, m), Point::new(h_intercept, m));
                 }
-            x if (rise_start..rise_end).contains(&x) => {
-                // positive
-                self.rising.y_at_x(h).unwrap()
-            },
-            x if (rise_end..).contains(&x) => {
-                // positive saturation
-                self.rising.y_at_x(rise_end).unwrap() + consts::MU0 * (h - rise_end)
             }
-            _ => unreachable!("invalid h: {}", h),
         };
-        let top = match h {
+
-            x if (..-rise_end).contains(&x) => {
+        // Find some m(h) on the active_segment such that sum(h) = h + m(h) = target_hm
-                // negative saturation
+        let sum_h = active_segment + Line::new(Point::new(0.0, 0.0), Point::new(1.0, 1.0));
-                -self.rising.y_at_x(-rise_end).unwrap() + consts::MU0 * (h + rise_end)
+        let new_h = sum_h.move_toward_y_unclamped(h, target_hm);
+
+        if sum_h.contains_x(new_h) {
+            // the segment contains a point with the target H+M
+            Ok(active_segment.at_x(new_h))
+        } else {
+            // the segment doesn't contain the desired point: clamp and try the next segment
+            Err(active_segment.clamp_by_x(new_h))
+        }
+    }
+    fn move_to(&self, mut h: f64, mut m: f64, target_hm: f64) -> (f64, f64) {
+        loop {
+            match self.step_toward(h, m, target_hm) {
+                Ok(p) => break (p.x, p.y),
+                Err(Point{ x: better_h, y: better_m}) => {
+                    h = better_h;
+                    m = better_m;
                 },
-            x if (-rise_end..-rise_start).contains(&x) => {
+            }
-                // negative
+        }
-                -self.rising.y_at_x(-h).unwrap()
-            },
-            x if (-rise_start..rise_end).contains(&x) => {
-                // positive
-                self.falling.y_at_x(h).unwrap()
-            },
-            x if (rise_end..).contains(&x) => {
-                // positive saturation
-                self.falling.y_at_x(rise_end).unwrap() + consts::MU0 * (h - rise_end)
-            },
-            _ => unreachable!("invalid h: {}", h),
-        };
-        (bottom, top)
     }
 }
 
-pub struct BHFerromagnet {
+#[derive(Clone)]
+pub struct PiecewiseLinearFerromagnet {
     /// Instantaneous magnetization in the z direction.
     pub mz: f64,
-    bh_curve: BHCurve,
+    mh_curve: MHCurve,
 }
 
-impl Material for BHFerromagnet {
+impl PiecewiseLinearFerromagnet {
+    /// Construct a B(H) curve from a sweep from H = 0 to Hs and back down to H = 0.
+    /// The curve below B = 0 is derived by symmetry.
+    /// Points are (H, B).
+    pub fn from_bh(points: &[(f64, f64)]) -> Self {
+        let mh_points: Vec<Point> = points.iter().cloned().map(|(h, b)| {
+            Point::new(h, b / consts::MU0 - h)
+        }).collect();
+
+        Self {
+            mz: 0.0,
+            mh_curve: MHCurve::new(&*mh_points),
+        }
+    }
+}
+
+impl Material for PiecewiseLinearFerromagnet {
+    fn mz(&self) -> f64 {
+        self.mz
+    }
     fn step_h(&mut self, context: &CellState, delta_bz: f64) -> f64 {
-        unimplemented!()
+        let (h, m) = (context.hz(), self.mz);
+        let target_hm = h + m + delta_bz/consts::MU0;
+
+        let (h, m) = self.mh_curve.move_to(h, m, target_hm);
+
+        self.mz = m;
+        h
     }
 }
 
@@ -148,7 +189,8 @@ impl Material for BHFerromagnet {
 #[derive(Clone)]
 pub enum GenericMaterial {
     Conductor(Conductor),
-    ComsolFerromagnet(ComsolFerromagnet),
+    // ComsolFerromagnet(ComsolFerromagnet),
+    PiecewiseLinearFerromagnet(PiecewiseLinearFerromagnet)
 }
 
 impl Default for GenericMaterial {

src/render.rs

@@ -50,14 +50,14 @@ impl ColorTermRenderer {
                 let cell = state.get(x, y);
                 //let r = norm_color(cell.bz() * consts::C);
                 //let r = 0;
-                let r = norm_color(cell.mat().mz());
+                let r = norm_color(cell.mat().mz()*1.0e-2);
                 let b = (55.0*cell.mat().conductivity()).min(255.0) as u8;
                 //let b = 0;
                 //let b = norm_color(cell.ey());
                 //let g = 0;
                 //let g = norm_color(cell.ex());
                 //let g = norm_color(curl(cell.ex(), cell.ey()));
-                let g = norm_color((cell.bz() * 3.0e8).into());
+                let g = norm_color((cell.bz() * 1.0e4).into());
                 //let g = norm_color(cell.ey().into());
                 write!(&mut buf, "{}", RGB(r, g, b).paint(square));
             }