conspire/constitutive/solid/elastic/internal_variables/
mod.rs

1//! Elastic solid constitutive models with internal variables.
2
3use crate::{
4    constitutive::{
5        ConstitutiveError,
6        solid::{Solid, elastic::AppliedLoad},
7    },
8    math::{
9        ContractFirstSecondIndicesWithSecondIndicesOf, ContractSecondIndexWithFirstIndexOf,
10        IDENTITY, Matrix, Rank2, Tensor, TensorArray, TensorTuple, Vector,
11        optimize::{EqualityConstraint, FirstOrderRootFinding, ZerothOrderRootFinding},
12    },
13    mechanics::{
14        CauchyStress, CauchyTangentStiffness, DeformationGradient, FirstPiolaKirchhoffStress,
15        FirstPiolaKirchhoffTangentStiffness, SecondPiolaKirchhoffStress,
16        SecondPiolaKirchhoffTangentStiffness,
17    },
18};
19
20/// Required methods for elastic solid constitutive models with internal variables.
21pub trait ElasticIV<V, T1, T2, T3>
22where
23    Self: Solid,
24{
25    /// Calculates and returns the Cauchy stress.
26    ///
27    /// ```math
28    /// \boldsymbol{\sigma} = J^{-1}\mathbf{P}\cdot\mathbf{F}^T
29    /// ```
30    fn cauchy_stress(
31        &self,
32        deformation_gradient: &DeformationGradient,
33        internal_variables: &V,
34    ) -> Result<CauchyStress, ConstitutiveError> {
35        Ok(deformation_gradient
36            * self.second_piola_kirchhoff_stress(deformation_gradient, internal_variables)?
37            * deformation_gradient.transpose()
38            / deformation_gradient.determinant())
39    }
40    /// Calculates and returns the tangent stiffness associated with the Cauchy stress.
41    ///
42    /// ```math
43    /// \mathcal{T}_{ijkL} = \frac{\partial\sigma_{ij}}{\partial F_{kL}} = J^{-1} \mathcal{G}_{MNkL} F_{iM} F_{jN} - \sigma_{ij} F_{kL}^{-T} + \left(\delta_{jk}\sigma_{is} + \delta_{ik}\sigma_{js}\right)F_{sL}^{-T}
44    /// ```
45    fn cauchy_tangent_stiffness(
46        &self,
47        deformation_gradient: &DeformationGradient,
48        internal_variables: &V,
49    ) -> Result<CauchyTangentStiffness, ConstitutiveError> {
50        let deformation_gradient_inverse_transpose = deformation_gradient.inverse_transpose();
51        let cauchy_stress = self.cauchy_stress(deformation_gradient, internal_variables)?;
52        let some_stress = &cauchy_stress * &deformation_gradient_inverse_transpose;
53        Ok(self
54            .second_piola_kirchhoff_tangent_stiffness(deformation_gradient, internal_variables)?
55            .contract_first_second_indices_with_second_indices_of(
56                deformation_gradient,
57                deformation_gradient,
58            )
59            / deformation_gradient.determinant()
60            - CauchyTangentStiffness::dyad_ij_kl(
61                &cauchy_stress,
62                &deformation_gradient_inverse_transpose,
63            )
64            + CauchyTangentStiffness::dyad_il_kj(&some_stress, &IDENTITY)
65            + CauchyTangentStiffness::dyad_ik_jl(&IDENTITY, &some_stress))
66    }
67    /// Calculates and returns the first Piola-Kirchhoff stress.
68    ///
69    /// ```math
70    /// \mathbf{P} = J\boldsymbol{\sigma}\cdot\mathbf{F}^{-T}
71    /// ```
72    fn first_piola_kirchhoff_stress(
73        &self,
74        deformation_gradient: &DeformationGradient,
75        internal_variables: &V,
76    ) -> Result<FirstPiolaKirchhoffStress, ConstitutiveError> {
77        Ok(
78            self.cauchy_stress(deformation_gradient, internal_variables)?
79                * deformation_gradient.inverse_transpose()
80                * deformation_gradient.determinant(),
81        )
82    }
83    /// Calculates and returns the tangent stiffness associated with the first Piola-Kirchhoff stress.
84    ///
85    /// ```math
86    /// \mathcal{C}_{iJkL} = \frac{\partial P_{iJ}}{\partial F_{kL}} = J \mathcal{T}_{iskL} F_{sJ}^{-T} + P_{iJ} F_{kL}^{-T} - P_{iL} F_{kJ}^{-T}
87    /// ```
88    fn first_piola_kirchhoff_tangent_stiffness(
89        &self,
90        deformation_gradient: &DeformationGradient,
91        internal_variables: &V,
92    ) -> Result<FirstPiolaKirchhoffTangentStiffness, ConstitutiveError> {
93        let deformation_gradient_inverse_transpose = deformation_gradient.inverse_transpose();
94        let first_piola_kirchhoff_stress =
95            self.first_piola_kirchhoff_stress(deformation_gradient, internal_variables)?;
96        Ok(self
97            .cauchy_tangent_stiffness(deformation_gradient, internal_variables)?
98            .contract_second_index_with_first_index_of(&deformation_gradient_inverse_transpose)
99            * deformation_gradient.determinant()
100            + FirstPiolaKirchhoffTangentStiffness::dyad_ij_kl(
101                &first_piola_kirchhoff_stress,
102                &deformation_gradient_inverse_transpose,
103            )
104            - FirstPiolaKirchhoffTangentStiffness::dyad_il_kj(
105                &first_piola_kirchhoff_stress,
106                &deformation_gradient_inverse_transpose,
107            ))
108    }
109    /// Calculates and returns the second Piola-Kirchhoff stress.
110    ///
111    /// ```math
112    /// \mathbf{S} = \mathbf{F}^{-1}\cdot\mathbf{P}
113    /// ```
114    fn second_piola_kirchhoff_stress(
115        &self,
116        deformation_gradient: &DeformationGradient,
117        internal_variables: &V,
118    ) -> Result<SecondPiolaKirchhoffStress, ConstitutiveError> {
119        Ok(deformation_gradient.inverse()
120            * self.first_piola_kirchhoff_stress(deformation_gradient, internal_variables)?)
121    }
122    /// Calculates and returns the tangent stiffness associated with the second Piola-Kirchhoff stress.
123    ///
124    /// ```math
125    /// \mathcal{G}_{IJkL} = \frac{\partial S_{IJ}}{\partial F_{kL}} = \mathcal{C}_{mJkL}F_{mI}^{-T} - S_{LJ}F_{kI}^{-T} = J \mathcal{T}_{mnkL} F_{mI}^{-T} F_{nJ}^{-T} + S_{IJ} F_{kL}^{-T} - S_{IL} F_{kJ}^{-T} -S_{LJ} F_{kI}^{-T}
126    /// ```
127    fn second_piola_kirchhoff_tangent_stiffness(
128        &self,
129        deformation_gradient: &DeformationGradient,
130        internal_variables: &V,
131    ) -> Result<SecondPiolaKirchhoffTangentStiffness, ConstitutiveError> {
132        let deformation_gradient_inverse_transpose = deformation_gradient.inverse_transpose();
133        let deformation_gradient_inverse = deformation_gradient_inverse_transpose.transpose();
134        let second_piola_kirchhoff_stress =
135            self.second_piola_kirchhoff_stress(deformation_gradient, internal_variables)?;
136        Ok(self
137            .cauchy_tangent_stiffness(deformation_gradient, internal_variables)?
138            .contract_first_second_indices_with_second_indices_of(
139                &deformation_gradient_inverse,
140                &deformation_gradient_inverse,
141            )
142            * deformation_gradient.determinant()
143            + SecondPiolaKirchhoffTangentStiffness::dyad_ij_kl(
144                &second_piola_kirchhoff_stress,
145                &deformation_gradient_inverse_transpose,
146            )
147            - SecondPiolaKirchhoffTangentStiffness::dyad_il_kj(
148                &second_piola_kirchhoff_stress,
149                &deformation_gradient_inverse_transpose,
150            )
151            - SecondPiolaKirchhoffTangentStiffness::dyad_ik_jl(
152                &deformation_gradient_inverse,
153                &second_piola_kirchhoff_stress,
154            ))
155    }
156    /// Returns the initial value for the internal variables.
157    fn internal_variables_initial(&self) -> V;
158    /// Calculates and returns the residual associated with the internal variables.
159    fn internal_variables_residual(
160        &self,
161        deformation_gradient: &DeformationGradient,
162        internal_variables: &V,
163    ) -> Result<V, ConstitutiveError>;
164    /// Calculates and returns the tangents associated with the internal variables.
165    fn internal_variables_tangents(
166        &self,
167        deformation_gradient: &DeformationGradient,
168        internal_variables: &V,
169    ) -> Result<(T1, T2, T3), ConstitutiveError>;
170    /// Returns the constraint indices for the internal variables.
171    fn internal_variables_constraints(&self) -> (&[usize], usize);
172}
173
174/// Zeroth-order root-finding methods for elastic solid constitutive models with internal variables.
175pub trait ZerothOrderRoot<V, T1, T2, T3>
176where
177    V: Tensor,
178{
179    /// Type representing all variables.
180    type Variables;
181    /// Solve for the unknown components of the deformation gradient under an applied load.
182    ///
183    /// ```math
184    /// \mathbf{P}(\mathbf{F}) - \boldsymbol{\lambda} - \mathbf{P}_0 = \mathbf{0}
185    /// ```
186    fn root(
187        &self,
188        applied_load: AppliedLoad,
189        solver: impl ZerothOrderRootFinding<Self::Variables>,
190    ) -> Result<(DeformationGradient, V), ConstitutiveError>;
191}
192
193/// First-order root-finding methods for elastic solid constitutive models with internal variables.
194pub trait FirstOrderRoot<V, T1, T2, T3>
195where
196    T1: Tensor,
197    T2: Tensor,
198    T3: Tensor,
199    V: Tensor,
200{
201    /// Type representing all variables.
202    type Variables;
203    /// Solve for the unknown components of the deformation gradient under an applied load.
204    ///
205    /// ```math
206    /// \mathbf{P}(\mathbf{F}) - \boldsymbol{\lambda} - \mathbf{P}_0 = \mathbf{0}
207    /// ```
208    fn root(
209        &self,
210        applied_load: AppliedLoad,
211        solver: impl FirstOrderRootFinding<
212            Self::Variables,
213            TensorTuple<FirstPiolaKirchhoffTangentStiffness, TensorTuple<T1, TensorTuple<T2, T3>>>,
214            Self::Variables,
215        >,
216    ) -> Result<(DeformationGradient, V), ConstitutiveError>;
217}
218
219impl<T, V, T1, T2, T3> ZerothOrderRoot<V, T1, T2, T3> for T
220where
221    T: ElasticIV<V, T1, T2, T3>,
222    V: Tensor,
223{
224    type Variables = TensorTuple<DeformationGradient, V>;
225    fn root(
226        &self,
227        applied_load: AppliedLoad,
228        solver: impl ZerothOrderRootFinding<Self::Variables>,
229    ) -> Result<(DeformationGradient, V), ConstitutiveError> {
230        let (matrix, vector) = bcs(self, applied_load);
231        match solver.root(
232            |variables: &Self::Variables| {
233                let (deformation_gradient, internal_variables) = variables.into();
234                Ok(TensorTuple::from((
235                    self.first_piola_kirchhoff_stress(deformation_gradient, internal_variables)?,
236                    self.internal_variables_residual(deformation_gradient, internal_variables)?,
237                )))
238            },
239            Self::Variables::from((
240                DeformationGradient::identity(),
241                self.internal_variables_initial(),
242            )),
243            EqualityConstraint::Linear(matrix, vector),
244        ) {
245            Ok(solution) => Ok(solution.into()),
246            Err(error) => Err(ConstitutiveError::Upstream(
247                format!("{error}"),
248                format!("{self:?}"),
249            )),
250        }
251    }
252}
253
254impl<T, V, T1, T2, T3> FirstOrderRoot<V, T1, T2, T3> for T
255where
256    T1: Tensor,
257    T2: Tensor,
258    T3: Tensor,
259    T: ElasticIV<V, T1, T2, T3>,
260    V: Tensor,
261{
262    type Variables = TensorTuple<DeformationGradient, V>;
263    fn root(
264        &self,
265        applied_load: AppliedLoad,
266        solver: impl FirstOrderRootFinding<
267            Self::Variables,
268            TensorTuple<FirstPiolaKirchhoffTangentStiffness, TensorTuple<T1, TensorTuple<T2, T3>>>,
269            Self::Variables,
270        >,
271    ) -> Result<(DeformationGradient, V), ConstitutiveError> {
272        let (matrix, vector) = bcs(self, applied_load);
273        match solver.root(
274            |variables: &Self::Variables| {
275                let (deformation_gradient, internal_variables) = variables.into();
276                Ok(TensorTuple::from((
277                    self.first_piola_kirchhoff_stress(deformation_gradient, internal_variables)?,
278                    self.internal_variables_residual(deformation_gradient, internal_variables)?,
279                )))
280            },
281            |variables: &Self::Variables| {
282                let (deformation_gradient, internal_variables) = variables.into();
283                let tangent_0 = self.first_piola_kirchhoff_tangent_stiffness(
284                    deformation_gradient,
285                    internal_variables,
286                )?;
287                let (tangent_1, tangent_2, tangent_3) =
288                    self.internal_variables_tangents(deformation_gradient, internal_variables)?;
289                Ok((tangent_0, (tangent_1, (tangent_2, tangent_3).into()).into()).into())
290            },
291            Self::Variables::from((
292                DeformationGradient::identity(),
293                self.internal_variables_initial(),
294            )),
295            EqualityConstraint::Linear(matrix, vector),
296        ) {
297            Ok(solution) => Ok(solution.into()),
298            Err(error) => Err(ConstitutiveError::Upstream(
299                format!("{error}"),
300                format!("{self:?}"),
301            )),
302        }
303    }
304}
305
306#[doc(hidden)]
307pub fn bcs<C, V, T1, T2, T3>(model: &C, applied_load: AppliedLoad) -> (Matrix, Vector)
308where
309    C: ElasticIV<V, T1, T2, T3>,
310{
311    let (extra_constraints, num_vars) = model.internal_variables_constraints();
312    match applied_load {
313        AppliedLoad::UniaxialStress(deformation_gradient_11) => {
314            let num_constraints = 4;
315            let num_constraints_vars = extra_constraints.len();
316            let mut matrix = Matrix::zero(num_constraints + num_constraints_vars, 9 + num_vars);
317            let mut vector = Vector::zero(num_constraints + num_constraints_vars);
318            matrix[0][0] = 1.0;
319            matrix[1][1] = 1.0;
320            matrix[2][2] = 1.0;
321            matrix[3][5] = 1.0;
322            let mut i = num_constraints;
323            extra_constraints.iter().for_each(|&j| {
324                matrix[i][j] = 1.0;
325                i += 1;
326            });
327            vector[0] = deformation_gradient_11;
328            (matrix, vector)
329        }
330        AppliedLoad::BiaxialStress(_, _) => {
331            todo!()
332        }
333    }
334}
conspire/constitutive/solid/elastic/internal_variables/mod.rs

conspire/constitutive/solid/elastic/internal_variables/
mod.rs
