conspire/domain/solid/elastic_viscoplastic/
mod.rs

1use crate::{
2    domain::{
3        Blocks, ElasticViscoplasticAndElastic, ElementModel, ElementModelError, Model,
4        NodalCoordinates, NodalCoordinatesHistory,
5        block::element::Elements,
6        solid::{NodalForcesSolid, NodalStiffnessesSolid, elastic::ElasticElements},
7    },
8    math::{
9        Derivative, Differentiable, Quantity, Scalar, Tensor, TensorTuple, TensorTupleVec,
10        TensorVec, TensorVector,
11        integrate::{
12            ButcherTableau, EmbeddedTableau, Integrable, IntegrationError, Product,
13            integrate_rkmk_dae_adaptive_first_order_root, rkmk_dae_step_first_order_root,
14        },
15        optimize::{EqualityConstraint, FirstOrderRootFinding},
16    },
17    mechanics::Times,
18    units::Time,
19};
20use std::ops::Mul;
21
22pub type ElasticViscoplasticBCs = fn(Quantity<Time>) -> EqualityConstraint;
23
24pub trait ElasticViscoplasticElements<S, const D: usize>
25where
26    Self: Elements,
27    S: Differentiable,
28{
29    fn initial_state(&self) -> S;
30    fn nodal_forces_into(
31        &self,
32        nodal_coordinates: &NodalCoordinates<D>,
33        state_variables: &S,
34        nodal_forces: &mut NodalForcesSolid<D>,
35    ) -> Result<(), ElementModelError>;
36    fn nodal_forces(
37        &self,
38        nodal_coordinates: &NodalCoordinates<D>,
39        state_variables: &S,
40    ) -> Result<NodalForcesSolid<D>, ElementModelError> {
41        let mut nodal_forces = NodalForcesSolid::zero(nodal_coordinates.len());
42        self.nodal_forces_into(nodal_coordinates, state_variables, &mut nodal_forces)?;
43        Ok(nodal_forces)
44    }
45    fn nodal_stiffnesses_into(
46        &self,
47        nodal_coordinates: &NodalCoordinates<D>,
48        state_variables: &S,
49        nodal_stiffnesses: &mut NodalStiffnessesSolid<D>,
50    ) -> Result<(), ElementModelError>;
51    fn nodal_stiffnesses(
52        &self,
53        nodal_coordinates: &NodalCoordinates<D>,
54        state_variables: &S,
55    ) -> Result<NodalStiffnessesSolid<D>, ElementModelError> {
56        let mut nodal_stiffnesses = NodalStiffnessesSolid::zero(nodal_coordinates.len());
57        self.nodal_stiffnesses_into(nodal_coordinates, state_variables, &mut nodal_stiffnesses)?;
58        Ok(nodal_stiffnesses)
59    }
60    fn state_variables_evolution(
61        &self,
62        nodal_coordinates: &NodalCoordinates<D>,
63        state_variables: &S,
64    ) -> Result<Derivative<S>, ElementModelError>;
65}
66
67impl<B, S, const D: usize> ElasticViscoplasticElements<S, D> for Model<B, D>
68where
69    B: ElasticViscoplasticElements<S, D>,
70    S: Differentiable,
71{
72    fn initial_state(&self) -> S {
73        self.blocks.initial_state()
74    }
75    fn nodal_forces_into(
76        &self,
77        nodal_coordinates: &NodalCoordinates<D>,
78        state_variables: &S,
79        nodal_forces: &mut NodalForcesSolid<D>,
80    ) -> Result<(), ElementModelError> {
81        self.blocks
82            .nodal_forces_into(nodal_coordinates, state_variables, nodal_forces)
83    }
84    fn nodal_stiffnesses_into(
85        &self,
86        nodal_coordinates: &NodalCoordinates<D>,
87        state_variables: &S,
88        nodal_stiffnesses: &mut NodalStiffnessesSolid<D>,
89    ) -> Result<(), ElementModelError> {
90        self.blocks
91            .nodal_stiffnesses_into(nodal_coordinates, state_variables, nodal_stiffnesses)
92    }
93    fn state_variables_evolution(
94        &self,
95        nodal_coordinates: &NodalCoordinates<D>,
96        state_variables: &S,
97    ) -> Result<Derivative<S>, ElementModelError> {
98        self.blocks
99            .state_variables_evolution(nodal_coordinates, state_variables)
100    }
101}
102
103impl<B1, B2, S, const D: usize> ElasticViscoplasticElements<S, D>
104    for ElasticViscoplasticAndElastic<B1, B2>
105where
106    B1: ElasticViscoplasticElements<S, D>,
107    B2: ElasticElements<D>,
108    S: Differentiable,
109{
110    fn initial_state(&self) -> S {
111        self.0.initial_state()
112    }
113    fn nodal_forces_into(
114        &self,
115        nodal_coordinates: &NodalCoordinates<D>,
116        state_variables: &S,
117        nodal_forces: &mut NodalForcesSolid<D>,
118    ) -> Result<(), ElementModelError> {
119        self.0
120            .nodal_forces_into(nodal_coordinates, state_variables, nodal_forces)?;
121        self.1.nodal_forces_into(nodal_coordinates, nodal_forces)
122    }
123    fn nodal_stiffnesses_into(
124        &self,
125        nodal_coordinates: &NodalCoordinates<D>,
126        state_variables: &S,
127        nodal_stiffnesses: &mut NodalStiffnessesSolid<D>,
128    ) -> Result<(), ElementModelError> {
129        self.0
130            .nodal_stiffnesses_into(nodal_coordinates, state_variables, nodal_stiffnesses)?;
131        self.1
132            .nodal_stiffnesses_into(nodal_coordinates, nodal_stiffnesses)
133    }
134    fn state_variables_evolution(
135        &self,
136        nodal_coordinates: &NodalCoordinates<D>,
137        state_variables: &S,
138    ) -> Result<Derivative<S>, ElementModelError> {
139        self.0
140            .state_variables_evolution(nodal_coordinates, state_variables)
141    }
142}
143
144impl<B1, B2, S1, S2, const D: usize> ElasticViscoplasticElements<TensorTuple<S1, S2>, D>
145    for Blocks<B1, B2>
146where
147    B1: ElasticViscoplasticElements<S1, D>,
148    B2: ElasticViscoplasticElements<S2, D>,
149    S1: Differentiable + Tensor,
150    S2: Differentiable + Tensor,
151    Derivative<S1>: Tensor,
152    Derivative<S2>: Tensor,
153{
154    fn initial_state(&self) -> TensorTuple<S1, S2> {
155        (self.0.initial_state(), self.1.initial_state()).into()
156    }
157    fn nodal_forces_into(
158        &self,
159        nodal_coordinates: &NodalCoordinates<D>,
160        state_variables: &TensorTuple<S1, S2>,
161        nodal_forces: &mut NodalForcesSolid<D>,
162    ) -> Result<(), ElementModelError> {
163        self.0
164            .nodal_forces_into(nodal_coordinates, &state_variables.0, nodal_forces)?;
165        self.1
166            .nodal_forces_into(nodal_coordinates, &state_variables.1, nodal_forces)
167    }
168    fn nodal_stiffnesses_into(
169        &self,
170        nodal_coordinates: &NodalCoordinates<D>,
171        state_variables: &TensorTuple<S1, S2>,
172        nodal_stiffnesses: &mut NodalStiffnessesSolid<D>,
173    ) -> Result<(), ElementModelError> {
174        self.0
175            .nodal_stiffnesses_into(nodal_coordinates, &state_variables.0, nodal_stiffnesses)?;
176        self.1
177            .nodal_stiffnesses_into(nodal_coordinates, &state_variables.1, nodal_stiffnesses)
178    }
179    fn state_variables_evolution(
180        &self,
181        nodal_coordinates: &NodalCoordinates<D>,
182        state_variables: &TensorTuple<S1, S2>,
183    ) -> Result<Derivative<TensorTuple<S1, S2>>, ElementModelError> {
184        Ok((
185            self.0
186                .state_variables_evolution(nodal_coordinates, &state_variables.0)?,
187            self.1
188                .state_variables_evolution(nodal_coordinates, &state_variables.1)?,
189        )
190            .into())
191    }
192}
193
194/// The RKMK-DAE stage-equilibrium step, mirroring the constitutive-level
195/// `RootRkmkDae`: `F` (here, nodal equilibrium) is resolved at every RK stage
196/// abscissa from every Gauss point's stage-consistent plastic state, rather
197/// than frozen across the load-step window — so the coupling is the tableau's
198/// own order instead of first order, at the cost of `Tab::STAGES` equilibrium
199/// solves per window instead of one. One blanket impl over any [`Model`] whose
200/// blocks are [`ElasticViscoplasticDaeElements`] — a single block, nested
201/// [`Blocks`] to any depth, or an [`ElasticViscoplasticAndElastic`] pairing.
202pub trait RootRkmkDae<const D: usize, Y = Quantity> {
203    /// The model's plastic-state history type.
204    type History;
205    /// Solve under an applied load, resolving nodal equilibrium at every RK
206    /// stage of every load-step window while every Gauss point's plastic
207    /// state advances on its group.
208    fn root_rkmk_dae<Tab: ButcherTableau>(
209        &self,
210        solver: impl FirstOrderRootFinding<
211            NodalForcesSolid<D>,
212            NodalStiffnessesSolid<D>,
213            NodalCoordinates<D>,
214        >,
215        time: &[Quantity<Time>],
216        bcs: ElasticViscoplasticBCs,
217    ) -> Result<(Times, NodalCoordinatesHistory<D>, Self::History), IntegrationError>;
218    /// As [`Self::root_rkmk_dae`], but the whole span is stepped under
219    /// embedded (`Tab::D`) error control rather than on the supplied load
220    /// grid.
221    ///
222    /// Two times in `time` give only the span, and the controller's own
223    /// accepted steps are reported. More than two are requested report
224    /// times — the convention of the flat DAE loop — and every Gauss point's
225    /// plastic state is served at each from the geodesic `HermiteSegment` of
226    /// the accepted step containing it, so nodal equilibrium is then
227    /// re-solved there.
228    fn root_rkmk_dae_adaptive<Tab: EmbeddedTableau>(
229        &self,
230        solver: impl FirstOrderRootFinding<
231            NodalForcesSolid<D>,
232            NodalStiffnessesSolid<D>,
233            NodalCoordinates<D>,
234        >,
235        time: &[Quantity<Time>],
236        bcs: ElasticViscoplasticBCs,
237        abs_tol: Scalar,
238        rel_tol: Scalar,
239    ) -> Result<(Times, NodalCoordinatesHistory<D>, Self::History), IntegrationError>;
240}
241
242/// Per-topology machinery behind [`RootRkmkDae`]: the whole-mesh
243/// [`Integrable`] every Gauss point's plastic state lives on (a
244/// per-Gauss-point [`List`](crate::math::integrate::List) for one block,
245/// a [`Product`] of those for [`Blocks`]), how to flatten/unflatten between
246/// it and the block's native per-element [`State`](Self::State), and how to
247/// evaluate every Gauss point's plastic rate at a stage's nodal coordinates.
248/// Composes recursively over nested [`Blocks`] and forwards over
249/// [`ElasticViscoplasticAndElastic`], the same way
250/// [`ElasticViscoplasticElements`] does.
251pub trait ElasticViscoplasticDaeElements<Y, const D: usize>
252where
253    Self: Elements,
254{
255    /// The whole-mesh field this topology's plastic state lives on.
256    type Field: Integrable<Increment: Differentiable<Time>>;
257    /// The per-element grouped state — a per-Gauss-point list for one block,
258    /// a [`TensorTuple`] of those for [`Blocks`].
259    type State: Clone + Differentiable + Tensor;
260    /// Time history of [`Self::State`].
261    type History: TensorVec<Item = Self::State>;
262    /// Flattens [`Self::State`] into [`Self::Field`]'s `Point`.
263    fn flatten(state: &Self::State) -> <Self::Field as Integrable>::Point;
264    /// The inverse of [`Self::flatten`].
265    fn unflatten(flat: &<Self::Field as Integrable>::Point) -> Self::State;
266    /// Evaluates every Gauss point's plastic rate at `nodal_coordinates`.
267    fn dae_rate(
268        &self,
269        t: Quantity<Time>,
270        nodal_coordinates: &NodalCoordinates<D>,
271        flat: &<Self::Field as Integrable>::Point,
272    ) -> Result<Derivative<<Self::Field as Integrable>::Increment, Time>, ElementModelError>;
273}
274
275impl<B1, B2, Y, const D: usize> ElasticViscoplasticDaeElements<Y, D> for Blocks<B1, B2>
276where
277    B1: ElasticViscoplasticDaeElements<Y, D>,
278    B2: ElasticViscoplasticDaeElements<Y, D>,
279    TensorTuple<<B1::Field as Integrable>::Point, <B2::Field as Integrable>::Point>: Tensor,
280    TensorTuple<<B1::Field as Integrable>::Increment, <B2::Field as Integrable>::Increment>: Tensor,
281    TensorTuple<B1::State, B2::State>: Clone + Differentiable + Tensor,
282    Derivative<B1::State>: Tensor,
283    Derivative<B2::State>: Tensor,
284    TensorTupleVec<B1::State, B2::State>: TensorVec<Item = TensorTuple<B1::State, B2::State>>,
285{
286    type Field = Product<B1::Field, B2::Field>;
287    type State = TensorTuple<B1::State, B2::State>;
288    type History = TensorTupleVec<B1::State, B2::State>;
289    fn flatten(state: &Self::State) -> <Self::Field as Integrable>::Point {
290        TensorTuple(B1::flatten(&state.0), B2::flatten(&state.1))
291    }
292    fn unflatten(flat: &<Self::Field as Integrable>::Point) -> Self::State {
293        TensorTuple(B1::unflatten(&flat.0), B2::unflatten(&flat.1))
294    }
295    fn dae_rate(
296        &self,
297        t: Quantity<Time>,
298        nodal_coordinates: &NodalCoordinates<D>,
299        flat: &<Self::Field as Integrable>::Point,
300    ) -> Result<Derivative<<Self::Field as Integrable>::Increment, Time>, ElementModelError> {
301        Ok((
302            self.0.dae_rate(t, nodal_coordinates, &flat.0)?,
303            self.1.dae_rate(t, nodal_coordinates, &flat.1)?,
304        )
305            .into())
306    }
307}
308
309impl<B1, B2, Y, const D: usize> ElasticViscoplasticDaeElements<Y, D>
310    for ElasticViscoplasticAndElastic<B1, B2>
311where
312    B1: ElasticViscoplasticDaeElements<Y, D>,
313    B2: ElasticElements<D>,
314{
315    type Field = B1::Field;
316    type State = B1::State;
317    type History = B1::History;
318    fn flatten(state: &Self::State) -> <Self::Field as Integrable>::Point {
319        B1::flatten(state)
320    }
321    fn unflatten(flat: &<Self::Field as Integrable>::Point) -> Self::State {
322        B1::unflatten(flat)
323    }
324    fn dae_rate(
325        &self,
326        t: Quantity<Time>,
327        nodal_coordinates: &NodalCoordinates<D>,
328        flat: &<Self::Field as Integrable>::Point,
329    ) -> Result<Derivative<<Self::Field as Integrable>::Increment, Time>, ElementModelError> {
330        self.0.dae_rate(t, nodal_coordinates, flat)
331    }
332}
333
334impl<B, Y> RootRkmkDae<3, Y> for Model<B, 3>
335where
336    B: ElasticViscoplasticDaeElements<Y, 3> + ElasticViscoplasticElements<B::State, 3>,
337    <B::Field as Integrable>::Point: Clone,
338    <B::Field as Integrable>::Increment: Clone + Differentiable<Time>,
339    for<'a> &'a Derivative<<B::Field as Integrable>::Increment, Time>:
340        Mul<Quantity<Time>, Output = <B::Field as Integrable>::Increment>,
341    Derivative<<B::Field as Integrable>::Increment, Time>:
342        Mul<Quantity<Time>, Output = <B::Field as Integrable>::Increment>,
343    B::State: Clone,
344    B::History: TensorVec<Item = B::State>,
345{
346    type History = B::History;
347    #[allow(clippy::type_complexity)]
348    fn root_rkmk_dae<Tab: ButcherTableau>(
349        &self,
350        solver: impl FirstOrderRootFinding<
351            NodalForcesSolid<3>,
352            NodalStiffnessesSolid<3>,
353            NodalCoordinates<3>,
354        >,
355        time: &[Quantity<Time>],
356        bcs: ElasticViscoplasticBCs,
357    ) -> Result<(Times, NodalCoordinatesHistory<3>, Self::History), IntegrationError> {
358        let blocks = self.blocks();
359        let function = |_: Quantity<Time>,
360                        state: &<B::Field as Integrable>::Point,
361                        nodal_coordinates: &NodalCoordinates<3>|
362         -> Result<NodalForcesSolid<3>, String> {
363            Ok(blocks.nodal_forces(nodal_coordinates, &B::unflatten(state))?)
364        };
365        let jacobian = |_: Quantity<Time>,
366                        state: &<B::Field as Integrable>::Point,
367                        nodal_coordinates: &NodalCoordinates<3>|
368         -> Result<NodalStiffnessesSolid<3>, String> {
369            Ok(blocks.nodal_stiffnesses(nodal_coordinates, &B::unflatten(state))?)
370        };
371        let rate =
372            |t: Quantity<Time>,
373             state: &<B::Field as Integrable>::Point,
374             nodal_coordinates: &NodalCoordinates<3>|
375             -> Result<Derivative<<B::Field as Integrable>::Increment, Time>, String> {
376                Ok(blocks.dae_rate(t, nodal_coordinates, state)?)
377            };
378        let equality_constraint = bcs;
379        let mut state = B::flatten(&ElasticViscoplasticElements::initial_state(blocks));
380        let guess: NodalCoordinates<3> = self.coordinates().clone().into();
381        let mut nodal_coordinates = solver
382            .root(
383                |x: &NodalCoordinates<3>| function(time[0], &state, x),
384                |x: &NodalCoordinates<3>| jacobian(time[0], &state, x),
385                guess,
386                equality_constraint(time[0]),
387                None,
388            )
389            .map_err(|error| IntegrationError::from(format!("{error:?}")))?;
390        let mut times = Times::new();
391        let mut nodal_coordinates_history = NodalCoordinatesHistory::new();
392        let mut state_variables_history = Self::History::new();
393        let mut scratch = Vec::new();
394        let mut carry = None;
395        times.push(time[0]);
396        nodal_coordinates_history.push(nodal_coordinates.clone());
397        state_variables_history.push(B::unflatten(&state));
398        for step in time.windows(2) {
399            let advanced = rkmk_dae_step_first_order_root::<
400                B::Field,
401                Tab,
402                NodalForcesSolid<3>,
403                NodalStiffnessesSolid<3>,
404                NodalCoordinates<3>,
405                Time,
406            >(
407                &mut |t, state, nodal_coordinates| rate(t, state, nodal_coordinates),
408                function,
409                jacobian,
410                &solver,
411                &state,
412                &nodal_coordinates,
413                step[0],
414                step[1] - step[0],
415                &mut scratch,
416                carry.as_ref(),
417                equality_constraint,
418            )
419            .map_err(|error| IntegrationError::from(format!("{error:?}")))?;
420            state = advanced.0;
421            nodal_coordinates = advanced.1;
422            carry = advanced.2;
423            times.push(step[1]);
424            nodal_coordinates_history.push(nodal_coordinates.clone());
425            state_variables_history.push(B::unflatten(&state));
426        }
427        Ok((times, nodal_coordinates_history, state_variables_history))
428    }
429    #[allow(clippy::type_complexity)]
430    fn root_rkmk_dae_adaptive<Tab: EmbeddedTableau>(
431        &self,
432        solver: impl FirstOrderRootFinding<
433            NodalForcesSolid<3>,
434            NodalStiffnessesSolid<3>,
435            NodalCoordinates<3>,
436        >,
437        time: &[Quantity<Time>],
438        bcs: ElasticViscoplasticBCs,
439        abs_tol: Scalar,
440        rel_tol: Scalar,
441    ) -> Result<(Times, NodalCoordinatesHistory<3>, Self::History), IntegrationError> {
442        let blocks = self.blocks();
443        let function = |_: Quantity<Time>,
444                        state: &<B::Field as Integrable>::Point,
445                        nodal_coordinates: &NodalCoordinates<3>|
446         -> Result<NodalForcesSolid<3>, String> {
447            Ok(blocks.nodal_forces(nodal_coordinates, &B::unflatten(state))?)
448        };
449        let jacobian = |_: Quantity<Time>,
450                        state: &<B::Field as Integrable>::Point,
451                        nodal_coordinates: &NodalCoordinates<3>|
452         -> Result<NodalStiffnessesSolid<3>, String> {
453            Ok(blocks.nodal_stiffnesses(nodal_coordinates, &B::unflatten(state))?)
454        };
455        let rate =
456            |t: Quantity<Time>,
457             state: &<B::Field as Integrable>::Point,
458             nodal_coordinates: &NodalCoordinates<3>|
459             -> Result<Derivative<<B::Field as Integrable>::Increment, Time>, String> {
460                Ok(blocks.dae_rate(t, nodal_coordinates, state)?)
461            };
462        let equality_constraint = bcs;
463        let state = B::flatten(&ElasticViscoplasticElements::initial_state(blocks));
464        let guess: NodalCoordinates<3> = self.coordinates().clone().into();
465        let nodal_coordinates = solver
466            .root(
467                |x: &NodalCoordinates<3>| function(time[0], &state, x),
468                |x: &NodalCoordinates<3>| jacobian(time[0], &state, x),
469                guess,
470                equality_constraint(time[0]),
471                None,
472            )
473            .map_err(|error| IntegrationError::from(format!("{error:?}")))?;
474        let (times, state_points_history, nodal_coordinates_history) =
475            integrate_rkmk_dae_adaptive_first_order_root::<
476                B::Field,
477                Tab,
478                NodalForcesSolid<3>,
479                NodalStiffnessesSolid<3>,
480                NodalCoordinates<3>,
481                TensorVector<<B::Field as Integrable>::Point>,
482                NodalCoordinatesHistory<3>,
483                Time,
484            >(
485                |t, state, nodal_coordinates| rate(t, state, nodal_coordinates),
486                function,
487                jacobian,
488                &solver,
489                time,
490                (state, nodal_coordinates),
491                abs_tol,
492                rel_tol,
493                equality_constraint,
494            )
495            .map_err(|error| IntegrationError::from(format!("{error:?}")))?;
496        let state_variables_history = state_points_history.iter().map(B::unflatten).collect();
497        Ok((times, nodal_coordinates_history, state_variables_history))
498    }
499}
conspire/domain/solid/elastic_viscoplastic/mod.rs

conspire/domain/solid/elastic_viscoplastic/
mod.rs
