Implement partial derivatives explicitly, for better precision

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;