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

1//! Elastic 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 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_foo(
31        &self,
32        deformation_gradient: &DeformationGradient,
33        internal_variables: &V,
34    ) -> Result<CauchyStress, ConstitutiveError> {
35        Ok(deformation_gradient
36            * self.second_piola_kirchhoff_stress_foo(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_foo(
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_foo(deformation_gradient, internal_variables)?;
52        let some_stress = &cauchy_stress * &deformation_gradient_inverse_transpose;
53        Ok(self
54            .second_piola_kirchhoff_tangent_stiffness_foo(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_foo(
73        &self,
74        deformation_gradient: &DeformationGradient,
75        internal_variables: &V,
76    ) -> Result<FirstPiolaKirchhoffStress, ConstitutiveError> {
77        Ok(
78            self.cauchy_stress_foo(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_foo(
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_foo(deformation_gradient, internal_variables)?;
96        Ok(self
97            .cauchy_tangent_stiffness_foo(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_foo(
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_foo(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_foo(
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_foo(deformation_gradient, internal_variables)?;
136        Ok(self
137            .cauchy_tangent_stiffness_foo(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 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 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 (extra_constraints, num_vars) = self.internal_variables_constraints();
231        let (matrix, vector) = match applied_load {
232            AppliedLoad::UniaxialStress(deformation_gradient_11) => {
233                let num_constraints = 4;
234                let num_constraints_vars = extra_constraints.len();
235                let mut matrix = Matrix::zero(num_constraints + num_constraints_vars, 9 + num_vars);
236                let mut vector = Vector::zero(num_constraints + num_constraints_vars);
237                matrix[0][0] = 1.0;
238                matrix[1][1] = 1.0;
239                matrix[2][2] = 1.0;
240                matrix[3][5] = 1.0;
241                let mut i = num_constraints;
242                extra_constraints.iter().for_each(|&j| {
243                    matrix[i][j] = 1.0;
244                    i += 1;
245                });
246                vector[0] = deformation_gradient_11;
247                (matrix, vector)
248            }
249            AppliedLoad::BiaxialStress(_, _) => {
250                todo!()
251            }
252        };
253        match solver.root(
254            |variables: &Self::Variables| {
255                let (deformation_gradient, internal_variables) = variables.into();
256                Ok(TensorTuple::from((
257                    self.first_piola_kirchhoff_stress_foo(
258                        deformation_gradient,
259                        internal_variables,
260                    )?,
261                    self.internal_variables_residual(deformation_gradient, internal_variables)?,
262                )))
263            },
264            Self::Variables::from((
265                DeformationGradient::identity(),
266                self.internal_variables_initial(),
267            )),
268            EqualityConstraint::Linear(matrix, vector),
269        ) {
270            Ok(solution) => Ok(solution.into()),
271            Err(error) => Err(ConstitutiveError::Upstream(
272                format!("{error}"),
273                format!("{self:?}"),
274            )),
275        }
276    }
277}
278
279impl<T, V, T1, T2, T3> FirstOrderRoot<V, T1, T2, T3> for T
280where
281    T1: Tensor,
282    T2: Tensor,
283    T3: Tensor,
284    T: ElasticIV<V, T1, T2, T3>,
285    V: Tensor,
286{
287    type Variables = TensorTuple<DeformationGradient, V>;
288    fn root(
289        &self,
290        applied_load: AppliedLoad,
291        solver: impl FirstOrderRootFinding<
292            Self::Variables,
293            TensorTuple<FirstPiolaKirchhoffTangentStiffness, TensorTuple<T1, TensorTuple<T2, T3>>>,
294            Self::Variables,
295        >,
296    ) -> Result<(DeformationGradient, V), ConstitutiveError> {
297        let (extra_constraints, num_vars) = self.internal_variables_constraints();
298        let (matrix, vector) = match applied_load {
299            AppliedLoad::UniaxialStress(deformation_gradient_11) => {
300                let num_constraints = 4;
301                let num_constraints_vars = extra_constraints.len();
302                let mut matrix = Matrix::zero(num_constraints + num_constraints_vars, 9 + num_vars);
303                let mut vector = Vector::zero(num_constraints + num_constraints_vars);
304                matrix[0][0] = 1.0;
305                matrix[1][1] = 1.0;
306                matrix[2][2] = 1.0;
307                matrix[3][5] = 1.0;
308                let mut i = num_constraints;
309                extra_constraints.iter().for_each(|&j| {
310                    matrix[i][j] = 1.0;
311                    i += 1;
312                });
313                vector[0] = deformation_gradient_11;
314                (matrix, vector)
315            }
316            AppliedLoad::BiaxialStress(_, _) => {
317                todo!()
318            }
319        };
320        match solver.root(
321            |variables: &Self::Variables| {
322                let (deformation_gradient, internal_variables) = variables.into();
323                Ok(TensorTuple::from((
324                    self.first_piola_kirchhoff_stress_foo(
325                        deformation_gradient,
326                        internal_variables,
327                    )?,
328                    self.internal_variables_residual(deformation_gradient, internal_variables)?,
329                )))
330            },
331            |variables: &Self::Variables| {
332                let (deformation_gradient, internal_variables) = variables.into();
333                let tangent_0 = self.first_piola_kirchhoff_tangent_stiffness_foo(
334                    deformation_gradient,
335                    internal_variables,
336                )?;
337                let (tangent_1, tangent_2, tangent_3) =
338                    self.internal_variables_tangents(deformation_gradient, internal_variables)?;
339                Ok((tangent_0, (tangent_1, (tangent_2, tangent_3).into()).into()).into())
340            },
341            Self::Variables::from((
342                DeformationGradient::identity(),
343                self.internal_variables_initial(),
344            )),
345            EqualityConstraint::Linear(matrix, vector),
346        ) {
347            Ok(solution) => Ok(solution.into()),
348            Err(error) => Err(ConstitutiveError::Upstream(
349                format!("{error}"),
350                format!("{self:?}"),
351            )),
352        }
353    }
354}
conspire/constitutive/solid/elastic/internal_variables/mod.rs

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