From 1de98d9550fcdc19cd35da8fafa58b9f8aa384a1 Mon Sep 17 00:00:00 2001 From: numzero Date: Mon, 10 Jun 2024 16:15:26 +0300 Subject: [PATCH] Extract a module --- src/bin/{flat.rs => flat/main.rs} | 154 +----------------------------- src/bin/flat/riemann.rs | 148 ++++++++++++++++++++++++++++ 2 files changed, 151 insertions(+), 151 deletions(-) rename src/bin/{flat.rs => flat/main.rs} (85%) create mode 100644 src/bin/flat/riemann.rs diff --git a/src/bin/flat.rs b/src/bin/flat/main.rs similarity index 85% rename from src/bin/flat.rs rename to src/bin/flat/main.rs index b5b870b..98db0bc 100644 --- a/src/bin/flat.rs +++ b/src/bin/flat/main.rs @@ -2,6 +2,9 @@ use std::f32::consts::{FRAC_PI_2, PI}; use flo_draw::*; use flo_canvas::*; use glam::*; + +mod riemann; + use riemann::{Tens2, Decomp2, Metric, trace_iter}; use shape::Shape; use Subspace::{Boundary, Inner, Outer}; @@ -686,154 +689,3 @@ fn test_rect_metric_derivs() { } } } - -mod riemann { - use glam::*; - - pub struct Decomp2 { - pub ortho: Mat2, - pub diag: Vec2, - } - - impl Decomp2 { - fn square(&self) -> Self { - Self { - ortho: self.ortho, - diag: self.diag * self.diag, - } - } - - pub(crate) fn inverse(&self) -> Self { - Self { - ortho: self.ortho, - diag: Vec2::splat(1.0) / self.diag, - } - } - } - - impl From for Mat2 { - fn from(value: Decomp2) -> Self { - value.ortho.transpose() * Mat2::from_diagonal(value.diag) * value.ortho - } - } - - pub type Tens2 = [Mat2; 2]; - - pub trait Metric { - fn sqrt_at(&self, pos: Vec2) -> Decomp2; - - fn at(&self, pos: Vec2) -> Mat2 { - self.sqrt_at(pos).square().into() - } - - fn inverse_at(&self, pos: Vec2) -> Mat2 { - self.sqrt_at(pos).square().inverse().into() - } - - fn part_derivs_at(&self, pos: Vec2) -> Tens2 { - part_deriv(|p| self.at(p), pos, 1.0 / 1024.0) // division by such eps is exact which is good for overall precision - } - - fn vec_length_at(&self, at: Vec2, v: Vec2) -> f32 { - v.dot(self.at(at) * v).sqrt() - } - - fn normalize_vec_at(&self, at: Vec2, v: Vec2) -> Vec2 { - v / self.vec_length_at(at, v) - } - - fn globalize(&self, at: Vec2, v: Vec2) -> Vec2 { - Mat2::from(self.sqrt_at(at).inverse()) * v - } - } - - pub struct TraceIter<'a, M: Metric> { - space: &'a M, - p: Vec2, - v: Vec2, - dt: f32, - } - - impl<'a, M: Metric> Iterator for TraceIter<'a, M> { - type Item = Vec2; - - fn next(&mut self) -> Option { - let a: Vec2 = -contract2(krist(self.space, self.p), self.v); - self.v = self.v + a * self.dt; - self.p = self.p + self.v * self.dt; - Some(self.p) - } - } - - pub fn trace_iter(space: &M, base: Vec2, dir: Vec2, dt: f32) -> TraceIter { - TraceIter { - space, - p: base, - v: space.normalize_vec_at(base, dir), - dt, - } - } - - pub fn krist(space: &impl Metric, pos: Vec2) -> Tens2 { - // Γ^i_k_l = .5 * g^i^m * (g_m_k,l + g_m_l,k - g_k_l,m) - let g = &space.inverse_at(pos); // с верхними индексами - let d = space.part_derivs_at(pos); - // ret[i][l][k] = sum((m) => .5f * g[m][i] * (d[k][l][m] + d[l][k][m] - d[m][k][l])) - make_tens2(|i, l, k| 0.5 * (0..2).map(|m| g.col(m)[i] * (d[l].col(k)[m] + d[k].col(m)[l] - d[m].col(k)[l])).sum::()) - } - - fn dir_deriv(f: impl Fn(Vec2) -> Mat2, pos: Vec2, delta: Vec2) -> Mat2 { - (f(pos + delta) - f(pos - delta)) / (2.0 * delta.length()) - } - - fn part_deriv(f: impl Fn(Vec2) -> Mat2, pos: Vec2, eps: f32) -> Tens2 { - [ - dir_deriv(&f, pos, vec2(eps, 0.0)), - dir_deriv(&f, pos, vec2(0.0, eps)), - ] - } - - /// Сворачивает тензор t с вектором u - pub fn contract(t: Tens2, u: Vec2) -> Mat2 { - mat2(t[0] * u, t[1] * u).transpose() - } - - /// Сворачивает тензор t с вектором v дважды, по второму и третьему индексам. - pub fn contract2(t: Tens2, v: Vec2) -> Vec2 { - contract(t, v) * v - } - - fn make_vec2(f: impl Fn(usize) -> f32) -> Vec2 { - Vec2::from_array(std::array::from_fn(|i| f(i))) - } - - fn make_mat2(f: impl Fn(usize, usize) -> f32) -> Mat2 { - Mat2::from_cols_array_2d(&std::array::from_fn(|i| std::array::from_fn(|j| f(i, j)))) - } - - fn make_tens2(f: impl Fn(usize, usize, usize) -> f32) -> Tens2 { - std::array::from_fn(|i| make_mat2(|j, k| f(i, j, k))) - } - - #[test] - fn m2() { - let m = make_mat2(|i, j| (i + 2 * j) as f32); - assert_eq!(m.col(0)[0], 0.0); - assert_eq!(m.col(1)[0], 1.0); - assert_eq!(m.col(0)[1], 2.0); - assert_eq!(m.col(1)[1], 3.0); - } - - #[test] - fn t2() { - let t = make_tens2(|i, j, k| (i + 2 * j + 4 * k) as f32); - assert_eq!(t[0].col(0)[0], 0.0); - assert_eq!(t[1].col(0)[0], 1.0); - assert_eq!(t[0].col(1)[0], 2.0); - assert_eq!(t[1].col(1)[0], 3.0); - assert_eq!(t[0].col(0)[1], 4.0); - assert_eq!(t[1].col(0)[1], 5.0); - assert_eq!(t[0].col(1)[1], 6.0); - assert_eq!(t[1].col(1)[1], 7.0); - } -} diff --git a/src/bin/flat/riemann.rs b/src/bin/flat/riemann.rs new file mode 100644 index 0000000..24b154f --- /dev/null +++ b/src/bin/flat/riemann.rs @@ -0,0 +1,148 @@ +use glam::*; + +pub struct Decomp2 { + pub ortho: Mat2, + pub diag: Vec2, +} + +impl Decomp2 { + fn square(&self) -> Self { + Self { + ortho: self.ortho, + diag: self.diag * self.diag, + } + } + + pub(crate) fn inverse(&self) -> Self { + Self { + ortho: self.ortho, + diag: Vec2::splat(1.0) / self.diag, + } + } +} + +impl From for Mat2 { + fn from(value: Decomp2) -> Self { + value.ortho.transpose() * Mat2::from_diagonal(value.diag) * value.ortho + } +} + +pub type Tens2 = [Mat2; 2]; + +pub trait Metric { + fn sqrt_at(&self, pos: Vec2) -> Decomp2; + + fn at(&self, pos: Vec2) -> Mat2 { + self.sqrt_at(pos).square().into() + } + + fn inverse_at(&self, pos: Vec2) -> Mat2 { + self.sqrt_at(pos).square().inverse().into() + } + + fn part_derivs_at(&self, pos: Vec2) -> Tens2 { + part_deriv(|p| self.at(p), pos, 1.0 / 1024.0) // division by such eps is exact which is good for overall precision + } + + fn vec_length_at(&self, at: Vec2, v: Vec2) -> f32 { + v.dot(self.at(at) * v).sqrt() + } + + fn normalize_vec_at(&self, at: Vec2, v: Vec2) -> Vec2 { + v / self.vec_length_at(at, v) + } + + fn globalize(&self, at: Vec2, v: Vec2) -> Vec2 { + Mat2::from(self.sqrt_at(at).inverse()) * v + } +} + +pub struct TraceIter<'a, M: Metric> { + space: &'a M, + p: Vec2, + v: Vec2, + dt: f32, +} + +impl<'a, M: Metric> Iterator for TraceIter<'a, M> { + type Item = Vec2; + + fn next(&mut self) -> Option { + let a: Vec2 = -contract2(krist(self.space, self.p), self.v); + self.v = self.v + a * self.dt; + self.p = self.p + self.v * self.dt; + Some(self.p) + } +} + +pub fn trace_iter(space: &M, base: Vec2, dir: Vec2, dt: f32) -> TraceIter { + TraceIter { + space, + p: base, + v: space.normalize_vec_at(base, dir), + dt, + } +} + +pub fn krist(space: &impl Metric, pos: Vec2) -> Tens2 { + // Γ^i_k_l = .5 * g^i^m * (g_m_k,l + g_m_l,k - g_k_l,m) + let g = &space.inverse_at(pos); // с верхними индексами + let d = space.part_derivs_at(pos); + // ret[i][l][k] = sum((m) => .5f * g[m][i] * (d[k][l][m] + d[l][k][m] - d[m][k][l])) + make_tens2(|i, l, k| 0.5 * (0..2).map(|m| g.col(m)[i] * (d[l].col(k)[m] + d[k].col(m)[l] - d[m].col(k)[l])).sum::()) +} + +fn dir_deriv(f: impl Fn(Vec2) -> Mat2, pos: Vec2, delta: Vec2) -> Mat2 { + (f(pos + delta) - f(pos - delta)) / (2.0 * delta.length()) +} + +fn part_deriv(f: impl Fn(Vec2) -> Mat2, pos: Vec2, eps: f32) -> Tens2 { + [ + dir_deriv(&f, pos, vec2(eps, 0.0)), + dir_deriv(&f, pos, vec2(0.0, eps)), + ] +} + +/// Сворачивает тензор t с вектором u +pub fn contract(t: Tens2, u: Vec2) -> Mat2 { + mat2(t[0] * u, t[1] * u).transpose() +} + +/// Сворачивает тензор t с вектором v дважды, по второму и третьему индексам. +pub fn contract2(t: Tens2, v: Vec2) -> Vec2 { + contract(t, v) * v +} + +fn make_vec2(f: impl Fn(usize) -> f32) -> Vec2 { + Vec2::from_array(std::array::from_fn(|i| f(i))) +} + +fn make_mat2(f: impl Fn(usize, usize) -> f32) -> Mat2 { + Mat2::from_cols_array_2d(&std::array::from_fn(|i| std::array::from_fn(|j| f(i, j)))) +} + +fn make_tens2(f: impl Fn(usize, usize, usize) -> f32) -> Tens2 { + std::array::from_fn(|i| make_mat2(|j, k| f(i, j, k))) +} + +#[test] +fn m2() { + let m = make_mat2(|i, j| (i + 2 * j) as f32); + assert_eq!(m.col(0)[0], 0.0); + assert_eq!(m.col(1)[0], 1.0); + assert_eq!(m.col(0)[1], 2.0); + assert_eq!(m.col(1)[1], 3.0); +} + +#[test] +fn t2() { + let t = make_tens2(|i, j, k| (i + 2 * j + 4 * k) as f32); + assert_eq!(t[0].col(0)[0], 0.0); + assert_eq!(t[1].col(0)[0], 1.0); + assert_eq!(t[0].col(1)[0], 2.0); + assert_eq!(t[1].col(1)[0], 3.0); + assert_eq!(t[0].col(0)[1], 4.0); + assert_eq!(t[1].col(0)[1], 5.0); + assert_eq!(t[0].col(1)[1], 6.0); + assert_eq!(t[1].col(1)[1], 7.0); +}