conspire/constitutive/solid/hyperelastic/
mod.rs

1//! Hyperelastic constitutive models.
2//!
3//! ---
4//!
5#![doc = include_str!("doc.md")]
6
7#[cfg(feature = "doc")]
8pub mod doc;
9
10#[cfg(test)]
11pub mod test;
12
13mod arruda_boyce;
14mod fung;
15mod gent;
16mod hencky;
17mod mooney_rivlin;
18mod neo_hookean;
19mod saint_venant_kirchhoff;
20mod yeoh;
21
22pub use self::{
23    arruda_boyce::ArrudaBoyce, fung::Fung, gent::Gent, hencky::Hencky, mooney_rivlin::MooneyRivlin,
24    neo_hookean::NeoHookean, saint_venant_kirchhoff::SaintVenantKirchhoff, yeoh::Yeoh,
25};
26use super::{
27    elastic::{AppliedLoad, Elastic},
28    *,
29};
30use crate::math::{
31    Matrix, TensorVec, Vector,
32    optimize::{EqualityConstraint, FirstOrderOptimization, SecondOrderOptimization},
33};
34
35/// Required methods for hyperelastic constitutive models.
36pub trait Hyperelastic
37where
38    Self: Elastic,
39{
40    /// Calculates and returns the Helmholtz free energy density.
41    ///
42    /// ```math
43    /// a = a(\mathbf{F})
44    /// ```
45    fn helmholtz_free_energy_density(
46        &self,
47        deformation_gradient: &DeformationGradient,
48    ) -> Result<Scalar, ConstitutiveError>;
49}
50
51/// First-order minimization methods for elastic constitutive models.
52pub trait FirstOrderMinimize {
53    /// Solve for the unknown components of the deformation gradient under an applied load.
54    ///
55    /// ```math
56    /// \Pi(\mathbf{F},\boldsymbol{\lambda}) = a(\mathbf{F}) - \boldsymbol{\lambda}:(\mathbf{F} - \mathbf{F}_0) - \mathbf{P}_0:\mathbf{F}
57    /// ```
58    fn minimize(
59        &self,
60        applied_load: AppliedLoad,
61        solver: impl FirstOrderOptimization<Scalar, DeformationGradient>,
62    ) -> Result<DeformationGradient, ConstitutiveError>;
63}
64
65/// Second-order minimization methods for elastic constitutive models.
66pub trait SecondOrderMinimize {
67    /// Solve for the unknown components of the deformation gradient under an applied load.
68    ///
69    /// ```math
70    /// \Pi(\mathbf{F},\boldsymbol{\lambda}) = a(\mathbf{F}) - \boldsymbol{\lambda}:(\mathbf{F} - \mathbf{F}_0) - \mathbf{P}_0:\mathbf{F}
71    /// ```
72    fn minimize(
73        &self,
74        applied_load: AppliedLoad,
75        solver: impl SecondOrderOptimization<
76            Scalar,
77            FirstPiolaKirchhoffStress,
78            FirstPiolaKirchhoffTangentStiffness,
79            DeformationGradient,
80        >,
81    ) -> Result<DeformationGradient, ConstitutiveError>;
82}
83
84impl<T> FirstOrderMinimize for T
85where
86    T: Hyperelastic,
87{
88    fn minimize(
89        &self,
90        applied_load: AppliedLoad,
91        solver: impl FirstOrderOptimization<Scalar, DeformationGradient>,
92    ) -> Result<DeformationGradient, ConstitutiveError> {
93        match match applied_load {
94            AppliedLoad::UniaxialStress(deformation_gradient_11) => {
95                let mut matrix = Matrix::zero(4, 9);
96                let mut vector = Vector::zero(4);
97                matrix[0][0] = 1.0;
98                matrix[1][1] = 1.0;
99                matrix[2][2] = 1.0;
100                matrix[3][5] = 1.0;
101                vector[0] = deformation_gradient_11;
102                solver.minimize(
103                    |deformation_gradient: &DeformationGradient| {
104                        Ok(self.helmholtz_free_energy_density(deformation_gradient)?)
105                    },
106                    |deformation_gradient: &DeformationGradient| {
107                        Ok(self.first_piola_kirchhoff_stress(deformation_gradient)?)
108                    },
109                    DeformationGradient::identity(),
110                    EqualityConstraint::Linear(matrix, vector),
111                )
112            }
113            AppliedLoad::BiaxialStress(deformation_gradient_11, deformation_gradient_22) => {
114                let mut matrix = Matrix::zero(5, 9);
115                let mut vector = Vector::zero(5);
116                matrix[0][0] = 1.0;
117                matrix[1][1] = 1.0;
118                matrix[2][2] = 1.0;
119                matrix[3][5] = 1.0;
120                matrix[4][4] = 1.0;
121                vector[0] = deformation_gradient_11;
122                vector[4] = deformation_gradient_22;
123                solver.minimize(
124                    |deformation_gradient: &DeformationGradient| {
125                        Ok(self.helmholtz_free_energy_density(deformation_gradient)?)
126                    },
127                    |deformation_gradient: &DeformationGradient| {
128                        Ok(self.first_piola_kirchhoff_stress(deformation_gradient)?)
129                    },
130                    DeformationGradient::identity(),
131                    EqualityConstraint::Linear(matrix, vector),
132                )
133            }
134        } {
135            Ok(deformation_gradient) => Ok(deformation_gradient),
136            Err(error) => Err(ConstitutiveError::Upstream(
137                format!("{error}"),
138                format!("{self:?}"),
139            )),
140        }
141    }
142}
143
144impl<T> SecondOrderMinimize for T
145where
146    T: Hyperelastic,
147{
148    fn minimize(
149        &self,
150        applied_load: AppliedLoad,
151        solver: impl SecondOrderOptimization<
152            Scalar,
153            FirstPiolaKirchhoffStress,
154            FirstPiolaKirchhoffTangentStiffness,
155            DeformationGradient,
156        >,
157    ) -> Result<DeformationGradient, ConstitutiveError> {
158        match match applied_load {
159            AppliedLoad::UniaxialStress(deformation_gradient_11) => {
160                let mut matrix = Matrix::zero(4, 9);
161                let mut vector = Vector::zero(4);
162                matrix[0][0] = 1.0;
163                matrix[1][1] = 1.0;
164                matrix[2][2] = 1.0;
165                matrix[3][5] = 1.0;
166                vector[0] = deformation_gradient_11;
167                solver.minimize(
168                    |deformation_gradient: &DeformationGradient| {
169                        Ok(self.helmholtz_free_energy_density(deformation_gradient)?)
170                    },
171                    |deformation_gradient: &DeformationGradient| {
172                        Ok(self.first_piola_kirchhoff_stress(deformation_gradient)?)
173                    },
174                    |deformation_gradient: &DeformationGradient| {
175                        Ok(self.first_piola_kirchhoff_tangent_stiffness(deformation_gradient)?)
176                    },
177                    DeformationGradient::identity(),
178                    EqualityConstraint::Linear(matrix, vector),
179                    None,
180                )
181            }
182            AppliedLoad::BiaxialStress(deformation_gradient_11, deformation_gradient_22) => {
183                let mut matrix = Matrix::zero(5, 9);
184                let mut vector = Vector::zero(5);
185                matrix[0][0] = 1.0;
186                matrix[1][1] = 1.0;
187                matrix[2][2] = 1.0;
188                matrix[3][5] = 1.0;
189                matrix[4][4] = 1.0;
190                vector[0] = deformation_gradient_11;
191                vector[4] = deformation_gradient_22;
192                solver.minimize(
193                    |deformation_gradient: &DeformationGradient| {
194                        Ok(self.helmholtz_free_energy_density(deformation_gradient)?)
195                    },
196                    |deformation_gradient: &DeformationGradient| {
197                        Ok(self.first_piola_kirchhoff_stress(deformation_gradient)?)
198                    },
199                    |deformation_gradient: &DeformationGradient| {
200                        Ok(self.first_piola_kirchhoff_tangent_stiffness(deformation_gradient)?)
201                    },
202                    DeformationGradient::identity(),
203                    EqualityConstraint::Linear(matrix, vector),
204                    None,
205                )
206            }
207        } {
208            Ok(deformation_gradient) => Ok(deformation_gradient),
209            Err(error) => Err(ConstitutiveError::Upstream(
210                format!("{error}"),
211                format!("{self:?}"),
212            )),
213        }
214    }
215}