conspire/math/integrate/field/
mod.rs

1mod adaptive;
2mod euler;
3mod hermite;
4mod rkmk;
5mod state;
6#[cfg(test)]
7mod test;
8
9use crate::math::{
10    Tensor, TensorError, TensorRank2, TensorTuple, TensorVector, integrate::IntegrationError,
11};
12use crate::units::Dimensionless;
13use std::{
14    marker::PhantomData,
15    ops::{Add, Mul},
16};
17
18pub use adaptive::{
19    integrate_rkmk, integrate_rkmk_adaptive, integrate_rkmk_dae_adaptive,
20    integrate_rkmk_dae_adaptive_first_order_root,
21    integrate_rkmk_dae_adaptive_second_order_minimize,
22};
23pub use euler::integrate_euler;
24pub use hermite::{HermiteSegment, interpolate_hermite};
25pub use rkmk::{
26    rkmk_dae_step, rkmk_dae_step_first_order_root, rkmk_dae_step_second_order_minimize, rkmk_step,
27};
28pub use state::{
29    EvolvedIncrement, EvolvedState, StateEvolution, integrate_rkmk_state,
30    integrate_rkmk_state_adaptive,
31};
32
33const RECONSTRUCT_FAILED: &str =
34    "the field increment has no reconstruction (matrix exponential undefined)";
35
36/// The geometry of one integrated state field: how an increment advances the state.
37///
38/// [`Self::Increment`] is an element of the field's tangent space (its Lie algebra
39/// for a group-valued field). It equals [`Self::Point`] for a flat field, but not
40/// in general — e.g. `F_p` is a `Reference → Intermediate` map while its algebra
41/// element `D_p Δt` maps `Intermediate → Intermediate`.
42pub trait Integrable {
43    /// The state value this field carries.
44    type Point: Tensor;
45    /// The tangent/algebra element that advances a [`Self::Point`].
46    type Increment: Tensor;
47    /// Advances `base` by `increment`.
48    fn reconstruct(
49        base: &Self::Point,
50        increment: &Self::Increment,
51    ) -> Result<Self::Point, TensorError>;
52    /// The RKMK correction: maps a rate-scaled increment to an algebra increment
53    /// at the accumulated algebra element `sigma`. Flat fields are the identity.
54    fn dexpinv(_sigma: &Self::Increment, increment: Self::Increment) -> Self::Increment {
55        increment
56    }
57}
58
59/// A state in a flat vector space: the increment simply adds.
60pub struct Flat<T>(PhantomData<T>);
61
62impl<T> Integrable for Flat<T>
63where
64    T: Clone + Tensor,
65    for<'a> T: Add<&'a T, Output = T>,
66{
67    type Point = T;
68    type Increment = T;
69    fn reconstruct(base: &T, increment: &T) -> Result<T, TensorError> {
70        Ok(base.clone() + increment)
71    }
72}
73
74/// A state acted on by the matrix exponential, `X_{n+1} = exp(increment) X_n`,
75/// staying on the unimodular group (`det = 1`) whenever the increment is
76/// trace-free. The state maps `B → A` while its algebra element maps `A → A`,
77/// so `F_p` (`Reference → Intermediate`) is `Unimodular<Intermediate, Reference>`.
78pub struct Unimodular<A, B = A>(PhantomData<(A, B)>);
79
80impl<A, B> Integrable for Unimodular<A, B>
81where
82    TensorRank2<3, A, B, Dimensionless>: Tensor,
83    TensorRank2<3, A, A, Dimensionless>: Tensor,
84    for<'a> TensorRank2<3, A, A, Dimensionless>:
85        Mul<&'a TensorRank2<3, A, B, Dimensionless>, Output = TensorRank2<3, A, B, Dimensionless>>,
86{
87    type Point = TensorRank2<3, A, B, Dimensionless>;
88    type Increment = TensorRank2<3, A, A, Dimensionless>;
89    fn reconstruct(
90        base: &Self::Point,
91        increment: &Self::Increment,
92    ) -> Result<Self::Point, TensorError> {
93        Ok(increment.expm()? * base)
94    }
95    fn dexpinv(sigma: &Self::Increment, increment: Self::Increment) -> Self::Increment {
96        sigma.dexpinv(&increment)
97    }
98}
99
100/// A composite of two fields; its state is the matching [`TensorTuple`], and an
101/// increment reconstructs component-wise. Nests right for three or more fields.
102pub struct Product<H, T>(PhantomData<(H, T)>);
103
104impl<H, T> Integrable for Product<H, T>
105where
106    H: Integrable,
107    T: Integrable,
108    TensorTuple<H::Point, T::Point>: Tensor,
109    TensorTuple<H::Increment, T::Increment>: Tensor,
110{
111    type Point = TensorTuple<H::Point, T::Point>;
112    type Increment = TensorTuple<H::Increment, T::Increment>;
113    fn reconstruct(
114        base: &Self::Point,
115        increment: &Self::Increment,
116    ) -> Result<Self::Point, TensorError> {
117        Ok(TensorTuple(
118            H::reconstruct(&base.0, &increment.0)?,
119            T::reconstruct(&base.1, &increment.1)?,
120        ))
121    }
122    fn dexpinv(sigma: &Self::Increment, increment: Self::Increment) -> Self::Increment {
123        TensorTuple(
124            H::dexpinv(&sigma.0, increment.0),
125            T::dexpinv(&sigma.1, increment.1),
126        )
127    }
128}
129
130/// A list of independent copies of one field, e.g. every Gauss point's plastic
131/// state across a mesh; an increment reconstructs entry-wise. Composes with
132/// [`Product`] for a multi-block mesh (`Product<List<Field1>, List<Field2>>`).
133pub struct List<Field>(PhantomData<Field>);
134
135impl<Field> Integrable for List<Field>
136where
137    Field: Integrable,
138    TensorVector<Field::Point>: Tensor<Item = Field::Point>,
139    TensorVector<Field::Increment>: Tensor<Item = Field::Increment>,
140{
141    type Point = TensorVector<Field::Point>;
142    type Increment = TensorVector<Field::Increment>;
143    fn reconstruct(
144        base: &Self::Point,
145        increment: &Self::Increment,
146    ) -> Result<Self::Point, TensorError> {
147        base.iter()
148            .zip(increment.iter())
149            .map(|(base, increment)| Field::reconstruct(base, increment))
150            .collect()
151    }
152    fn dexpinv(sigma: &Self::Increment, increment: Self::Increment) -> Self::Increment {
153        sigma
154            .iter()
155            .zip(increment)
156            .map(|(sigma, increment)| Field::dexpinv(sigma, increment))
157            .collect()
158    }
159}
160
161fn reconstruct_or_err<Field: Integrable>(
162    base: &Field::Point,
163    increment: &Field::Increment,
164) -> Result<Field::Point, IntegrationError> {
165    Field::reconstruct(base, increment)
166        .map_err(|_| IntegrationError::from(RECONSTRUCT_FAILED.to_string()))
167}
conspire/math/integrate/field/mod.rs

conspire/math/integrate/field/
mod.rs
