From 84068a5a1379ad5c2a723a0da613fd264e72ca37 Mon Sep 17 00:00:00 2001 From: numzero Date: Tue, 25 Jun 2024 12:16:31 +0300 Subject: [PATCH] Move Tube to a file --- src/bin/flat/main.rs | 78 ++----------------------------------- src/bin/flat/tube/metric.rs | 77 ++++++++++++++++++++++++++++++++++++ src/bin/flat/tube/mod.rs | 1 + 3 files changed, 81 insertions(+), 75 deletions(-) create mode 100644 src/bin/flat/tube/metric.rs create mode 100644 src/bin/flat/tube/mod.rs diff --git a/src/bin/flat/main.rs b/src/bin/flat/main.rs index 2f65d2a..b2c6acb 100644 --- a/src/bin/flat/main.rs +++ b/src/bin/flat/main.rs @@ -6,10 +6,12 @@ use glam::*; mod riemann; mod fns; mod float_fun; +mod tube; -use riemann::{Tens2, Decomp2, Metric, trace_iter}; +use riemann::{Metric, trace_iter}; use shape::Shape; use Subspace::{Boundary, Inner, Outer}; +use tube::metric::Tube; const DT: f32 = 0.1; @@ -607,77 +609,3 @@ impl FlatCell for TubeInside { (vec2(-self.tube.inner_radius, -self.tube.internal_halflength), vec2(self.tube.inner_radius, self.tube.internal_halflength)) } } - -#[derive(Copy, Clone, Debug)] -struct Tube { - outer_radius: f32, - inner_radius: f32, - external_halflength: f32, - internal_halflength: f32, -} - -impl Tube { - fn fx(&self) -> fns::LinearLimiter { fns::LinearLimiter { min: self.inner_radius, max: self.outer_radius } } - fn fy(&self) -> fns::QuadraticAccelerator { fns::QuadraticAccelerator { internal: self.internal_halflength, external: self.external_halflength } } - - pub fn y(&self, v: f32) -> f32 { self.fy().x(v) } - pub fn v(&self, y: f32) -> f32 { self.fy().u(y) } - pub fn dy(&self, v: f32) -> f32 { self.fy().dx(v) } - pub fn dv(&self, y: f32) -> f32 { self.fy().du(y) } -} - -impl Metric for Tube { - fn sqrt_at(&self, pos: Vec2) -> Decomp2 { - let sx = self.fx().value(pos.x); - let sy = self.fy().du(pos.y); - let s = sx + sy - sx * sy; - assert!(sx.is_finite()); - assert!(sy.is_finite()); - assert!(sy > 0.0); - Decomp2 { - ortho: Mat2::IDENTITY, - diag: vec2(1.0, s), - } - } - - fn part_derivs_at(&self, pos: Vec2) -> Tens2 { - let sx = self.fx().value(pos.x); - let sy = self.fy().du(pos.y); - let s = sx + sy - sx * sy; - let dsx_dx = self.fx().derivative(pos.x); - let dsy_dy = self.fy().d2u(pos.y); - 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_tube_metric_derivs() { - struct Approx(Tube); - impl Metric for Approx { - fn sqrt_at(&self, pos: Vec2) -> Decomp2 { self.0.sqrt_at(pos) } - } - let testee = Tube { - 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, "Bad derivative computation at {pos}:\n explicit: {computed:?}\n numerical: {reference:?}\n"); - } - } -} diff --git a/src/bin/flat/tube/metric.rs b/src/bin/flat/tube/metric.rs new file mode 100644 index 0000000..b64603b --- /dev/null +++ b/src/bin/flat/tube/metric.rs @@ -0,0 +1,77 @@ +use glam::{f32, Mat2, Vec2, vec2}; +use crate::fns; +use crate::riemann::{Decomp2, Metric, Tens2}; + +#[derive(Copy, Clone, Debug)] +pub struct Tube { + pub outer_radius: f32, + pub inner_radius: f32, + pub external_halflength: f32, + pub internal_halflength: f32, +} + +impl Tube { + fn fx(&self) -> fns::LinearLimiter { fns::LinearLimiter { min: self.inner_radius, max: self.outer_radius } } + fn fy(&self) -> fns::QuadraticAccelerator { fns::QuadraticAccelerator { internal: self.internal_halflength, external: self.external_halflength } } + + pub fn y(&self, v: f32) -> f32 { self.fy().x(v) } + pub fn v(&self, y: f32) -> f32 { self.fy().u(y) } + pub fn dy(&self, v: f32) -> f32 { self.fy().dx(v) } + pub fn dv(&self, y: f32) -> f32 { self.fy().du(y) } +} + +impl Metric for Tube { + fn sqrt_at(&self, pos: Vec2) -> Decomp2 { + let sx = self.fx().value(pos.x); + let sy = self.fy().du(pos.y); + let s = sx + sy - sx * sy; + assert!(sx.is_finite()); + assert!(sy.is_finite()); + assert!(sy > 0.0); + Decomp2 { + ortho: Mat2::IDENTITY, + diag: vec2(1.0, s), + } + } + + fn part_derivs_at(&self, pos: Vec2) -> Tens2 { + let sx = self.fx().value(pos.x); + let sy = self.fy().du(pos.y); + let s = sx + sy - sx * sy; + let dsx_dx = self.fx().derivative(pos.x); + let dsy_dy = self.fy().d2u(pos.y); + 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_tube_metric_derivs() { + struct Approx(Tube); + impl Metric for Approx { + fn sqrt_at(&self, pos: Vec2) -> Decomp2 { self.0.sqrt_at(pos) } + } + let testee = Tube { + 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, "Bad derivative computation at {pos}:\n explicit: {computed:?}\n numerical: {reference:?}\n"); + } + } +} diff --git a/src/bin/flat/tube/mod.rs b/src/bin/flat/tube/mod.rs new file mode 100644 index 0000000..e0026c2 --- /dev/null +++ b/src/bin/flat/tube/mod.rs @@ -0,0 +1 @@ +pub mod metric;