From e57692587a4dfee5c72969f862c5d662bf964cbf Mon Sep 17 00:00:00 2001 From: numzero Date: Thu, 7 Nov 2024 02:42:52 +0300 Subject: [PATCH] Make the tube actually cylindrical --- src/tube/coords.rs | 43 ++++++++++++++++++++----------------------- src/tube/metric.rs | 25 ++++++++++++++++--------- src/tube/mod.rs | 8 +++++--- 3 files changed, 41 insertions(+), 35 deletions(-) diff --git a/src/tube/coords.rs b/src/tube/coords.rs index 6d68127..657420b 100644 --- a/src/tube/coords.rs +++ b/src/tube/coords.rs @@ -3,7 +3,7 @@ use glam::{vec3, Mat3, Vec3}; use crate::riemann::Metric; use crate::types::{Location, Ray}; -use super::{Rect, Tube}; +use super::{Tube, YCylinder}; pub trait FlatCoordinateSystem { fn flat_to_global(&self, v: T) -> T; @@ -80,8 +80,9 @@ impl FlatCoordinateSystem for InnerCS { impl FlatRegion for InnerCS { fn distance_to_boundary(&self, ray: Ray) -> Option { - Rect { - size: vec3(self.0.inner_radius, self.0.internal_halflength, self.0.inner_radius), + YCylinder { + radius: self.0.inner_radius, + half_length: self.0.internal_halflength, } .trace_out_of(ray) } @@ -97,12 +98,9 @@ impl MetricCS for OuterCS { impl FlatCoordinateSystem for OuterCS { fn flat_to_global(&self, pos: Vec3) -> Vec3 { - let inner = Rect { - size: vec3( - self.0.inner_radius + 1.0, - self.0.external_halflength, - self.0.inner_radius + 1.0, - ), + let inner = YCylinder { + radius: self.0.inner_radius + 1.0, + half_length: self.0.external_halflength, }; if inner.is_inside(pos) { let Vec3 { x, y: v, z } = pos; @@ -116,12 +114,9 @@ impl FlatCoordinateSystem for OuterCS { } fn global_to_flat(&self, pos: Vec3) -> Vec3 { - let inner = Rect { - size: vec3( - self.0.inner_radius + 1.0, - self.0.external_halflength, - self.0.inner_radius + 1.0, - ), + let inner = YCylinder { + radius: self.0.inner_radius + 1.0, + half_length: self.0.external_halflength, }; if inner.is_inside(pos) { let Vec3 { x: u, y, z: w } = pos; // в основной СК @@ -135,8 +130,9 @@ impl FlatCoordinateSystem for OuterCS { impl FlatRegion for OuterCS { fn distance_to_boundary(&self, ray: Ray) -> Option { - Rect { - size: vec3(self.0.outer_radius, self.0.external_halflength, self.0.outer_radius), + YCylinder { + radius: self.0.outer_radius, + half_length: self.0.external_halflength, } .trace_into(ray) } @@ -346,15 +342,16 @@ mod tests { external_halflength: 300.0, }); // TODO replace 200.20016 with something sane + // NOTE it can’t test −30..30 as that area intersects the boundary test_flat_region( &mapper, - (vec3(-30.0, -300.0, -30.0), vec3(30.0, -1.0, 30.0)), - (vec3(-30.0, -300.0, -30.0), vec3(30.0, -200.20016, 30.0)), + (vec3(-20.0, -300.0, -20.0), vec3(20.0, -1.0, 20.0)), + (vec3(-20.0, -300.0, -20.0), vec3(20.0, -200.20016, 20.0)), ); test_flat_region( &mapper, - (vec3(-30.0, 1.0, -30.0), vec3(30.0, 300.0, 30.0)), - (vec3(-30.0, 200.20016, -30.0), vec3(30.0, 300.0, 30.0)), + (vec3(-20.0, 1.0, -20.0), vec3(20.0, 300.0, 20.0)), + (vec3(-20.0, 200.20016, -20.0), vec3(20.0, 300.0, 20.0)), ); test_flat_region( &mapper, @@ -451,9 +448,9 @@ mod tests { } } // accelerating - for x in linspace(-29., 29., 20) { + for x in linspace(-21., 21., 20) { for y in linspace(1., 299., 20) { - for z in linspace(-29., 29., 20) { + for z in linspace(-21., 21., 20) { let v = mapper .global_to_flat(Location { pos: vec3(x, y, z), diff --git a/src/tube/metric.rs b/src/tube/metric.rs index ba65a81..9a303e8 100644 --- a/src/tube/metric.rs +++ b/src/tube/metric.rs @@ -42,10 +42,11 @@ impl Tube { impl Metric for Tube { fn sqrt_at(&self, pos: Vec3) -> Decomp3 { - let sx = self.fx().value(pos.x); + let r = f32::hypot(pos.x, pos.z); + let sr = self.fx().value(r); let sy = self.fy().du(pos.y); - let s = sx + sy - sx * sy; - assert!(sx.is_finite()); + let s = sr + sy - sr * sy; + assert!(sr.is_finite()); assert!(sy.is_finite()); assert!(sy > 0.0); Decomp3 { @@ -55,17 +56,23 @@ impl Metric for Tube { } fn part_derivs_at(&self, pos: Vec3) -> Tens3 { - let sx = self.fx().value(pos.x); + let r = f32::hypot(pos.x, pos.z); + let sr = self.fx().value(r); let sy = self.fy().du(pos.y); - let s = sx + sy - sx * sy; - let dsx_dx = self.fx().derivative(pos.x); + let s = sr + sy - sr * sy; + let dsr_dr = self.fx().derivative(r); 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; + let ds2_dr = 2.0 * s * (1.0 - sy) * dsr_dr; + let ds2_dy = 2.0 * s * (1.0 - sr) * dsy_dy; + let (ds2_dx, ds2_dz) = if ds2_dr.abs() < 1e-5 { + (0., 0.) + } else { + (ds2_dr * pos.x / r, ds2_dr * pos.z / r) + }; [ Mat3::from_cols_array(&[0., 0., 0., 0., ds2_dx, 0., 0., 0., 0.]), Mat3::from_cols_array(&[0., 0., 0., 0., ds2_dy, 0., 0., 0., 0.]), - Mat3::from_cols_array(&[0., 0., 0., 0., 0., 0., 0., 0., 0.]), + Mat3::from_cols_array(&[0., 0., 0., 0., ds2_dz, 0., 0., 0., 0.]), ] } } diff --git a/src/tube/mod.rs b/src/tube/mod.rs index cdd13da..7bd5a70 100644 --- a/src/tube/mod.rs +++ b/src/tube/mod.rs @@ -28,10 +28,12 @@ pub enum Subspace { impl Space { fn which_subspace(&self, pt: Vec3) -> Subspace { if pt.y.abs() > self.tube.external_halflength { + return Outer; + } + let r = f32::hypot(pt.x, pt.z); + if r > self.tube.outer_radius { Outer - } else if pt.x.abs() > self.tube.outer_radius { - Outer - } else if pt.x.abs() > self.tube.inner_radius { + } else if r > self.tube.inner_radius { Boundary } else { Inner