conspire/constitutive/solid/hyperelastic/
mod.rs

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