diff --git a/src/bin/flat.rs b/src/bin/flat.rs index 692e708..0f55015 100644 --- a/src/bin/flat.rs +++ b/src/bin/flat.rs @@ -3,7 +3,7 @@ use std::f32::consts::PI; use flo_draw::*; use flo_canvas::*; use glam::*; -use riemann::{Decomp2, Metric, trace_iter}; +use riemann::{Tens2, Decomp2, Metric, trace_iter}; use shape::Shape; use Subspace::{Boundary, Inner, Outer}; @@ -591,6 +591,49 @@ impl Metric for Rect { diag: vec2(1.0, sy.lerp(1.0, sx)), } } + + fn part_derivs_at(&self, pos: Vec2) -> Tens2 { + let x = pos.x.abs(); + let y = pos.y.abs(); + let sx = ((x - self.inner_radius) / (self.outer_radius - self.inner_radius)).clamp(0.0, 1.0); + let sy = if y <= self.external_halflength { 1.0 / self.root(y) } else { 1.0 }; + let s = sy.lerp(1.0, sx); + let dsx_dx = if x > self.inner_radius && x < self.outer_radius { pos.x.signum() / (self.outer_radius - self.inner_radius) } else { 0.0 }; + let dsy_dy = if y <= self.external_halflength { -2.0 * (1.0 - self.γ()) / self.internal_halflength * sy.powi(3) * pos.y.signum() } else { 0.0 }; + let ds2_dx = 2.0 * s * (1.0 - sy) * dsx_dx; + let ds2_dy = 2.0 * s * (1.0 - sx) * dsy_dy; + [ + Mat2::from_cols_array(&[0.0, 0.0, 0.0, ds2_dx]), + Mat2::from_cols_array(&[0.0, 0.0, 0.0, ds2_dy]), + ] + } +} + +#[test] +fn test_rect_metric_derivs() { + struct Approx(Rect); + impl Metric for Approx { + fn sqrt_at(&self, pos: Vec2) -> Decomp2 { self.0.sqrt_at(pos) } + } + let testee = Rect { + inner_radius: 30.0, + outer_radius: 50.0, + internal_halflength: 100.0, + external_halflength: 300.0, + }; + let approx = Approx(testee); + let epsilon = 1.0e-3; + let margin = 1.0 / 16.0; + let mul = 1.0 + margin; + for x in itertools_num::linspace(-mul * testee.outer_radius, mul * testee.outer_radius, 100) { + for y in itertools_num::linspace(-mul * testee.external_halflength, mul * testee.external_halflength, 100) { + let pos = vec2(x, y); + let computed = testee.part_derivs_at(pos); + let reference = approx.part_derivs_at(pos); + let eq = (0..2).all(|coord| computed[coord].abs_diff_eq(reference[coord], epsilon)); + assert!(eq, "At {}:\nactual: {:?}\nexpected: {:?}\n", pos, computed, reference); + } + } } mod riemann { @@ -623,7 +666,7 @@ mod riemann { } } - type Tens2 = [Mat2; 2]; + pub type Tens2 = [Mat2; 2]; pub trait Metric { fn sqrt_at(&self, pos: Vec2) -> Decomp2;