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        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(deformation_gradient, internal_variables)?,
258                    self.internal_variables_residual(deformation_gradient, internal_variables)?,
259                )))
260            },
261            Self::Variables::from((
262                DeformationGradient::identity(),
263                self.internal_variables_initial(),
264            )),
265            EqualityConstraint::Linear(matrix, vector),
266        ) {
267            Ok(solution) => Ok(solution.into()),
268            Err(error) => Err(ConstitutiveError::Upstream(
269                format!("{error}"),
270                format!("{self:?}"),
271            )),
272        }
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        match solver.root_block(
298            |deformation_gradient: &DeformationGradient, internal_variables: &V| {
299                Ok(self.first_piola_kirchhoff_stress(deformation_gradient, internal_variables)?)
300            },
301            |deformation_gradient: &DeformationGradient, internal_variables: &V| {
302                Ok(self.internal_variables_residual(deformation_gradient, internal_variables)?)
303            },
304            |deformation_gradient: &DeformationGradient, internal_variables: &V| {
305                Ok(self.tangents(deformation_gradient, internal_variables)?)
306            },
307            (
308                DeformationGradient::identity(),
309                self.internal_variables_initial(),
310            ),
311            constraint_global,
312            constraint_local,
313            None,
314            strategy,
315        ) {
316            Ok(solution) => Ok(solution),
317            Err(error) => Err(ConstitutiveError::Upstream(
318                format!("{error}"),
319                format!("{self:?}"),
320            )),
321        }
322    }
323}
324
325#[doc(hidden)]
326pub fn bcs_block<C, V>(
327    model: &C,
328    applied_load: AppliedLoad,
329) -> ((CscMatrix, Vector), (CscMatrix, Vector))
330where
331    C: ElasticIV<V>,
332    V: Tensor,
333{
334    let fixed = model.internal_variables_fixed();
335    let num_internal_variables = model.internal_variables_initial().size();
336    let pattern_vars = fixed.iter().enumerate().map(|(i, &j)| (i, j)).collect();
337    let mut matrix_vars =
338        CscMatrix::from_pattern(fixed.len(), num_internal_variables, pattern_vars);
339    matrix_vars.fill(|_, _| 1.0);
340    let local = (matrix_vars, Vector::zero(fixed.len()));
341    let (vector, pattern) = match applied_load {
342        AppliedLoad::UniaxialStress(deformation_gradient_11) => {
343            let mut vector = Vector::zero(4);
344            vector[0] = deformation_gradient_11;
345            (vector, vec![(0, 0), (1, 1), (2, 2), (3, 5)])
346        }
347        AppliedLoad::BiaxialStress(deformation_gradient_11, deformation_gradient_22) => {
348            let mut vector = Vector::zero(5);
349            vector[0] = deformation_gradient_11;
350            vector[4] = deformation_gradient_22;
351            (vector, vec![(0, 0), (1, 1), (2, 2), (3, 5), (4, 4)])
352        }
353    };
354    let mut matrix = CscMatrix::from_pattern(vector.len(), 9, pattern);
355    matrix.fill(|_, _| 1.0);
356    ((matrix, vector), local)
357}
358
359#[doc(hidden)]
360pub fn bcs<C, V>(model: &C, applied_load: AppliedLoad) -> (Matrix, Vector)
361where
362    C: ElasticIV<V>,
363    V: Tensor,
364{
365    let fixed = model.internal_variables_fixed();
366    let num_internal_variables = model.internal_variables_initial().size();
367    let num_deformation_gradient = 9;
368    let (num_constraints, prescribed) = match applied_load {
369        AppliedLoad::UniaxialStress(deformation_gradient_11) => (
370            4,
371            vec![
372                (0, 0, deformation_gradient_11),
373                (1, 1, 0.0),
374                (2, 2, 0.0),
375                (3, 5, 0.0),
376            ],
377        ),
378        AppliedLoad::BiaxialStress(deformation_gradient_11, deformation_gradient_22) => (
379            5,
380            vec![
381                (0, 0, deformation_gradient_11),
382                (1, 1, 0.0),
383                (2, 2, 0.0),
384                (3, 5, 0.0),
385                (4, 4, deformation_gradient_22),
386            ],
387        ),
388    };
389    let mut matrix = Matrix::zero(
390        num_constraints + fixed.len(),
391        num_deformation_gradient + num_internal_variables,
392    );
393    let mut vector = Vector::zero(num_constraints + fixed.len());
394    prescribed.iter().for_each(|&(row, column, value)| {
395        matrix[row][column] = 1.0;
396        vector[row] = value
397    });
398    fixed
399        .iter()
400        .enumerate()
401        .for_each(|(i, &j)| matrix[num_constraints + i][num_deformation_gradient + j] = 1.0);
402    (matrix, vector)
403}
conspire/constitutive/solid/elastic/internal_variables/mod.rs

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