conspire/constitutive/solid/hyperelastic_viscoplastic/
mod.rs

1//! Hyperelastic-viscoplastic solid constitutive models.
2//!
3//! ---
4//!
5#![doc = include_str!("doc.md")]
6
7#[cfg(feature = "doc")]
8pub mod doc;
9
10mod canonical;
11
12use crate::{
13    constitutive::{
14        ConstitutiveError,
15        fluid::viscoplastic::{
16            ViscoplasticEvolutionHistory, ViscoplasticStateVariables,
17            ViscoplasticStateVariablesHistory,
18        },
19        solid::{
20            elastic_plastic::bcs,
21            elastic_viscoplastic::{AppliedLoad, ElasticViscoplastic},
22        },
23    },
24    math::{
25        Derivative, Differentiable, Quantity, Scalar, Tensor, TensorArray, TensorVec, Vector,
26        integrate::{
27            ButcherTableau, EmbeddedTableau, EvolvedIncrement, ExplicitDaeFirstOrderMinimize,
28            ExplicitDaeSecondOrderMinimize, Integrable, StateEvolution,
29            integrate_rkmk_dae_adaptive_second_order_minimize, rkmk_dae_step_second_order_minimize,
30        },
31        optimize::{EqualityConstraint, FirstOrderOptimization, SecondOrderOptimization},
32    },
33    mechanics::{
34        DeformationGradient, DeformationGradientPlastic, DeformationGradients,
35        FirstPiolaKirchhoffStress, FirstPiolaKirchhoffTangentStiffness, Times,
36    },
37    units::{EnergyDensity, Time},
38};
39use std::ops::Mul;
40
41/// Required methods for hyperelastic-viscoplastic solid constitutive models.
42pub trait HyperelasticViscoplastic<Y>
43where
44    Self: ElasticViscoplastic<Y>,
45    Y: Differentiable + Tensor,
46{
47    /// Calculates and returns the Helmholtz free energy density.
48    ///
49    /// ```math
50    /// a = a(\mathbf{F}_\mathrm{e})
51    /// ```
52    fn helmholtz_free_energy_density(
53        &self,
54        deformation_gradient: &DeformationGradient,
55        deformation_gradient_p: &DeformationGradientPlastic,
56    ) -> Result<Quantity<EnergyDensity>, ConstitutiveError>;
57}
58
59/// First-order minimization methods for hyperelastic-viscoplastic solid constitutive models.
60pub trait FirstOrderMinimize<Y>
61where
62    Y: Differentiable + Tensor,
63{
64    /// Solve for the unknown components of the deformation gradients under an applied load.
65    ///
66    /// ```math
67    /// \Pi(\mathbf{F},\mathbf{F}_\mathrm{p},\boldsymbol{\lambda}) = a(\mathbf{F},\mathbf{F}_\mathrm{p}) - \boldsymbol{\lambda}:(\mathbf{F} - \mathbf{F}_0) - \mathbf{P}_0:\mathbf{F}
68    /// ```
69    fn minimize(
70        &self,
71        applied_load: AppliedLoad,
72        integrator: impl ExplicitDaeFirstOrderMinimize<
73            Quantity<EnergyDensity>,
74            FirstPiolaKirchhoffStress,
75            ViscoplasticStateVariables<Y>,
76            DeformationGradient,
77            ViscoplasticStateVariablesHistory<Y>,
78            DeformationGradients,
79            ViscoplasticEvolutionHistory<Y>,
80        >,
81        solver: impl FirstOrderOptimization<
82            Quantity<EnergyDensity>,
83            FirstPiolaKirchhoffStress,
84            DeformationGradient,
85        >,
86    ) -> Result<
87        (
88            Times,
89            DeformationGradients,
90            ViscoplasticStateVariablesHistory<Y>,
91        ),
92        ConstitutiveError,
93    >;
94}
95
96/// Second-order minimization methods for hyperelastic-viscoplastic solid constitutive models.
97pub trait SecondOrderMinimize<Y>
98where
99    Y: Differentiable + Tensor,
100{
101    /// Solve for the unknown components of the deformation gradients under an applied load.
102    ///
103    /// ```math
104    /// \Pi(\mathbf{F},\mathbf{F}_\mathrm{p},\boldsymbol{\lambda}) = a(\mathbf{F},\mathbf{F}_\mathrm{p}) - \boldsymbol{\lambda}:(\mathbf{F} - \mathbf{F}_0) - \mathbf{P}_0:\mathbf{F}
105    /// ```
106    fn minimize(
107        &self,
108        applied_load: AppliedLoad,
109        integrator: impl ExplicitDaeSecondOrderMinimize<
110            Quantity<EnergyDensity>,
111            FirstPiolaKirchhoffStress,
112            FirstPiolaKirchhoffTangentStiffness,
113            ViscoplasticStateVariables<Y>,
114            DeformationGradient,
115            ViscoplasticStateVariablesHistory<Y>,
116            DeformationGradients,
117            ViscoplasticEvolutionHistory<Y>,
118        >,
119        solver: impl SecondOrderOptimization<
120            Quantity<EnergyDensity>,
121            FirstPiolaKirchhoffStress,
122            FirstPiolaKirchhoffTangentStiffness,
123            DeformationGradient,
124        >,
125    ) -> Result<
126        (
127            Times,
128            DeformationGradients,
129            ViscoplasticStateVariablesHistory<Y>,
130        ),
131        ConstitutiveError,
132    >;
133}
134
135impl<C, Y> FirstOrderMinimize<Y> for C
136where
137    C: HyperelasticViscoplastic<Y>,
138    Y: Differentiable + Tensor,
139{
140    fn minimize(
141        &self,
142        applied_load: AppliedLoad,
143        integrator: impl ExplicitDaeFirstOrderMinimize<
144            Quantity<EnergyDensity>,
145            FirstPiolaKirchhoffStress,
146            ViscoplasticStateVariables<Y>,
147            DeformationGradient,
148            ViscoplasticStateVariablesHistory<Y>,
149            DeformationGradients,
150            ViscoplasticEvolutionHistory<Y>,
151        >,
152        solver: impl FirstOrderOptimization<
153            Quantity<EnergyDensity>,
154            FirstPiolaKirchhoffStress,
155            DeformationGradient,
156        >,
157    ) -> Result<
158        (
159            Times,
160            DeformationGradients,
161            ViscoplasticStateVariablesHistory<Y>,
162        ),
163        ConstitutiveError,
164    > {
165        let (matrix, prescribed, time) = bcs(applied_load);
166        let mut vector = Vector::zero(matrix.len());
167        let (times, state_variables, _, deformation_gradients) = integrator
168            .integrate(
169                |_: Quantity<Time>,
170                 state_variables: &ViscoplasticStateVariables<Y>,
171                 deformation_gradient: &DeformationGradient| {
172                    Ok(self.state_variables_evolution(deformation_gradient, state_variables)?)
173                },
174                |_: Quantity<Time>,
175                 state_variables: &ViscoplasticStateVariables<Y>,
176                 deformation_gradient: &DeformationGradient| {
177                    let deformation_gradient_p = &state_variables.0;
178                    Ok(self.helmholtz_free_energy_density(
179                        deformation_gradient,
180                        deformation_gradient_p,
181                    )?)
182                },
183                |_: Quantity<Time>,
184                 state_variables: &ViscoplasticStateVariables<Y>,
185                 deformation_gradient: &DeformationGradient| {
186                    let deformation_gradient_p = &state_variables.0;
187                    Ok(self.first_piola_kirchhoff_stress(
188                        deformation_gradient,
189                        deformation_gradient_p,
190                    )?)
191                },
192                solver,
193                time,
194                (self.initial_state(), DeformationGradient::identity()),
195                |t: Quantity<Time>| {
196                    prescribed
197                        .iter()
198                        .for_each(|(index, function)| vector[*index] = function(t));
199                    EqualityConstraint::Linear(matrix.clone(), vector.clone())
200                },
201            )
202            .map_err(|error| ConstitutiveError::upstream(error, self))?;
203        Ok((times, deformation_gradients, state_variables))
204    }
205}
206
207impl<C, Y> SecondOrderMinimize<Y> for C
208where
209    C: HyperelasticViscoplastic<Y>,
210    Y: Differentiable + Tensor,
211{
212    fn minimize(
213        &self,
214        applied_load: AppliedLoad,
215        integrator: impl ExplicitDaeSecondOrderMinimize<
216            Quantity<EnergyDensity>,
217            FirstPiolaKirchhoffStress,
218            FirstPiolaKirchhoffTangentStiffness,
219            ViscoplasticStateVariables<Y>,
220            DeformationGradient,
221            ViscoplasticStateVariablesHistory<Y>,
222            DeformationGradients,
223            ViscoplasticEvolutionHistory<Y>,
224        >,
225        solver: impl SecondOrderOptimization<
226            Quantity<EnergyDensity>,
227            FirstPiolaKirchhoffStress,
228            FirstPiolaKirchhoffTangentStiffness,
229            DeformationGradient,
230        >,
231    ) -> Result<
232        (
233            Times,
234            DeformationGradients,
235            ViscoplasticStateVariablesHistory<Y>,
236        ),
237        ConstitutiveError,
238    > {
239        let (matrix, prescribed, time) = bcs(applied_load);
240        let mut vector = Vector::zero(matrix.len());
241        let (times, state_variables, _, deformation_gradients) = integrator
242            .integrate(
243                |_: Quantity<Time>,
244                 state_variables: &ViscoplasticStateVariables<Y>,
245                 deformation_gradient: &DeformationGradient| {
246                    Ok(self.state_variables_evolution(deformation_gradient, state_variables)?)
247                },
248                |_: Quantity<Time>,
249                 state_variables: &ViscoplasticStateVariables<Y>,
250                 deformation_gradient: &DeformationGradient| {
251                    let deformation_gradient_p = &state_variables.0;
252                    Ok(self.helmholtz_free_energy_density(
253                        deformation_gradient,
254                        deformation_gradient_p,
255                    )?)
256                },
257                |_: Quantity<Time>,
258                 state_variables: &ViscoplasticStateVariables<Y>,
259                 deformation_gradient: &DeformationGradient| {
260                    let deformation_gradient_p = &state_variables.0;
261                    Ok(self.first_piola_kirchhoff_stress(
262                        deformation_gradient,
263                        deformation_gradient_p,
264                    )?)
265                },
266                |_: Quantity<Time>,
267                 state_variables: &ViscoplasticStateVariables<Y>,
268                 deformation_gradient: &DeformationGradient| {
269                    let deformation_gradient_p = &state_variables.0;
270                    Ok(self.first_piola_kirchhoff_tangent_stiffness(
271                        deformation_gradient,
272                        deformation_gradient_p,
273                    )?)
274                },
275                solver,
276                time,
277                (self.initial_state(), DeformationGradient::identity()),
278                |t: Quantity<Time>| {
279                    prescribed
280                        .iter()
281                        .for_each(|(index, function)| vector[*index] = function(t));
282                    EqualityConstraint::Linear(matrix.clone(), vector.clone())
283                },
284                None,
285            )
286            .map_err(|error| ConstitutiveError::upstream(error, self))?;
287        Ok((times, deformation_gradients, state_variables))
288    }
289}
290
291/// RKMK-DAE return-map methods for hyperelastic-viscoplastic solid
292/// constitutive models — the minimize-based sibling of
293/// [`crate::constitutive::solid::elastic_viscoplastic::RootRkmkDae`], for
294/// models whose equilibrium is naturally posed as a potential minimization
295/// rather than a stress-residual root. `F_p` still advances on its group at
296/// the tableau's own order; only how `F` is resolved within a stage differs.
297/// Blanket over any [`HyperelasticViscoplastic`] model, same as
298/// [`SecondOrderMinimize`] itself.
299pub trait RootRkmkDaeMinimize<Y>
300where
301    Y: Differentiable + Tensor,
302{
303    /// `F` is re-solved by potential minimization at every stage abscissa of
304    /// the window while `F_p` advances on its group, generic over the
305    /// hardening variable `Y`.
306    fn root_rkmk_dae_minimize<Tab: ButcherTableau>(
307        &self,
308        applied_load: AppliedLoad,
309        solver: impl SecondOrderOptimization<
310            Quantity<EnergyDensity>,
311            FirstPiolaKirchhoffStress,
312            FirstPiolaKirchhoffTangentStiffness,
313            DeformationGradient,
314        >,
315    ) -> Result<
316        (
317            Times,
318            DeformationGradients,
319            ViscoplasticStateVariablesHistory<Y>,
320        ),
321        ConstitutiveError,
322    >;
323    /// As [`Self::root_rkmk_dae_minimize`], but the whole span is stepped
324    /// under embedded (`Tab::D`) error control rather than on the supplied
325    /// load grid.
326    fn root_rkmk_dae_adaptive_minimize<Tab: EmbeddedTableau>(
327        &self,
328        applied_load: AppliedLoad,
329        solver: impl SecondOrderOptimization<
330            Quantity<EnergyDensity>,
331            FirstPiolaKirchhoffStress,
332            FirstPiolaKirchhoffTangentStiffness,
333            DeformationGradient,
334        >,
335        abs_tol: Scalar,
336        rel_tol: Scalar,
337    ) -> Result<
338        (
339            Times,
340            DeformationGradients,
341            ViscoplasticStateVariablesHistory<Y>,
342        ),
343        ConstitutiveError,
344    >;
345}
346
347impl<C, Y> RootRkmkDaeMinimize<Y> for C
348where
349    C: HyperelasticViscoplastic<Y>
350        + StateEvolution<
351            Time,
352            Y,
353            Drive = DeformationGradient,
354            Field: Integrable<Point = ViscoplasticStateVariables<Y>>,
355        >,
356    Y: Differentiable + Tensor,
357    EvolvedIncrement<C, Time, Y>: Clone + Differentiable<Time>,
358    for<'a> &'a Derivative<EvolvedIncrement<C, Time, Y>, Time>:
359        Mul<Quantity<Time>, Output = EvolvedIncrement<C, Time, Y>>,
360{
361    #[allow(clippy::type_complexity)]
362    fn root_rkmk_dae_minimize<Tab: ButcherTableau>(
363        &self,
364        applied_load: AppliedLoad,
365        solver: impl SecondOrderOptimization<
366            Quantity<EnergyDensity>,
367            FirstPiolaKirchhoffStress,
368            FirstPiolaKirchhoffTangentStiffness,
369            DeformationGradient,
370        >,
371    ) -> Result<
372        (
373            Times,
374            DeformationGradients,
375            ViscoplasticStateVariablesHistory<Y>,
376        ),
377        ConstitutiveError,
378    > {
379        let (matrix, prescribed, time) = bcs(applied_load);
380        let mut state = <Self as StateEvolution<Time, Y>>::initial_state(self);
381        let mut scratch = Vec::new();
382        let equality_constraint = |t: Quantity<Time>| {
383            let mut vector = Vector::zero(matrix.len());
384            prescribed
385                .iter()
386                .for_each(|(index, function)| vector[*index] = function(t));
387            EqualityConstraint::Linear(matrix.clone(), vector)
388        };
389        let function = |_: Quantity<Time>,
390                        state: &ViscoplasticStateVariables<Y>,
391                        deformation_gradient: &DeformationGradient|
392         -> Result<Quantity<EnergyDensity>, String> {
393            Ok(self.helmholtz_free_energy_density(deformation_gradient, &state.0)?)
394        };
395        let jacobian = |_: Quantity<Time>,
396                        state: &ViscoplasticStateVariables<Y>,
397                        deformation_gradient: &DeformationGradient|
398         -> Result<FirstPiolaKirchhoffStress, String> {
399            Ok(self.first_piola_kirchhoff_stress(deformation_gradient, &state.0)?)
400        };
401        let hessian = |_: Quantity<Time>,
402                       state: &ViscoplasticStateVariables<Y>,
403                       deformation_gradient: &DeformationGradient|
404         -> Result<FirstPiolaKirchhoffTangentStiffness, String> {
405            Ok(self.first_piola_kirchhoff_tangent_stiffness(deformation_gradient, &state.0)?)
406        };
407        let mut deformation_gradient = solver
408            .minimize(
409                |deformation_gradient: &DeformationGradient| {
410                    function(time[0], &state, deformation_gradient)
411                },
412                |deformation_gradient: &DeformationGradient| {
413                    jacobian(time[0], &state, deformation_gradient)
414                },
415                |deformation_gradient: &DeformationGradient| {
416                    hessian(time[0], &state, deformation_gradient)
417                },
418                DeformationGradient::identity(),
419                equality_constraint(time[0]),
420                None,
421            )
422            .map_err(|error| ConstitutiveError::upstream(String::from(error), self))?;
423        let mut times = Times::new();
424        let mut deformation_gradients = DeformationGradients::new();
425        let mut state_variables = ViscoplasticStateVariablesHistory::new();
426        let mut carry = None;
427        times.push(time[0]);
428        deformation_gradients.push(deformation_gradient.clone());
429        state_variables.push(state.clone());
430        for step in time.windows(2) {
431            let advanced = rkmk_dae_step_second_order_minimize::<
432                <Self as StateEvolution<Time, Y>>::Field,
433                Tab,
434                Quantity<EnergyDensity>,
435                FirstPiolaKirchhoffStress,
436                FirstPiolaKirchhoffTangentStiffness,
437                DeformationGradient,
438                Time,
439            >(
440                &mut |t, state, deformation_gradient| {
441                    self.state_rate(t, deformation_gradient, state)
442                },
443                function,
444                jacobian,
445                hessian,
446                &solver,
447                &state,
448                &deformation_gradient,
449                step[0],
450                step[1] - step[0],
451                &mut scratch,
452                carry.as_ref(),
453                equality_constraint,
454                None,
455            )
456            .map_err(|error| ConstitutiveError::upstream(error, self))?;
457            state = advanced.0;
458            deformation_gradient = advanced.1;
459            carry = advanced.2;
460            times.push(step[1]);
461            deformation_gradients.push(deformation_gradient.clone());
462            state_variables.push(state.clone());
463        }
464        Ok((times, deformation_gradients, state_variables))
465    }
466    #[allow(clippy::type_complexity)]
467    fn root_rkmk_dae_adaptive_minimize<Tab: EmbeddedTableau>(
468        &self,
469        applied_load: AppliedLoad,
470        solver: impl SecondOrderOptimization<
471            Quantity<EnergyDensity>,
472            FirstPiolaKirchhoffStress,
473            FirstPiolaKirchhoffTangentStiffness,
474            DeformationGradient,
475        >,
476        abs_tol: Scalar,
477        rel_tol: Scalar,
478    ) -> Result<
479        (
480            Times,
481            DeformationGradients,
482            ViscoplasticStateVariablesHistory<Y>,
483        ),
484        ConstitutiveError,
485    > {
486        let (matrix, prescribed, time) = bcs(applied_load);
487        let state = <Self as StateEvolution<Time, Y>>::initial_state(self);
488        let equality_constraint = |t: Quantity<Time>| {
489            let mut vector = Vector::zero(matrix.len());
490            prescribed
491                .iter()
492                .for_each(|(index, function)| vector[*index] = function(t));
493            EqualityConstraint::Linear(matrix.clone(), vector)
494        };
495        let function = |_: Quantity<Time>,
496                        state: &ViscoplasticStateVariables<Y>,
497                        deformation_gradient: &DeformationGradient|
498         -> Result<Quantity<EnergyDensity>, String> {
499            Ok(self.helmholtz_free_energy_density(deformation_gradient, &state.0)?)
500        };
501        let jacobian = |_: Quantity<Time>,
502                        state: &ViscoplasticStateVariables<Y>,
503                        deformation_gradient: &DeformationGradient|
504         -> Result<FirstPiolaKirchhoffStress, String> {
505            Ok(self.first_piola_kirchhoff_stress(deformation_gradient, &state.0)?)
506        };
507        let hessian = |_: Quantity<Time>,
508                       state: &ViscoplasticStateVariables<Y>,
509                       deformation_gradient: &DeformationGradient|
510         -> Result<FirstPiolaKirchhoffTangentStiffness, String> {
511            Ok(self.first_piola_kirchhoff_tangent_stiffness(deformation_gradient, &state.0)?)
512        };
513        let deformation_gradient = solver
514            .minimize(
515                |deformation_gradient: &DeformationGradient| {
516                    function(time[0], &state, deformation_gradient)
517                },
518                |deformation_gradient: &DeformationGradient| {
519                    jacobian(time[0], &state, deformation_gradient)
520                },
521                |deformation_gradient: &DeformationGradient| {
522                    hessian(time[0], &state, deformation_gradient)
523                },
524                DeformationGradient::identity(),
525                equality_constraint(time[0]),
526                None,
527            )
528            .map_err(|error| ConstitutiveError::upstream(String::from(error), self))?;
529        let (times, state_variables, deformation_gradients) =
530            integrate_rkmk_dae_adaptive_second_order_minimize::<
531                <Self as StateEvolution<Time, Y>>::Field,
532                Tab,
533                Quantity<EnergyDensity>,
534                FirstPiolaKirchhoffStress,
535                FirstPiolaKirchhoffTangentStiffness,
536                DeformationGradient,
537                ViscoplasticStateVariablesHistory<Y>,
538                DeformationGradients,
539                Time,
540            >(
541                |t, state, deformation_gradient| self.state_rate(t, deformation_gradient, state),
542                function,
543                jacobian,
544                hessian,
545                &solver,
546                time,
547                (state, deformation_gradient),
548                abs_tol,
549                rel_tol,
550                equality_constraint,
551                None,
552            )
553            .map_err(|error| ConstitutiveError::upstream(error, self))?;
554        Ok((times, deformation_gradients, state_variables))
555    }
556}
conspire/constitutive/solid/hyperelastic_viscoplastic/mod.rs

conspire/constitutive/solid/hyperelastic_viscoplastic/
mod.rs
