diffgeom/
metric.rs

1use super::coordinates::{CoordinateSystem, Point};
2use super::tensors::{
3    ContravariantIndex, CovariantIndex, InnerProduct, InvTwoForm, Tensor, TwoForm,
4};
5use crate::inner;
6use crate::typenum::consts::{U0, U1, U2, U3};
7use crate::typenum::{Exp, Pow, Unsigned};
8use generic_array::ArrayLength;
9
10/// Trait representing the metric properties of the coordinate system
11pub trait MetricSystem: CoordinateSystem
12where
13    <Self as CoordinateSystem>::Dimension: Pow<U2> + Pow<U3>,
14    Exp<<Self as CoordinateSystem>::Dimension, U2>: ArrayLength<f64>,
15    Exp<<Self as CoordinateSystem>::Dimension, U3>: ArrayLength<f64>,
16{
17    /// Returns the metric tensor at a given point.
18    fn g(point: &Point<Self>) -> TwoForm<Self>;
19
20    /// Returns the inverse metric tensor at a given point.
21    ///
22    /// The default implementation calculates the metric and then inverts it. A direct
23    /// implementation may be desirable for more performance.
24    fn inv_g(point: &Point<Self>) -> InvTwoForm<Self> {
25        Self::g(point).inverse().unwrap()
26    }
27
28    /// Returns the partial derivatives of the metric at a given point.
29    ///
30    /// The default implementation calculates them numerically. A direct implementation
31    /// may be desirable for performance.
32    fn dg(point: &Point<Self>) -> Tensor<Self, (CovariantIndex, (CovariantIndex, CovariantIndex))> {
33        let d = Self::dimension();
34        let mut result = Tensor::zero(point.clone());
35        let h = Self::small(point);
36
37        for j in 0..d {
38            let mut x = point.clone();
39            x[j] = x[j] - h;
40            let g1 = Self::g(&x);
41
42            x[j] = x[j] + h * 2.0;
43            let g2 = Self::g(&x);
44
45            for coord in g1.iter_coords() {
46                // calculate dg_i/dx^j
47                let index = [coord[0], coord[1], j];
48                result[&index[..]] = (g2[&*coord] - g1[&*coord]) / (2.0 * h);
49            }
50        }
51
52        result
53    }
54
55    /// Returns the covariant Christoffel symbols (with three lower indices).
56    ///
57    /// The default implementation calculates them from the metric. A direct implementation
58    /// may be desirable for performance.
59    fn covariant_christoffel(
60        point: &Point<Self>,
61    ) -> Tensor<Self, (CovariantIndex, (CovariantIndex, CovariantIndex))> {
62        let dg = Self::dg(point);
63        let mut result =
64            Tensor::<Self, (CovariantIndex, (CovariantIndex, CovariantIndex))>::zero(point.clone());
65
66        for i in result.iter_coords() {
67            result[&*i] =
68                0.5 * (dg[&*i] + dg[&[i[0], i[2], i[1]][..]] - dg[&[i[1], i[2], i[0]][..]]);
69        }
70
71        result
72    }
73
74    /// Returns the Christoffel symbols.
75    ///
76    /// The default implementation calculates them from the metric. A direct implementation
77    /// may be desirable for performance.
78    fn christoffel(
79        point: &Point<Self>,
80    ) -> Tensor<Self, (ContravariantIndex, (CovariantIndex, CovariantIndex))> {
81        let ig = Self::inv_g(point);
82        let gamma = Self::covariant_christoffel(point);
83
84        <InvTwoForm<Self> as InnerProduct<
85            Tensor<Self, (CovariantIndex, (CovariantIndex, CovariantIndex))>,
86            U1,
87            U2,
88        >>::inner_product(ig, gamma)
89    }
90}
91
92impl<T> Tensor<T, ContravariantIndex>
93where
94    T: MetricSystem,
95    T::Dimension: Pow<U1> + Pow<U2> + Pow<U3> + Unsigned,
96    Exp<T::Dimension, U1>: ArrayLength<f64>,
97    Exp<T::Dimension, U2>: ArrayLength<f64>,
98    Exp<T::Dimension, U3>: ArrayLength<f64>,
99{
100    pub fn square(&self) -> f64 {
101        let g = T::g(self.get_point());
102        let temp = inner!(_, _; U1, U2; g, self.clone());
103        *inner!(_, _; U0, U1; temp, self.clone())
104    }
105
106    pub fn normalize(&mut self) {
107        let len = self.square().abs().sqrt();
108        for i in 0..T::Dimension::to_usize() {
109            self[i] /= len;
110        }
111    }
112}
113
114impl<T> Tensor<T, CovariantIndex>
115where
116    T: MetricSystem,
117    T::Dimension: Pow<U1> + Pow<U2> + Pow<U3> + Unsigned,
118    Exp<T::Dimension, U1>: ArrayLength<f64>,
119    Exp<T::Dimension, U2>: ArrayLength<f64>,
120    Exp<T::Dimension, U3>: ArrayLength<f64>,
121{
122    pub fn square(&self) -> f64 {
123        let g = T::inv_g(self.get_point());
124        let temp = inner!(_, _; U1, U2; g, self.clone());
125        *inner!(_, _; U0, U1; temp, self.clone())
126    }
127
128    pub fn normalize(&mut self) {
129        let len = self.square().abs().sqrt();
130        for i in 0..T::Dimension::to_usize() {
131            self[i] /= len;
132        }
133    }
134}
diffgeom/metric.rs

diffgeom/
metric.rs
