Add this new UpdateMethod enum and PmlState to allow for future pml implementation

I don't like the mutability peculiarities around it, but it should do
the job for now.

src/mat.rs

@@ -1,6 +1,7 @@
 use crate::{CellState, consts};
 use crate::flt::{Flt, Real};
 use crate::geom::{Line2d, Vec2, Vec3, Polygon2d};
+use crate::sim::{UpdateMethod, UpdateMethodMut};
 
 use lazy_static::lazy_static;
 use log::{debug, trace};
@@ -10,10 +11,9 @@ use std::cmp::Ordering;
 
 #[enum_dispatch]
 pub trait Material {
-    /// Return \sigma, the electrical conductivity.
-    /// For a vacuum, this is zero. For a perfect conductor, \inf.
-    fn conductivity(&self) -> Vec3 {
-        Default::default()
+    fn update_method_mut<'a>(&'a mut self) -> UpdateMethodMut<'a> {
+        // by default, behave as a vacuum
+        UpdateMethodMut::Maxwell { conductivity: Vec3::zero() }
     }
     /// Return the magnetization.
     fn m(&self) -> Vec3 {
@@ -27,12 +27,27 @@ pub trait Material {
 }
 
 pub trait MaterialExt {
+    fn update_method<'a>(&'a self) -> UpdateMethod<'a>;
+    fn conductivity(&self) -> Vec3;
     fn is_vacuum(&self) -> bool;
 }
 
 impl<M: Material> MaterialExt for M {
+    fn update_method<'a>(&'a self) -> UpdateMethod<'a> {
+        match unsafe { &mut *(self as *const M as *mut M) }.update_method_mut() {
+            UpdateMethodMut::Maxwell { conductivity } => UpdateMethod::Maxwell { conductivity},
+            UpdateMethodMut::Pml { pml, pseudo_conductivity } =>
+                UpdateMethod::Pml { pml, pseudo_conductivity },
+        }
+    }
+    fn conductivity(&self) -> Vec3 {
+        match self.update_method() {
+            UpdateMethod::Maxwell { conductivity } => conductivity,
+            UpdateMethod::Pml { pseudo_conductivity, .. } => pseudo_conductivity,
+        }
+    }
     fn is_vacuum(&self) -> bool {
-        self.conductivity() == Default::default()
+        self.conductivity() == Vec3::zero()
     }
 }
 
@@ -54,8 +69,8 @@ impl Static {
 }
 
 impl Material for Static {
-    fn conductivity(&self) -> Vec3 {
-        self.conductivity
+    fn update_method_mut<'a>(&'a mut self) -> UpdateMethodMut<'a> {
+        UpdateMethodMut::Maxwell { conductivity: self.conductivity }
     }
     fn m(&self) -> Vec3 {
         self.m
@@ -89,8 +104,8 @@ impl Conductor {
 }
 
 impl Material for Conductor {
-    fn conductivity(&self) -> Vec3 {
-        self.conductivity
+    fn update_method_mut<'a>(&'a mut self) -> UpdateMethodMut<'a> {
+        UpdateMethodMut::Maxwell { conductivity: self.conductivity }
     }
 }
 
@@ -102,9 +117,10 @@ pub trait PiecewiseLinearFerromagnet {
 }
 
 impl<T: PiecewiseLinearFerromagnet> Material for T {
-    fn conductivity(&self) -> Vec3 {
+    fn update_method_mut<'a>(&'a mut self) -> UpdateMethodMut<'a> {
         let c = T::conductivity();
-        Vec3::new(c, c, c)
+        let conductivity = Vec3::new(c, c, c);
+        UpdateMethodMut::Maxwell { conductivity }
     }
     fn m(&self) -> Vec3 {
         self.m()

src/sim.rs

@@ -1,6 +1,6 @@
 use crate::{flt::{Flt, Real}, consts};
 use crate::geom::{Coord, Cube, Index, InvertedRegion, Meters, Region, Vec3, Vec3u};
-use crate::mat::{self, GenericMaterial, Material};
+use crate::mat::{self, GenericMaterial, Material, MaterialExt as _};
 use crate::stim::AbstractStimulus;
 use dyn_clone::{self, DynClone};
 use log::trace;
@@ -13,6 +13,43 @@ use std::iter::Sum;
 
 pub type StaticSim = SimState<mat::Static>;
 
+#[derive(Default, Copy, Clone, PartialEq)]
+pub struct PmlState {
+    /// Auxiliary field used when stepping E
+    aux_e: Vec3,
+    /// Auxiliary field used when stepping B
+    aux_b: Vec3,
+}
+
+impl PmlState {
+    pub fn new() -> Self {
+        Self::default()
+    }
+}
+
+/// Used to instruct the Sim how to advance the state at a given location.
+///
+/// `UpdateMethod::Maxwell` means to apply Maxwell's equations as usual, with
+/// the provided parameters.
+///
+/// `UpdateMethod::Pml` means to apply a modified formulation of Maxwell's equations
+/// which are reflectionless for normally incident waves. This allows the material to behave
+/// as a "Perfectly Matched Layer" (PML), dissipating energy AS IF it was a conductor
+/// with the provided `pseudo_conductivity`.
+#[derive(PartialEq)]
+pub enum UpdateMethodMut<'a> {
+    /// `conductivity` is traditionally represented by $\sigma$
+    Maxwell { conductivity: Vec3 },
+    Pml { pml: &'a mut PmlState, pseudo_conductivity: Vec3 },
+}
+
+#[derive(PartialEq)]
+pub enum UpdateMethod<'a> {
+    /// `conductivity` is traditionally represented by $\sigma$
+    Maxwell { conductivity: Vec3 },
+    Pml { pml: &'a PmlState, pseudo_conductivity: Vec3 },
+}
+
 pub trait GenericSim: Send + Sync + DynClone {
     fn sample(&self, pos: Meters) -> Cell<mat::Static>;
     fn impulse_e_meters(&mut self, pos: Meters, amount: Vec3);
@@ -580,6 +617,7 @@ impl<M: Material> Cell<M> {
     }
 
     pub fn current_density(&self) -> Vec3 {
+        // TODO: does this make sense for Pml?
         let conductivity = self.mat.conductivity();
         self.e().elem_mul(conductivity)
     }
@@ -664,7 +702,6 @@ impl<M: Material> Cell<M> {
         let twice_delta_t = TWO() * delta_t;
         let half_delta_t = HALF() * delta_t;
 
-        let h_prev = self.h();
         let b_prev = self.b();
 
         let (delta_ey_x, delta_ez_x) = match (neighbors.right, neighbors.left) {
@@ -691,6 +728,8 @@ impl<M: Material> Cell<M> {
         // \nabla x E
         let nabla_e = Vec3::from((delta_ez_y - delta_ey_z, delta_ex_z - delta_ez_x, delta_ey_x - delta_ex_y));
 
+        match self.mat.update_method() {
+            UpdateMethod::Maxwell { .. } => {
                 // From: $\nabla x E = -dB/dt$
                 let db_dt = solve_step_diff_eq(
                     nabla_e,
@@ -715,6 +754,9 @@ impl<M: Material> Cell<M> {
                     mat: self.mat,
                 }
             }
+            UpdateMethod::Pml { .. } => unimplemented!(),
+        }
+    }
 
     /// delta_t = timestep covered by this function. i.e. it should be half the timestep of the simulation
     /// feature_size = how many units apart is the center of each adjacent cell on the grid.
@@ -751,9 +793,6 @@ impl<M: Material> Cell<M> {
         use consts::real::{EPS0, HALF, THIRD, TWO, ZERO};
         let twice_delta_t = TWO() * delta_t;
 
-        let sigma = self.mat.conductivity();
-
-
         let e_prev = self.e();
 
         let (delta_hy_x, delta_hz_x) = match (neighbors.left, neighbors.right) {
@@ -780,6 +819,9 @@ impl<M: Material> Cell<M> {
         // \nabla x H
         let nabla_h = Vec3::from((delta_hz_y - delta_hy_z, delta_hx_z - delta_hz_x, delta_hy_x - delta_hx_y));
 
+        match self.mat.update_method() {
+            UpdateMethod::Maxwell { conductivity } => {
+                let sigma = conductivity;
 
                 // From: $\nabla x H = \epsilon_0 dE/dt + \sigma E$
                 let de_dt = solve_step_diff_eq(
@@ -802,6 +844,9 @@ impl<M: Material> Cell<M> {
                     mat: self.mat,
                 }
             }
+            UpdateMethod::Pml { .. } => unimplemented!(),
+        }
+    }
 }
 
 #[derive(Copy, Clone, Debug, Default, PartialEq, Serialize, Deserialize)]