diff --git a/examples/boundary_conditions.rs b/examples/boundary_conditions.rs index 9c463fa..a0609fc 100644 --- a/examples/boundary_conditions.rs +++ b/examples/boundary_conditions.rs @@ -3,9 +3,9 @@ use coremem::{Driver, mat}; fn main() { coremem::init_logging(); let width = 201; - let boundary = 20; + let boundary = 50; let mut driver = Driver::new(width, width, 1e-3 /* feature size */); - driver.add_y4m_renderer(&*format!("boundary_conditions_upml-{}.y4m", boundary)); + driver.add_y4m_renderer(&*format!("boundary_conditions_upml_round-{}.y4m", boundary)); // driver.add_term_renderer(); //driver.add_boundary(boundary, 0.1); diff --git a/src/driver.rs b/src/driver.rs index 58e1132..b369a5e 100644 --- a/src/driver.rs +++ b/src/driver.rs @@ -76,22 +76,32 @@ impl Driver { pub fn add_upml_boundary(&mut self, thickness: u32) { // Based on explanation here (slide 63): https://empossible.net/wp-content/uploads/2020/01/Lecture-The-Perfectly-Matched-Layer.pdf - for inset in 0..thickness { - let depth = thickness - inset; - let conductivity = 0.5 * consts::EPS0 / self.state.timestep() * (depth as Flt/thickness as Flt).powf(3.0); - for x in inset..self.state.width() - inset { - let conductor = mat::Static::anisotropic_conductor(Vec3::new(0.0, conductivity, 0.0)); - // left - *self.state.get_mut(x, inset).mat_mut() = conductor.clone().into(); - // right - *self.state.get_mut(x, self.state.height() - 1 - inset).mat_mut() = conductor.into(); - } - for y in inset..self.state.height() - inset { - let conductor = mat::Static::anisotropic_conductor(Vec3::new(conductivity, 0.0, 0.0)); - // top - *self.state.get_mut(inset, y).mat_mut() = conductor.clone().into(); - // bottom - *self.state.get_mut(self.state.width() - 1 - inset, y).mat_mut() = conductor.into(); + let h = self.state.height(); + let w = self.state.width(); + for y in 0..h { + for x in 0..w { + let depth_x = if x < thickness { + thickness - x + } else if x >= w - thickness { + w - x + } else { + 0 + }; + let depth_y = if y < thickness { + thickness - y + } else if y >= h - thickness { + h - y + } else { + 0 + }; + + if depth_x > 0 || depth_y > 0 { + let scale = 0.5 * consts::EPS0 / self.state.timestep(); + let cond_x = scale * (depth_x as Flt/thickness as Flt).powf(3.0); + let cond_y = scale * (depth_y as Flt/thickness as Flt).powf(3.0); + let conductor = mat::Static::anisotropic_conductor(Vec3::new(cond_x, cond_y, 0.0)); + *self.state.get_mut(x, y).mat_mut() = conductor.into(); + } } } }