conspire/domain/solid/hyperelastic_viscoplastic/
mod.rs

1use crate::{
2    domain::{
3        Blocks, ElasticViscoplasticAndElastic, ElementModel, ElementModelError, Model,
4        NodalCoordinates, NodalCoordinatesHistory,
5        block::{element::Elements, finalize_node_neighbors, solver_from_neighbors},
6        solid::{
7            NodalForcesSolid, NodalStiffnessesSolid,
8            elastic_viscoplastic::{
9                ElasticViscoplasticBCs, ElasticViscoplasticDaeElements, ElasticViscoplasticElements,
10            },
11            hyperelastic::HyperelasticElements,
12        },
13    },
14    math::{
15        Derivative, Differentiable, Quantity, Scalar, Tensor, TensorTuple, TensorVec, TensorVector,
16        integrate::{
17            ButcherTableau, EmbeddedTableau, Integrable, IntegrationError,
18            integrate_rkmk_dae_adaptive_second_order_minimize, rkmk_dae_step_second_order_minimize,
19        },
20        optimize::SecondOrderOptimization,
21    },
22    mechanics::Times,
23    units::{Energy, Time},
24};
25use std::ops::Mul;
26
27pub trait HyperelasticViscoplasticElements<S, const D: usize>
28where
29    Self: ElasticViscoplasticElements<S, D>,
30    S: Differentiable,
31{
32    fn helmholtz_free_energy(
33        &self,
34        nodal_coordinates: &NodalCoordinates<D>,
35        state_variables: &S,
36    ) -> Result<Quantity<Energy>, ElementModelError>;
37}
38
39impl<B, S, const D: usize> HyperelasticViscoplasticElements<S, D> for Model<B, D>
40where
41    B: HyperelasticViscoplasticElements<S, D>,
42    S: Differentiable,
43{
44    fn helmholtz_free_energy(
45        &self,
46        nodal_coordinates: &NodalCoordinates<D>,
47        state_variables: &S,
48    ) -> Result<Quantity<Energy>, ElementModelError> {
49        self.blocks
50            .helmholtz_free_energy(nodal_coordinates, state_variables)
51    }
52}
53
54impl<B1, B2, S, const D: usize> HyperelasticViscoplasticElements<S, D>
55    for ElasticViscoplasticAndElastic<B1, B2>
56where
57    B1: HyperelasticViscoplasticElements<S, D>,
58    B2: HyperelasticElements<D>,
59    S: Differentiable,
60{
61    fn helmholtz_free_energy(
62        &self,
63        nodal_coordinates: &NodalCoordinates<D>,
64        state_variables: &S,
65    ) -> Result<Quantity<Energy>, ElementModelError> {
66        Ok(self
67            .0
68            .helmholtz_free_energy(nodal_coordinates, state_variables)?
69            + self.1.helmholtz_free_energy(nodal_coordinates)?)
70    }
71}
72
73impl<B1, B2, S1, S2, const D: usize> HyperelasticViscoplasticElements<TensorTuple<S1, S2>, D>
74    for Blocks<B1, B2>
75where
76    B1: HyperelasticViscoplasticElements<S1, D>,
77    B2: HyperelasticViscoplasticElements<S2, D>,
78    S1: Differentiable + Tensor,
79    S2: Differentiable + Tensor,
80    Derivative<S1>: Tensor,
81    Derivative<S2>: Tensor,
82{
83    fn helmholtz_free_energy(
84        &self,
85        nodal_coordinates: &NodalCoordinates<D>,
86        state_variables: &TensorTuple<S1, S2>,
87    ) -> Result<Quantity<Energy>, ElementModelError> {
88        Ok(self
89            .0
90            .helmholtz_free_energy(nodal_coordinates, &state_variables.0)?
91            + self
92                .1
93                .helmholtz_free_energy(nodal_coordinates, &state_variables.1)?)
94    }
95}
96
97/// The minimize-based sibling of
98/// [`crate::domain::solid::elastic_viscoplastic::RootRkmkDae`]: `F` (here,
99/// nodal equilibrium) is resolved by potential minimization at every RK stage
100/// abscissa from every Gauss point's stage-consistent plastic state, for
101/// models whose equilibrium is naturally posed that way rather than as a
102/// stress residual. One blanket impl over any [`Model`] whose blocks are
103/// [`ElasticViscoplasticDaeElements`] (the topology/flatten/unflatten/rate
104/// machinery, shared with the root-finding sibling since the field shape
105/// doesn't care how a stage is solved) and [`HyperelasticViscoplasticElements`]
106/// (for the potential itself) — a single [`Block`](crate::fem::block::Block),
107/// nested [`Blocks`] to any depth, or an [`ElasticViscoplasticAndElastic`]
108/// pairing.
109pub trait RootRkmkDaeMinimize<const D: usize, Y = Quantity> {
110    /// The model's plastic-state history type.
111    type History;
112    /// Solve under an applied load, resolving nodal equilibrium by potential
113    /// minimization at every RK stage of every load-step window while every
114    /// Gauss point's plastic state advances on its group.
115    fn root_rkmk_dae_minimize<Tab: ButcherTableau>(
116        &self,
117        solver: impl SecondOrderOptimization<
118            Quantity<Energy>,
119            NodalForcesSolid<D>,
120            NodalStiffnessesSolid<D>,
121            NodalCoordinates<D>,
122        >,
123        time: &[Quantity<Time>],
124        bcs: ElasticViscoplasticBCs,
125    ) -> Result<(Times, NodalCoordinatesHistory<D>, Self::History), IntegrationError>;
126    /// As [`Self::root_rkmk_dae_minimize`], but the whole span is stepped
127    /// under embedded (`Tab::D`) error control rather than on the supplied
128    /// load grid.
129    ///
130    /// Two times in `time` give only the span, and the controller's own
131    /// accepted steps are reported. More than two are requested report
132    /// times — the convention of the flat DAE loop — and every Gauss point's
133    /// plastic state is served at each from the geodesic `HermiteSegment` of
134    /// the accepted step containing it, so nodal equilibrium is then
135    /// re-solved (by minimization) there.
136    fn root_rkmk_dae_adaptive_minimize<Tab: EmbeddedTableau>(
137        &self,
138        solver: impl SecondOrderOptimization<
139            Quantity<Energy>,
140            NodalForcesSolid<D>,
141            NodalStiffnessesSolid<D>,
142            NodalCoordinates<D>,
143        >,
144        time: &[Quantity<Time>],
145        bcs: ElasticViscoplasticBCs,
146        abs_tol: Scalar,
147        rel_tol: Scalar,
148    ) -> Result<(Times, NodalCoordinatesHistory<D>, Self::History), IntegrationError>;
149}
150
151impl<B, Y> RootRkmkDaeMinimize<3, Y> for Model<B, 3>
152where
153    B: ElasticViscoplasticDaeElements<Y, 3> + HyperelasticViscoplasticElements<B::State, 3>,
154    <B::Field as Integrable>::Point: Clone,
155    <B::Field as Integrable>::Increment: Clone + Differentiable<Time>,
156    for<'a> &'a Derivative<<B::Field as Integrable>::Increment, Time>:
157        Mul<Quantity<Time>, Output = <B::Field as Integrable>::Increment>,
158    Derivative<<B::Field as Integrable>::Increment, Time>:
159        Mul<Quantity<Time>, Output = <B::Field as Integrable>::Increment>,
160    B::State: Clone,
161    B::History: TensorVec<Item = B::State>,
162{
163    type History = B::History;
164    #[allow(clippy::type_complexity)]
165    fn root_rkmk_dae_minimize<Tab: ButcherTableau>(
166        &self,
167        solver: impl SecondOrderOptimization<
168            Quantity<Energy>,
169            NodalForcesSolid<3>,
170            NodalStiffnessesSolid<3>,
171            NodalCoordinates<3>,
172        >,
173        time: &[Quantity<Time>],
174        bcs: ElasticViscoplasticBCs,
175    ) -> Result<(Times, NodalCoordinatesHistory<3>, Self::History), IntegrationError> {
176        let blocks = self.blocks();
177        let mut neighbors = vec![Vec::new(); self.coordinates().len()];
178        self.node_neighbors(&mut neighbors);
179        finalize_node_neighbors(&mut neighbors);
180        let sparse = solver_from_neighbors(&neighbors, &bcs(time[0]), 3, true);
181        let function = |_: Quantity<Time>,
182                        state: &<B::Field as Integrable>::Point,
183                        nodal_coordinates: &NodalCoordinates<3>|
184         -> Result<Quantity<Energy>, String> {
185            Ok(blocks.helmholtz_free_energy(nodal_coordinates, &B::unflatten(state))?)
186        };
187        let jacobian = |_: Quantity<Time>,
188                        state: &<B::Field as Integrable>::Point,
189                        nodal_coordinates: &NodalCoordinates<3>|
190         -> Result<NodalForcesSolid<3>, String> {
191            Ok(blocks.nodal_forces(nodal_coordinates, &B::unflatten(state))?)
192        };
193        let hessian = |_: Quantity<Time>,
194                       state: &<B::Field as Integrable>::Point,
195                       nodal_coordinates: &NodalCoordinates<3>|
196         -> Result<NodalStiffnessesSolid<3>, String> {
197            Ok(blocks.nodal_stiffnesses(nodal_coordinates, &B::unflatten(state))?)
198        };
199        let rate =
200            |t: Quantity<Time>,
201             state: &<B::Field as Integrable>::Point,
202             nodal_coordinates: &NodalCoordinates<3>|
203             -> Result<Derivative<<B::Field as Integrable>::Increment, Time>, String> {
204                Ok(blocks.dae_rate(t, nodal_coordinates, state)?)
205            };
206        let equality_constraint = bcs;
207        let mut state = B::flatten(&ElasticViscoplasticElements::initial_state(blocks));
208        let guess: NodalCoordinates<3> = self.coordinates().clone().into();
209        let mut nodal_coordinates = solver
210            .minimize(
211                |x: &NodalCoordinates<3>| function(time[0], &state, x),
212                |x: &NodalCoordinates<3>| jacobian(time[0], &state, x),
213                |x: &NodalCoordinates<3>| hessian(time[0], &state, x),
214                guess,
215                equality_constraint(time[0]),
216                Some(sparse.clone()),
217            )
218            .map_err(|error| IntegrationError::from(format!("{error:?}")))?;
219        let mut times = Times::new();
220        let mut nodal_coordinates_history = NodalCoordinatesHistory::new();
221        let mut state_variables_history = Self::History::new();
222        let mut scratch = Vec::new();
223        let mut carry = None;
224        times.push(time[0]);
225        nodal_coordinates_history.push(nodal_coordinates.clone());
226        state_variables_history.push(B::unflatten(&state));
227        for step in time.windows(2) {
228            let advanced = rkmk_dae_step_second_order_minimize::<
229                B::Field,
230                Tab,
231                Quantity<Energy>,
232                NodalForcesSolid<3>,
233                NodalStiffnessesSolid<3>,
234                NodalCoordinates<3>,
235                Time,
236            >(
237                &mut |t, state, nodal_coordinates| rate(t, state, nodal_coordinates),
238                function,
239                jacobian,
240                hessian,
241                &solver,
242                &state,
243                &nodal_coordinates,
244                step[0],
245                step[1] - step[0],
246                &mut scratch,
247                carry.as_ref(),
248                equality_constraint,
249                Some(sparse.clone()),
250            )
251            .map_err(|error| IntegrationError::from(format!("{error:?}")))?;
252            state = advanced.0;
253            nodal_coordinates = advanced.1;
254            carry = advanced.2;
255            times.push(step[1]);
256            nodal_coordinates_history.push(nodal_coordinates.clone());
257            state_variables_history.push(B::unflatten(&state));
258        }
259        Ok((times, nodal_coordinates_history, state_variables_history))
260    }
261    #[allow(clippy::type_complexity)]
262    fn root_rkmk_dae_adaptive_minimize<Tab: EmbeddedTableau>(
263        &self,
264        solver: impl SecondOrderOptimization<
265            Quantity<Energy>,
266            NodalForcesSolid<3>,
267            NodalStiffnessesSolid<3>,
268            NodalCoordinates<3>,
269        >,
270        time: &[Quantity<Time>],
271        bcs: ElasticViscoplasticBCs,
272        abs_tol: Scalar,
273        rel_tol: Scalar,
274    ) -> Result<(Times, NodalCoordinatesHistory<3>, Self::History), IntegrationError> {
275        let blocks = self.blocks();
276        let mut neighbors = vec![Vec::new(); self.coordinates().len()];
277        self.node_neighbors(&mut neighbors);
278        finalize_node_neighbors(&mut neighbors);
279        let sparse = solver_from_neighbors(&neighbors, &bcs(time[0]), 3, true);
280        let function = |_: Quantity<Time>,
281                        state: &<B::Field as Integrable>::Point,
282                        nodal_coordinates: &NodalCoordinates<3>|
283         -> Result<Quantity<Energy>, String> {
284            Ok(blocks.helmholtz_free_energy(nodal_coordinates, &B::unflatten(state))?)
285        };
286        let jacobian = |_: Quantity<Time>,
287                        state: &<B::Field as Integrable>::Point,
288                        nodal_coordinates: &NodalCoordinates<3>|
289         -> Result<NodalForcesSolid<3>, String> {
290            Ok(blocks.nodal_forces(nodal_coordinates, &B::unflatten(state))?)
291        };
292        let hessian = |_: Quantity<Time>,
293                       state: &<B::Field as Integrable>::Point,
294                       nodal_coordinates: &NodalCoordinates<3>|
295         -> Result<NodalStiffnessesSolid<3>, String> {
296            Ok(blocks.nodal_stiffnesses(nodal_coordinates, &B::unflatten(state))?)
297        };
298        let rate =
299            |t: Quantity<Time>,
300             state: &<B::Field as Integrable>::Point,
301             nodal_coordinates: &NodalCoordinates<3>|
302             -> Result<Derivative<<B::Field as Integrable>::Increment, Time>, String> {
303                Ok(blocks.dae_rate(t, nodal_coordinates, state)?)
304            };
305        let equality_constraint = bcs;
306        let state = B::flatten(&ElasticViscoplasticElements::initial_state(blocks));
307        let guess: NodalCoordinates<3> = self.coordinates().clone().into();
308        let nodal_coordinates = solver
309            .minimize(
310                |x: &NodalCoordinates<3>| function(time[0], &state, x),
311                |x: &NodalCoordinates<3>| jacobian(time[0], &state, x),
312                |x: &NodalCoordinates<3>| hessian(time[0], &state, x),
313                guess,
314                equality_constraint(time[0]),
315                Some(sparse.clone()),
316            )
317            .map_err(|error| IntegrationError::from(format!("{error:?}")))?;
318        let (times, state_points_history, nodal_coordinates_history) =
319            integrate_rkmk_dae_adaptive_second_order_minimize::<
320                B::Field,
321                Tab,
322                Quantity<Energy>,
323                NodalForcesSolid<3>,
324                NodalStiffnessesSolid<3>,
325                NodalCoordinates<3>,
326                TensorVector<<B::Field as Integrable>::Point>,
327                NodalCoordinatesHistory<3>,
328                Time,
329            >(
330                |t, state, nodal_coordinates| rate(t, state, nodal_coordinates),
331                function,
332                jacobian,
333                hessian,
334                &solver,
335                time,
336                (state, nodal_coordinates),
337                abs_tol,
338                rel_tol,
339                equality_constraint,
340                Some(sparse),
341            )
342            .map_err(|error| IntegrationError::from(format!("{error:?}")))?;
343        let state_variables_history = state_points_history.iter().map(B::unflatten).collect();
344        Ok((times, nodal_coordinates_history, state_variables_history))
345    }
346}
conspire/domain/solid/hyperelastic_viscoplastic/mod.rs

conspire/domain/solid/hyperelastic_viscoplastic/
mod.rs
