Split into modules
This commit is contained in:
parent
4c29c8de0b
commit
10563d8d56
225
src/bin/flat.rs
225
src/bin/flat.rs
|
|
@ -2,6 +2,7 @@ use flo_draw::*;
|
|||
use flo_canvas::*;
|
||||
use glm::*;
|
||||
use num_traits::identities::Zero;
|
||||
use riemann::{Decomp2, Metric, trace_iter};
|
||||
|
||||
pub fn main() {
|
||||
let space = Coil {
|
||||
|
|
@ -75,126 +76,138 @@ impl Metric for Coil {
|
|||
}
|
||||
}
|
||||
|
||||
struct Decomp2 {
|
||||
ortho: Mat2,
|
||||
diag: Vec2,
|
||||
}
|
||||
mod riemann {
|
||||
use glm::*;
|
||||
use num_traits::identities::Zero;
|
||||
|
||||
type Tens2 = [Mat2; 2];
|
||||
|
||||
trait Metric {
|
||||
fn halfmetric(&self, pos: Vec2) -> Decomp2;
|
||||
|
||||
fn metric(&self, pos: Vec2) -> Mat2 {
|
||||
let h = self.halfmetric(pos);
|
||||
transpose(&h.ortho) * diagonal(h.diag * h.diag) * h.ortho
|
||||
pub struct Decomp2 {
|
||||
pub ortho: Mat2,
|
||||
pub diag: Vec2,
|
||||
}
|
||||
|
||||
fn dmetric(&self, pos: Vec2) -> Tens2 {
|
||||
part_deriv(|p| self.metric(p), pos, 1.0e-3)
|
||||
}
|
||||
type Tens2 = [Mat2; 2];
|
||||
|
||||
fn length(&self, at: Vec2, v: Vec2) -> f32 {
|
||||
sqrt(dot(v, self.metric(at) * v))
|
||||
}
|
||||
pub trait Metric {
|
||||
fn halfmetric(&self, pos: Vec2) -> Decomp2;
|
||||
|
||||
fn normalize(&self, at: Vec2, v: Vec2) -> Vec2 {
|
||||
v / self.length(at, v)
|
||||
}
|
||||
fn globalize(&self, at: Vec2, v: Vec2) -> Vec2 {
|
||||
let h = self.halfmetric(at);
|
||||
transpose(&h.ortho) * diagonal( Vec2::from_s(1.0) / h.diag) * h.ortho * v
|
||||
}
|
||||
}
|
||||
fn metric(&self, pos: Vec2) -> Mat2 {
|
||||
let h = self.halfmetric(pos);
|
||||
transpose(&h.ortho) * diagonal(h.diag * h.diag) * h.ortho
|
||||
}
|
||||
|
||||
struct TraceIter<'a, M: Metric> {
|
||||
space: &'a M,
|
||||
p: Vec2,
|
||||
v: Vec2,
|
||||
dt: f32,
|
||||
}
|
||||
fn dmetric(&self, pos: Vec2) -> Tens2 {
|
||||
part_deriv(|p| self.metric(p), pos, 1.0e-3)
|
||||
}
|
||||
|
||||
impl<'a, M: Metric> Iterator for TraceIter<'a, M> {
|
||||
type Item = Vec2;
|
||||
fn length(&self, at: Vec2, v: Vec2) -> f32 {
|
||||
sqrt(dot(v, self.metric(at) * v))
|
||||
}
|
||||
|
||||
fn next(&mut self) -> Option<Self::Item> {
|
||||
let a: Vec2 = -convolute(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)
|
||||
}
|
||||
}
|
||||
|
||||
fn trace(space: &impl Metric, base: Vec2, dir: Vec2, distance: f32, dt: f32) -> Vec<Vec2> {
|
||||
trace_iter(space, base, dir, dt).take((distance / dt) as usize).collect()
|
||||
}
|
||||
|
||||
fn trace_iter<M: Metric>(space: &M, base: Vec2, dir: Vec2, dt: f32) -> TraceIter<M> {
|
||||
TraceIter{
|
||||
space: space,
|
||||
p: base,
|
||||
v: space.normalize(base, dir),
|
||||
dt: dt,
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn t_iter() {
|
||||
let space = Coil {
|
||||
coil_scale: 2.0,
|
||||
coil_r: 300.0,
|
||||
coil_w: 50.0,
|
||||
coil_m: 10.0,
|
||||
};
|
||||
let base = vec2(-500.0, 0.0);
|
||||
let dir = vec2(1.0, 0.3);
|
||||
let dt = 1.0;
|
||||
let steps = 1000;
|
||||
let a = trace(&space, base, dir, dt * (steps as f32), dt);
|
||||
let b: Vec<Vec2> = trace_iter(&space, base, dir, dt).take(steps).collect();
|
||||
assert_eq!(a, b);
|
||||
}
|
||||
|
||||
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 = inverse(&space.metric(pos)); // с верхними индексами
|
||||
let d = space.dmetric(pos);
|
||||
let mut ret: Tens2 = [Mat2::zero(); 2];
|
||||
// ret[i][l][k] = sum((m) => .5f * g[m][i] * (d[k][l][m] + d[l][k][m] - d[m][k][l]))
|
||||
for i in 0..2 {
|
||||
for l in 0..2 {
|
||||
for k in 0..2 {
|
||||
let mut v = 0.0;
|
||||
for m in 0..2 {
|
||||
v += g[m][i] * (d[l][k][m] + d[k][m][l] - d[m][k][l]);
|
||||
}
|
||||
ret[i][l][k] = 0.5 * v;
|
||||
}
|
||||
fn normalize(&self, at: Vec2, v: Vec2) -> Vec2 {
|
||||
v / self.length(at, v)
|
||||
}
|
||||
fn globalize(&self, at: Vec2, v: Vec2) -> Vec2 {
|
||||
let h = self.halfmetric(at);
|
||||
transpose(&h.ortho) * diagonal(Vec2::from_s(1.0) / h.diag) * h.ortho * v
|
||||
}
|
||||
}
|
||||
ret
|
||||
}
|
||||
|
||||
fn dir_deriv(f: impl Fn(Vec2) -> Mat2, pos: Vec2, delta: Vec2) -> Mat2 {
|
||||
(f(pos + delta) - f(pos - delta)) / (2.0 * length(delta))
|
||||
}
|
||||
pub struct TraceIter<'a, M: Metric> {
|
||||
space: &'a M,
|
||||
p: Vec2,
|
||||
v: Vec2,
|
||||
dt: f32,
|
||||
}
|
||||
|
||||
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)),
|
||||
]
|
||||
}
|
||||
impl<'a, M: Metric> Iterator for TraceIter<'a, M> {
|
||||
type Item = Vec2;
|
||||
|
||||
fn convolute(G: Tens2, v: Vec2) -> Vec2 {
|
||||
vec2(
|
||||
dot(v, G[0] * v),
|
||||
dot(v, G[1] * v)
|
||||
)
|
||||
}
|
||||
fn next(&mut self) -> Option<Self::Item> {
|
||||
let a: Vec2 = -convolute(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)
|
||||
}
|
||||
}
|
||||
|
||||
fn diagonal(v: Vec2) -> Mat2 {
|
||||
mat2(v.x, 0.0, 0.0, v.y)
|
||||
fn trace(space: &impl Metric, base: Vec2, dir: Vec2, distance: f32, dt: f32) -> Vec<Vec2> {
|
||||
trace_iter(space, base, dir, dt).take((distance / dt) as usize).collect()
|
||||
}
|
||||
|
||||
pub fn trace_iter<M: Metric>(space: &M, base: Vec2, dir: Vec2, dt: f32) -> TraceIter<M> {
|
||||
TraceIter {
|
||||
space: space,
|
||||
p: base,
|
||||
v: space.normalize(base, dir),
|
||||
dt: dt,
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod test {
|
||||
use glm::*;
|
||||
use crate::riemann::{trace, trace_iter};
|
||||
use crate::Coil;
|
||||
|
||||
#[test]
|
||||
fn t_iter() {
|
||||
let space = Coil {
|
||||
coil_scale: 2.0,
|
||||
coil_r: 300.0,
|
||||
coil_w: 50.0,
|
||||
coil_m: 10.0,
|
||||
};
|
||||
let base = vec2(-500.0, 0.0);
|
||||
let dir = vec2(1.0, 0.3);
|
||||
let dt = 1.0;
|
||||
let steps = 1000;
|
||||
let a = trace(&space, base, dir, dt * (steps as f32), dt);
|
||||
let b: Vec<Vec2> = trace_iter(&space, base, dir, dt).take(steps).collect();
|
||||
assert_eq!(a, b);
|
||||
}
|
||||
}
|
||||
|
||||
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 = inverse(&space.metric(pos)); // с верхними индексами
|
||||
let d = space.dmetric(pos);
|
||||
let mut ret: Tens2 = [Mat2::zero(); 2];
|
||||
// ret[i][l][k] = sum((m) => .5f * g[m][i] * (d[k][l][m] + d[l][k][m] - d[m][k][l]))
|
||||
for i in 0..2 {
|
||||
for l in 0..2 {
|
||||
for k in 0..2 {
|
||||
let mut v = 0.0;
|
||||
for m in 0..2 {
|
||||
v += g[m][i] * (d[l][k][m] + d[k][m][l] - d[m][k][l]);
|
||||
}
|
||||
ret[i][l][k] = 0.5 * v;
|
||||
}
|
||||
}
|
||||
}
|
||||
ret
|
||||
}
|
||||
|
||||
fn dir_deriv(f: impl Fn(Vec2) -> Mat2, pos: Vec2, delta: Vec2) -> Mat2 {
|
||||
(f(pos + delta) - f(pos - delta)) / (2.0 * length(delta))
|
||||
}
|
||||
|
||||
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)),
|
||||
]
|
||||
}
|
||||
|
||||
fn convolute(G: Tens2, v: Vec2) -> Vec2 {
|
||||
vec2(
|
||||
dot(v, G[0] * v),
|
||||
dot(v, G[1] * v)
|
||||
)
|
||||
}
|
||||
|
||||
fn diagonal(v: Vec2) -> Mat2 {
|
||||
mat2(v.x, 0.0, 0.0, v.y)
|
||||
}
|
||||
}
|
||||
|
||||
fn sqr(x: f32) -> f32 {
|
||||
|
|
|
|||
Loading…
Reference in New Issue
Block a user