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        ContractFirstSecondWithSecond, ContractSecondWithFirst, Hessian, HessianBlock, IDENTITY,
10        Jacobian, Matrix, Rank2, Tensor, TensorArray, TensorTuple, Vector,
11        optimize::{
12            EqualityConstraint, FirstOrderRootFindingBlock, SolveStrategy, ZerothOrderRootFinding,
13        },
14        sparse::CscMatrix,
15    },
16    mechanics::{
17        CauchyStress, CauchyTangentStiffness, DeformationGradient, FirstPiolaKirchhoffStress,
18        FirstPiolaKirchhoffTangentStiffness, SecondPiolaKirchhoffStress,
19        SecondPiolaKirchhoffTangentStiffness,
20    },
21};
22
23/// The tangents of the coupled system, in the order the block solver takes them.
24pub type Tangents<C, V> = (
25    FirstPiolaKirchhoffTangentStiffness,
26    <C as ElasticIV<V>>::TangentVu,
27    <C as ElasticIV<V>>::TangentUv,
28    <C as ElasticIV<V>>::TangentVv,
29);
30
31/// Required methods for elastic solid constitutive models with internal variables.
32pub trait ElasticIV<V>
33where
34    Self: Solid,
35{
36    /// The residual associated with the internal variables.
37    /// the internal variables are in equilibrium over, typically a stress.
38    type Residual: Jacobian;
39    /// The tangent of the internal variables residual with the deformation gradient.
40    type TangentVu: HessianBlock;
41    /// The tangent of the deformation gradient residual with the internal variables.
42    type TangentUv: HessianBlock;
43    /// The tangent of the internal variables residual with the internal variables.
44    type TangentVv: Hessian + HessianBlock;
45    /// Calculates and returns the Cauchy stress.
46    ///
47    /// ```math
48    /// \boldsymbol{\sigma} = J^{-1}\mathbf{P}\cdot\mathbf{F}^T
49    /// ```
50    fn cauchy_stress(
51        &self,
52        deformation_gradient: &DeformationGradient,
53        internal_variables: &V,
54    ) -> Result<CauchyStress, ConstitutiveError> {
55        Ok(deformation_gradient
56            * self.second_piola_kirchhoff_stress(deformation_gradient, internal_variables)?
57            * deformation_gradient.transpose()
58            / deformation_gradient.determinant())
59    }
60    /// Calculates and returns the tangent stiffness associated with the Cauchy stress.
61    ///
62    /// ```math
63    /// \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}
64    /// ```
65    fn cauchy_tangent_stiffness(
66        &self,
67        deformation_gradient: &DeformationGradient,
68        internal_variables: &V,
69    ) -> Result<CauchyTangentStiffness, ConstitutiveError> {
70        let deformation_gradient_inverse_transpose = deformation_gradient.inverse_transpose();
71        let cauchy_stress = self.cauchy_stress(deformation_gradient, internal_variables)?;
72        let some_stress = &cauchy_stress * &deformation_gradient_inverse_transpose;
73        Ok(self
74            .second_piola_kirchhoff_tangent_stiffness(deformation_gradient, internal_variables)?
75            .contract_first_second_with_second(deformation_gradient, deformation_gradient)
76            / deformation_gradient.determinant()
77            - CauchyTangentStiffness::dyad_ij_kl(
78                &cauchy_stress,
79                &deformation_gradient_inverse_transpose,
80            )
81            + CauchyTangentStiffness::dyad_il_kj(&some_stress, &IDENTITY)
82            + CauchyTangentStiffness::dyad_ik_jl(&IDENTITY, &some_stress))
83    }
84    /// Calculates and returns the first Piola-Kirchhoff stress.
85    ///
86    /// ```math
87    /// \mathbf{P} = J\boldsymbol{\sigma}\cdot\mathbf{F}^{-T}
88    /// ```
89    fn first_piola_kirchhoff_stress(
90        &self,
91        deformation_gradient: &DeformationGradient,
92        internal_variables: &V,
93    ) -> Result<FirstPiolaKirchhoffStress, ConstitutiveError> {
94        Ok(
95            self.cauchy_stress(deformation_gradient, internal_variables)?
96                * deformation_gradient.inverse_transpose()
97                * deformation_gradient.determinant(),
98        )
99    }
100    /// Calculates and returns the tangent stiffness associated with the first Piola-Kirchhoff stress.
101    ///
102    /// ```math
103    /// \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}
104    /// ```
105    fn first_piola_kirchhoff_tangent_stiffness(
106        &self,
107        deformation_gradient: &DeformationGradient,
108        internal_variables: &V,
109    ) -> Result<FirstPiolaKirchhoffTangentStiffness, ConstitutiveError> {
110        let deformation_gradient_inverse_transpose = deformation_gradient.inverse_transpose();
111        let first_piola_kirchhoff_stress =
112            self.first_piola_kirchhoff_stress(deformation_gradient, internal_variables)?;
113        Ok(self
114            .cauchy_tangent_stiffness(deformation_gradient, internal_variables)?
115            .contract_second_with_first(&deformation_gradient_inverse_transpose)
116            * deformation_gradient.determinant()
117            + FirstPiolaKirchhoffTangentStiffness::dyad_ij_kl(
118                &first_piola_kirchhoff_stress,
119                &deformation_gradient_inverse_transpose,
120            )
121            - FirstPiolaKirchhoffTangentStiffness::dyad_il_kj(
122                &first_piola_kirchhoff_stress,
123                &deformation_gradient_inverse_transpose,
124            ))
125    }
126    /// Calculates and returns the second Piola-Kirchhoff stress.
127    ///
128    /// ```math
129    /// \mathbf{S} = \mathbf{F}^{-1}\cdot\mathbf{P}
130    /// ```
131    fn second_piola_kirchhoff_stress(
132        &self,
133        deformation_gradient: &DeformationGradient,
134        internal_variables: &V,
135    ) -> Result<SecondPiolaKirchhoffStress, ConstitutiveError> {
136        Ok(deformation_gradient.inverse()
137            * self.first_piola_kirchhoff_stress(deformation_gradient, internal_variables)?)
138    }
139    /// Calculates and returns the tangent stiffness associated with the second Piola-Kirchhoff stress.
140    ///
141    /// ```math
142    /// \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}
143    /// ```
144    fn second_piola_kirchhoff_tangent_stiffness(
145        &self,
146        deformation_gradient: &DeformationGradient,
147        internal_variables: &V,
148    ) -> Result<SecondPiolaKirchhoffTangentStiffness, ConstitutiveError> {
149        let deformation_gradient_inverse_transpose = deformation_gradient.inverse_transpose();
150        let deformation_gradient_inverse = deformation_gradient_inverse_transpose.transpose();
151        let second_piola_kirchhoff_stress =
152            self.second_piola_kirchhoff_stress(deformation_gradient, internal_variables)?;
153        Ok(self
154            .cauchy_tangent_stiffness(deformation_gradient, internal_variables)?
155            .contract_first_second_with_second(
156                &deformation_gradient_inverse,
157                &deformation_gradient_inverse,
158            )
159            * deformation_gradient.determinant()
160            + SecondPiolaKirchhoffTangentStiffness::dyad_ij_kl(
161                &second_piola_kirchhoff_stress,
162                &deformation_gradient_inverse_transpose,
163            )
164            - SecondPiolaKirchhoffTangentStiffness::dyad_il_kj(
165                &second_piola_kirchhoff_stress,
166                &deformation_gradient_inverse_transpose,
167            )
168            - SecondPiolaKirchhoffTangentStiffness::dyad_ik_jl(
169                &deformation_gradient_inverse,
170                &second_piola_kirchhoff_stress,
171            ))
172    }
173    /// Returns the initial value for the internal variables.
174    fn internal_variables_initial(&self) -> V;
175    /// Calculates and returns the residual associated with the internal variables.
176    fn internal_variables_residual(
177        &self,
178        deformation_gradient: &DeformationGradient,
179        internal_variables: &V,
180    ) -> Result<Self::Residual, ConstitutiveError>;
181    /// Returns the indices of the internal variables held at zero.
182    fn internal_variables_fixed(&self) -> &[usize];
183    /// Calculates and returns the tangents of the coupled system.
184    fn tangents(
185        &self,
186        deformation_gradient: &DeformationGradient,
187        internal_variables: &V,
188    ) -> Result<Tangents<Self, V>, ConstitutiveError>;
189}
190
191/// Zeroth-order root-finding methods for elastic solid constitutive models with internal variables.
192pub trait ZerothOrderRoot<V>
193where
194    V: Tensor,
195{
196    /// Type representing all residuals.
197    type Residuals;
198    /// Type representing all variables.
199    type Variables;
200    /// Solve for the unknown components of the deformation gradient under an applied load.
201    ///
202    /// ```math
203    /// \mathbf{P}(\mathbf{F}) - \boldsymbol{\lambda} - \mathbf{P}_0 = \mathbf{0}
204    /// ```
205    fn root(
206        &self,
207        applied_load: AppliedLoad,
208        solver: impl ZerothOrderRootFinding<Self::Residuals, Self::Variables>,
209    ) -> Result<(DeformationGradient, V), ConstitutiveError>;
210}
211
212/// First-order root-finding methods for elastic solid constitutive models with internal variables.
213pub trait FirstOrderRoot<V>
214where
215    Self: ElasticIV<V>,
216    V: Tensor,
217{
218    /// Solve for the unknown components of the deformation gradient under an applied load.
219    ///
220    /// ```math
221    /// \mathbf{P}(\mathbf{F}) - \boldsymbol{\lambda} - \mathbf{P}_0 = \mathbf{0}
222    /// ```
223    fn root(
224        &self,
225        applied_load: AppliedLoad,
226        solver: impl FirstOrderRootFindingBlock<
227            DeformationGradient,
228            V,
229            FirstPiolaKirchhoffStress,
230            <Self as ElasticIV<V>>::Residual,
231            FirstPiolaKirchhoffTangentStiffness,
232            Self::TangentVu,
233            Self::TangentUv,
234            Self::TangentVv,
235        >,
236        strategy: SolveStrategy,
237    ) -> Result<(DeformationGradient, V), ConstitutiveError>;
238}
239
240impl<T, V> ZerothOrderRoot<V> for T
241where
242    T: ElasticIV<V>,
243    V: Tensor,
244{
245    type Residuals = TensorTuple<FirstPiolaKirchhoffStress, <T as ElasticIV<V>>::Residual>;
246    type Variables = TensorTuple<DeformationGradient, V>;
247    fn root(
248        &self,
249        applied_load: AppliedLoad,
250        solver: impl ZerothOrderRootFinding<Self::Residuals, Self::Variables>,
251    ) -> Result<(DeformationGradient, V), ConstitutiveError> {
252        let (matrix, vector) = bcs(self, applied_load);
253        let solution = solver
254            .root(
255                |variables: &Self::Variables| {
256                    let (deformation_gradient, internal_variables) = variables.into();
257                    Ok(TensorTuple::from((
258                        self.first_piola_kirchhoff_stress(
259                            deformation_gradient,
260                            internal_variables,
261                        )?,
262                        self.internal_variables_residual(deformation_gradient, internal_variables)?,
263                    )))
264                },
265                Self::Variables::from((
266                    DeformationGradient::identity(),
267                    self.internal_variables_initial(),
268                )),
269                EqualityConstraint::Linear(matrix, vector),
270            )
271            .map_err(|error| ConstitutiveError::upstream(error, self))?;
272        Ok(solution.into())
273    }
274}
275
276impl<T, V> FirstOrderRoot<V> for T
277where
278    T: ElasticIV<V>,
279    V: Tensor,
280{
281    fn root(
282        &self,
283        applied_load: AppliedLoad,
284        solver: impl FirstOrderRootFindingBlock<
285            DeformationGradient,
286            V,
287            FirstPiolaKirchhoffStress,
288            <Self as ElasticIV<V>>::Residual,
289            FirstPiolaKirchhoffTangentStiffness,
290            Self::TangentVu,
291            Self::TangentUv,
292            Self::TangentVv,
293        >,
294        strategy: SolveStrategy,
295    ) -> Result<(DeformationGradient, V), ConstitutiveError> {
296        let (constraint_global, constraint_local) = bcs_block(self, applied_load);
297        solver
298            .root_block(
299                |deformation_gradient: &DeformationGradient, internal_variables: &V| {
300                    Ok(self
301                        .first_piola_kirchhoff_stress(deformation_gradient, internal_variables)?)
302                },
303                |deformation_gradient: &DeformationGradient, internal_variables: &V| {
304                    Ok(self.internal_variables_residual(deformation_gradient, internal_variables)?)
305                },
306                |deformation_gradient: &DeformationGradient, internal_variables: &V| {
307                    Ok(self.tangents(deformation_gradient, internal_variables)?)
308                },
309                (
310                    DeformationGradient::identity(),
311                    self.internal_variables_initial(),
312                ),
313                constraint_global,
314                constraint_local,
315                None,
316                strategy,
317            )
318            .map_err(|error| ConstitutiveError::upstream(error, self))
319    }
320}
321
322#[doc(hidden)]
323pub fn bcs_block<C, V>(
324    model: &C,
325    applied_load: AppliedLoad,
326) -> ((CscMatrix, Vector), (CscMatrix, Vector))
327where
328    C: ElasticIV<V>,
329    V: Tensor,
330{
331    let fixed = model.internal_variables_fixed();
332    let num_internal_variables = model.internal_variables_initial().size();
333    let pattern_vars = fixed.iter().enumerate().map(|(i, &j)| (i, j)).collect();
334    let mut matrix_vars =
335        CscMatrix::from_pattern(fixed.len(), num_internal_variables, pattern_vars);
336    matrix_vars.fill(|_, _| 1.0);
337    let local = (matrix_vars, Vector::zero(fixed.len()));
338    let (vector, pattern) = match applied_load {
339        AppliedLoad::UniaxialStress(deformation_gradient_11) => {
340            let mut vector = Vector::zero(4);
341            vector[0] = deformation_gradient_11;
342            (vector, vec![(0, 0), (1, 1), (2, 2), (3, 5)])
343        }
344        AppliedLoad::BiaxialStress(deformation_gradient_11, deformation_gradient_22) => {
345            let mut vector = Vector::zero(5);
346            vector[0] = deformation_gradient_11;
347            vector[4] = deformation_gradient_22;
348            (vector, vec![(0, 0), (1, 1), (2, 2), (3, 5), (4, 4)])
349        }
350    };
351    let mut matrix = CscMatrix::from_pattern(vector.len(), 9, pattern);
352    matrix.fill(|_, _| 1.0);
353    ((matrix, vector), local)
354}
355
356#[doc(hidden)]
357pub fn bcs<C, V>(model: &C, applied_load: AppliedLoad) -> (Matrix, Vector)
358where
359    C: ElasticIV<V>,
360    V: Tensor,
361{
362    let fixed = model.internal_variables_fixed();
363    let num_internal_variables = model.internal_variables_initial().size();
364    let num_deformation_gradient = 9;
365    let (num_constraints, prescribed) = match applied_load {
366        AppliedLoad::UniaxialStress(deformation_gradient_11) => (
367            4,
368            vec![
369                (0, 0, deformation_gradient_11),
370                (1, 1, 0.0),
371                (2, 2, 0.0),
372                (3, 5, 0.0),
373            ],
374        ),
375        AppliedLoad::BiaxialStress(deformation_gradient_11, deformation_gradient_22) => (
376            5,
377            vec![
378                (0, 0, deformation_gradient_11),
379                (1, 1, 0.0),
380                (2, 2, 0.0),
381                (3, 5, 0.0),
382                (4, 4, deformation_gradient_22),
383            ],
384        ),
385    };
386    let mut matrix = Matrix::zero(
387        num_constraints + fixed.len(),
388        num_deformation_gradient + num_internal_variables,
389    );
390    let mut vector = Vector::zero(num_constraints + fixed.len());
391    prescribed.iter().for_each(|&(row, column, value)| {
392        matrix[row][column] = 1.0;
393        vector[row] = value
394    });
395    fixed
396        .iter()
397        .enumerate()
398        .for_each(|(i, &j)| matrix[num_constraints + i][num_deformation_gradient + j] = 1.0);
399    (matrix, vector)
400}
conspire/constitutive/solid/elastic/internal_variables/mod.rs

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