diff --git a/src/mat.rs b/src/mat.rs
index 9218538..41abff7 100644
--- a/src/mat.rs
+++ b/src/mat.rs
@@ -1,5 +1,7 @@
-use crate::{CellState, consts, flt::Flt};
+use crate::{CellState, consts};
+use crate::flt::Flt;
 use crate::geom::{Line, Point, Polygon};
+use crate::vec3::Vec3;
 use enum_dispatch::enum_dispatch;
 use std::cmp::Ordering;
 
@@ -10,21 +12,21 @@ pub trait Material {
     fn conductivity(&self) -> Flt {
         0.0
     }
-    /// Return the magnetization in the z direction (M_z).
-    fn mz(&self) -> Flt {
-        0.0
+    /// Return the magnetization.
+    fn m(&self) -> Vec3 {
+        Vec3::zero()
     }
-    /// Return the new h_z, and optionally change any internal state (e.g. magnetization).
-    fn step_h(&mut self, context: &CellState, delta_bz: Flt) -> Flt {
+    /// Return the new h, and optionally change any internal state (e.g. magnetization).
+    fn step_h(&mut self, context: &CellState, delta_b: Vec3) -> Vec3 {
         // B = mu0*(H+M), and in a default material M=0.
-        context.hz() + delta_bz / consts::MU0
+        context.h() + delta_b / consts::real::MU0()
     }
 }
 
 #[derive(Clone, Default)]
 pub struct Static {
     pub conductivity: Flt,
-    pub mz: Flt,
+    pub m: Vec3,
 }
 
 impl Static {
@@ -40,57 +42,8 @@ impl Material for Static {
     fn conductivity(&self) -> Flt {
         self.conductivity
     }
-    fn mz(&self) -> Flt {
-        self.mz
-    }
-}
-
-#[derive(Clone, Default)]
-pub struct ComsolFerromagnet {
-    /// Instantaneous magnetization in the z direction.
-    pub mz: Flt,
-    // Ferromagnetic parameters:
-    /// Magnetic saturation (symmetric, so saturates to either +Ms or -Ms).
-    pub ms: Flt,
-    /// Rate at which M changes w.r.t. H.
-    pub xi: Flt,
-}
-
-impl Material for ComsolFerromagnet {
-    fn mz(&self) -> Flt {
-        self.mz
-    }
-    fn step_h(&mut self, context: &CellState, delta_bz: Flt) -> Flt {
-        // COMSOL ferromagnetic model: https://www.comsol.com/blogs/accessing-external-material-models-for-magnetic-simulations/
-        //   delta(M) = xi * delta*(H)            # delta(M) = M(t+1) - M(t)
-        //   new_M = clamp(M+delta(M), -Ms, Ms))  # i.e. clamp to saturation
-        //   Where delta*(H) is H(t+1) - H(t), but clamping H(t) to <= 0 if positively saturated
-        //   and clamping H(t) to >= 0 if negatively saturated
-        let expected_delta_h = delta_bz / consts::MU0;
-        let expected_new_h = context.hz() + expected_delta_h;
-        // let unsaturated_new_m = self.mz + self.xi * unsaturated_delta_h;
-        let unsaturated_new_m = if self.mz == self.ms {
-            // positively saturated
-            let delta_h = expected_new_h.min(0.0) - context.hz().min(0.0);
-            // assert!(self.ms == 0.0, "pos sat");
-            self.mz + self.xi * delta_h
-        } else if self.mz == -self.ms {
-            // negatively saturated
-            let delta_h = expected_new_h.max(0.0) - context.hz().max(0.0);
-            //assert!(self.ms == 0.0, "neg sat");
-            self.mz + self.xi * delta_h
-        } else {
-            // not saturated
-            self.mz + self.xi * expected_delta_h
-        };
-        let new_m = unsaturated_new_m.min(self.ms).max(-self.ms);
-
-        let delta_m = new_m - self.mz;
-        self.mz = new_m;
-        // B = mu0*(H + M)
-        // \Delta B = mu0*(\Delta H + \Delta M)
-        let delta_h = delta_bz / consts::MU0 - delta_m;
-        context.hz() + delta_h
+    fn m(&self) -> Vec3 {
+        self.m
     }
 }
 
@@ -169,9 +122,11 @@ impl MHCurve {
 #[derive(Clone)]
 pub struct PiecewiseLinearFerromagnet {
     /// Instantaneous magnetization in the z direction.
-    pub mz: Flt,
+    pub m: Vec3,
     pub conductivity: Flt,
-    mh_curve: MHCurve,
+    mh_curve_x: MHCurve,
+    mh_curve_y: MHCurve,
+    mh_curve_z: MHCurve,
 }
 
 impl PiecewiseLinearFerromagnet {
@@ -184,28 +139,33 @@ impl PiecewiseLinearFerromagnet {
         }).collect();
 
         Self {
-            mz: 0.0,
+            m: Vec3::default(),
             conductivity: 0.0,
-            mh_curve: MHCurve::new(&*mh_points),
+            mh_curve_x: MHCurve::new(&*mh_points),
+            mh_curve_y: MHCurve::new(&*mh_points),
+            mh_curve_z: MHCurve::new(&*mh_points),
         }
     }
 }
 
 impl Material for PiecewiseLinearFerromagnet {
-    fn mz(&self) -> Flt {
-        self.mz
+    fn m(&self) -> Vec3 {
+        self.m
     }
     fn conductivity(&self) -> Flt {
         self.conductivity
     }
-    fn step_h(&mut self, context: &CellState, delta_bz: Flt) -> Flt {
-        let (h, m) = (context.hz(), self.mz);
-        let target_hm = h + m + delta_bz/consts::MU0;
+    fn step_h(&mut self, context: &CellState, delta_b: Vec3) -> Vec3 {
+        let (h, m) = (context.h(), self.m);
+        let target_hm = h + m + delta_b/consts::real::MU0();
 
-        let (h, m) = self.mh_curve.move_to(h, m, target_hm);
+        // TODO: this is probably not the best way to generalize a BH curve into 3d.
+        let (hx, mx) = self.mh_curve_x.move_to(h.x(), m.x(), target_hm.x());
+        let (hy, my) = self.mh_curve_x.move_to(h.y(), m.y(), target_hm.y());
+        let (hz, mz) = self.mh_curve_x.move_to(h.z(), m.z(), target_hm.z());
 
-        self.mz = m;
-        h
+        self.m = Vec3::new(mx, my, mz);
+        Vec3::new(hx, hy, hz)
     }
 }
 
diff --git a/src/meas.rs b/src/meas.rs
index 9b7346c..75cb206 100644
--- a/src/meas.rs
+++ b/src/meas.rs
@@ -49,8 +49,8 @@ pub struct Magnetization(pub u32, pub u32);
 
 impl AbstractMeasurement for Magnetization {
     fn eval(&self, state: &SimSnapshot) -> String {
-        let mz = state.get(self.0, self.1).mat().mz();
-        format!("Mz({}, {}): {:.2e}", self.0, self.1, mz)
+        let m = state.get(self.0, self.1).mat().m();
+        format!("M({}, {}): ({:.2e}, {:.2e}, {:.2e})", self.0, self.1, m.x(), m.y(), m.z())
     }
 }
 
@@ -59,8 +59,8 @@ pub struct MagneticFlux(pub u32, pub u32);
 
 impl AbstractMeasurement for MagneticFlux {
     fn eval(&self, state: &SimSnapshot) -> String {
-        let bz = state.get(self.0, self.1).bz();
-        format!("Bz({}, {}): {:.2e}", self.0, self.1, bz)
+        let b = state.get(self.0, self.1).b();
+        format!("Bz({}, {}): ({:.2e}, {:.2e}, {:.2e})", self.0, self.1, b.x(), b.y(), b.z())
     }
 }
 
@@ -69,8 +69,8 @@ pub struct MagneticStrength(pub u32, pub u32);
 
 impl AbstractMeasurement for MagneticStrength {
     fn eval(&self, state: &SimSnapshot) -> String {
-        let bz = state.get(self.0, self.1).hz();
-        format!("Hz({}, {}): {:.2e}", self.0, self.1, bz)
+        let h = state.get(self.0, self.1).h();
+        format!("Hz({}, {}): ({:.2e}, {:.2e}, {:.2e})", self.0, self.1, h.x(), h.y(), h.z())
     }
 }
 
@@ -82,9 +82,10 @@ impl AbstractMeasurement for Energy {
         #[allow(non_snake_case)]
         let dV = f*f*f;
         let e = state.map_sum(|cell| {
-            0.5 * cell.hz() * cell.bz() * dV
+            0.5 * cell.h().z() * cell.b().z() * dV
         });
-        format!("U: {:.2e}", e)
+        // TODO: update to 3 dimensions!
+        format!("U(TODO): {:.2e}", e)
     }
 }
 
diff --git a/src/render.rs b/src/render.rs
index 8e3d465..6516ef3 100644
--- a/src/render.rs
+++ b/src/render.rs
@@ -81,9 +81,9 @@ impl<'a> RenderSteps<'a> {
         for y in 0..h {
             for x in 0..w {
                 let cell = self.sim.get(x, y);
-                let r = scale_signed_to_u8(cell.mat().mz(), 100.0);
+                let r = scale_signed_to_u8(cell.mat().m().z(), 100.0);
                 let b = scale_unsigned_to_u8(cell.mat().conductivity(), 10.0);
-                let g = scale_signed_to_u8(cell.bz(), 1.0e-4);
+                let g = scale_signed_to_u8(cell.b().z(), 1.0e-4);
                 self.im.put_pixel(x, y, Rgb([r, g, b]));
             }
         }
diff --git a/src/sim.rs b/src/sim.rs
index 4a13790..5147ba8 100644
--- a/src/sim.rs
+++ b/src/sim.rs
@@ -75,14 +75,14 @@ impl<M: Material + Default + Send + Sync> SimState<M> {
                 let mut cell: Cell<mat::Static> = Default::default();
                 for sample in cells {
                     cell.state.e += sample.state.e();
-                    cell.state.hz += sample.hz();
+                    cell.state.h += sample.h();
                     cell.mat.conductivity += sample.mat.conductivity();
-                    cell.mat.mz += sample.mat.mz();
+                    cell.mat.m += sample.mat.m();
                 }
                 cell.state.e /= Real::from_inner(sample_area);
-                cell.state.hz /= sample_area;
+                cell.state.h /= Real::from_inner(sample_area);
                 cell.mat.conductivity /= sample_area;
-                cell.mat.mz /= sample_area;
+                cell.mat.m /= Real::from_inner(sample_area);
                 cell
             });
 
@@ -181,6 +181,9 @@ impl<M> Cell<M> {
     pub fn hz(&self) -> Flt {
         self.state.hz()
     }
+    pub fn h(&self) -> Vec3 {
+        self.state.h()
+    }
     pub fn mat(&self) -> &M {
         &self.mat
     }
@@ -190,8 +193,11 @@ impl<M> Cell<M> {
 }
 
 impl<M: Material> Cell<M> {
+    pub fn b(&self) -> Vec3 {
+        (self.h() + self.mat.m()) * consts::real::MU0()
+    }
     pub fn bz(&self) -> Flt {
-        consts::MU0 * (self.hz() + self.mat.mz())
+        self.b().z()
     }
 
     pub fn current(&self) -> Point {
@@ -200,7 +206,7 @@ impl<M: Material> Cell<M> {
     }
 
     fn impulse_bz(&mut self, delta_bz: Flt) {
-        self.state.hz = self.mat.step_h(&self.state, delta_bz).into();
+        self.state.h = self.mat.step_h(&self.state, Vec3::new(0.0, 0.0, delta_bz));
     }
 
     fn step_h(mut self, right: &Self, down: &Self, delta_t: Real, feature_size: Real) -> Self {
@@ -227,10 +233,10 @@ impl<M: Material> Cell<M> {
         let nabla_e = Vec3::from((delta_ez_y - delta_ey_z, delta_ex_z - delta_ez_x, delta_ey_x - delta_ex_y));
 
         let delta_b = -nabla_e * delta_t;
-        let hz = self.mat.step_h(&self.state, delta_b.z());
+        let h = self.mat.step_h(&self.state, delta_b);
         Cell {
             state: CellState {
-                hz: hz.into(),
+                h: h.into(),
                 ..self.state
             },
             mat: self.mat,
@@ -284,7 +290,7 @@ impl<M: Material> Cell<M> {
         Cell {
             state: CellState {
                 e: e_next,
-                hz: self.state.hz,
+                h: self.state.h,
             },
             mat: self.mat,
         }
@@ -294,7 +300,7 @@ impl<M: Material> Cell<M> {
 #[derive(Copy, Clone, Default)]
 pub struct CellState {
     e: Vec3,
-    hz: Real,
+    h: Vec3,
 }
 
 impl CellState {
@@ -305,7 +311,10 @@ impl CellState {
         self.e.y()
     }
     pub fn hz(&self) -> Flt {
-        self.hz.into()
+        self.h.z().into()
+    }
+    pub fn h(&self) -> Vec3 {
+        self.h
     }
     pub fn e(&self) -> Vec3 {
         self.e
diff --git a/src/vec3.rs b/src/vec3.rs
index 51bbc95..ad06979 100644
--- a/src/vec3.rs
+++ b/src/vec3.rs
@@ -17,6 +17,12 @@ impl Vec3 {
             z: z.into()
         }
     }
+    pub fn unit() -> Self {
+        Self::new(1.0, 1.0, 1.0)
+    }
+    pub fn zero() -> Self {
+        Self::new(0.0, 0.0, 0.0)
+    }
     pub fn x(&self) -> Flt {
         self.x.into()
     }