conspire/constitutive/hybrid/elastic/multiplicative/
mod.rs

1#[cfg(test)]
2mod test;
3
4use crate::{
5    constitutive::{
6        ConstitutiveError,
7        hybrid::{Multiplicative, MultiplicativeTrait},
8        solid::{
9            Solid,
10            elastic::{Elastic, internal_variables::ElasticIV},
11        },
12    },
13    math::{
14        ContractThirdIndexWithFirstIndexOf, IDENTITY_10, Rank2, TensorArray, TensorRank2,
15        TensorRank4,
16        optimize::{EqualityConstraint, GradientDescent, ZerothOrderRootFinding},
17    },
18    mechanics::{
19        CauchyStress, CauchyTangentStiffness, CauchyTangentStiffness1, DeformationGradient,
20        DeformationGradient2, FirstPiolaKirchhoffStress, FirstPiolaKirchhoffStress1,
21        FirstPiolaKirchhoffStress2, FirstPiolaKirchhoffTangentStiffness,
22        FirstPiolaKirchhoffTangentStiffness2, Scalar, SecondPiolaKirchhoffStress,
23        SecondPiolaKirchhoffTangentStiffness,
24    },
25};
26
27impl<C1, C2> Solid for Multiplicative<C1, C2>
28where
29    // C1: Elastic,
30    // C2: Elastic,
31    C1: Solid,
32    C2: Clone + std::fmt::Debug,
33{
34    fn bulk_modulus(&self) -> Scalar {
35        todo!() // can do right when switch to the _foo methods
36    }
37    fn shear_modulus(&self) -> Scalar {
38        todo!() // can do right when switch to the _foo methods
39    }
40}
41
42impl<C1, C2>
43    ElasticIV<
44        DeformationGradient2,
45        TensorRank4<3, 2, 0, 1, 0>,
46        TensorRank4<3, 1, 0, 2, 0>,
47        FirstPiolaKirchhoffTangentStiffness2,
48    > for Multiplicative<C1, C2>
49where
50    C1: Elastic,
51    C2: Elastic,
52{
53    /// Calculates and returns the Cauchy stress.
54    ///
55    /// ```math
56    /// \boldsymbol{\sigma} = \frac{1}{J_2}\,\boldsymbol{\sigma}_1
57    /// ```
58    fn cauchy_stress_foo(
59        &self,
60        deformation_gradient: &DeformationGradient,
61        deformation_gradient_2: &DeformationGradient2,
62    ) -> Result<CauchyStress, ConstitutiveError> {
63        let (deformation_gradient_2_inverse, jacobian_2) =
64            deformation_gradient_2.inverse_and_determinant();
65        let deformation_gradient_1 = deformation_gradient * &deformation_gradient_2_inverse;
66        Ok(self.0.cauchy_stress(&deformation_gradient_1.into())? / jacobian_2)
67    }
68    /// Calculates and returns the tangent stiffness associated with the Cauchy stress.
69    ///
70    /// ```math
71    /// \boldsymbol{\mathcal{T}} = \frac{1}{J_2}\,\boldsymbol{\mathcal{T}}_1\cdot\mathbf{F}_2^{-T}
72    /// ```
73    fn cauchy_tangent_stiffness_foo(
74        &self,
75        deformation_gradient: &DeformationGradient,
76        deformation_gradient_2: &DeformationGradient2,
77    ) -> Result<CauchyTangentStiffness, ConstitutiveError> {
78        let (deformation_gradient_2_inverse, jacobian_2) =
79            deformation_gradient_2.inverse_and_determinant();
80        let deformation_gradient_1 = deformation_gradient * &deformation_gradient_2_inverse;
81        Ok(CauchyTangentStiffness1::from(
82            self.0
83                .cauchy_tangent_stiffness(&deformation_gradient_1.into())?,
84        ) * (deformation_gradient_2_inverse.transpose() / jacobian_2))
85    }
86    /// Calculates and returns the first Piola-Kirchhoff stress.
87    ///
88    /// ```math
89    /// \mathbf{P} = \mathbf{P}_1\cdot\mathbf{F}_2^{-T}
90    /// ```
91    fn first_piola_kirchhoff_stress_foo(
92        &self,
93        deformation_gradient: &DeformationGradient,
94        deformation_gradient_2: &DeformationGradient2,
95    ) -> Result<FirstPiolaKirchhoffStress, ConstitutiveError> {
96        Ok(
97            self.cauchy_stress_foo(deformation_gradient, deformation_gradient_2)?
98                * deformation_gradient.inverse_transpose()
99                * deformation_gradient.determinant(),
100        )
101    }
102    /// Calculates and returns the second Piola-Kirchhoff stress.
103    ///
104    /// ```math
105    /// \mathbf{S} = \mathbf{F}_2^{-1}\cdot\mathbf{S}_1\cdot\mathbf{F}_2^{-T}
106    /// ```
107    fn second_piola_kirchhoff_stress_foo(
108        &self,
109        deformation_gradient: &DeformationGradient,
110        deformation_gradient_2: &DeformationGradient2,
111    ) -> Result<SecondPiolaKirchhoffStress, ConstitutiveError> {
112        Ok(deformation_gradient.inverse()
113            * self
114                .first_piola_kirchhoff_stress_foo(deformation_gradient, deformation_gradient_2)?)
115    }
116    fn internal_variables_initial(&self) -> DeformationGradient2 {
117        DeformationGradient2::identity()
118    }
119    /// Calculates and returns the residual associated with the second deformation gradient.
120    ///
121    /// ```math
122    /// \mathbf{R} = \mathbf{P}_2 - \mathbf{M}_1\cdot\mathbf{F}_2^{-T}
123    /// ```
124    fn internal_variables_residual(
125        &self,
126        deformation_gradient: &DeformationGradient,
127        deformation_gradient_2: &DeformationGradient2,
128    ) -> Result<DeformationGradient2, ConstitutiveError> {
129        let deformation_gradient_2_inverse = deformation_gradient_2.inverse();
130        let deformation_gradient_1 = deformation_gradient * &deformation_gradient_2_inverse;
131        Ok(FirstPiolaKirchhoffStress2::from(
132            self.1
133                .first_piola_kirchhoff_stress(deformation_gradient_2.into())?,
134        ) - deformation_gradient_1.transpose()
135            * FirstPiolaKirchhoffStress1::from(
136                self.0
137                    .first_piola_kirchhoff_stress(&deformation_gradient_1.into())?,
138            )
139            * deformation_gradient_2_inverse.transpose())
140    }
141    /// Calculates and returns the tangents associated with the internal variables.
142    ///
143    /// ```math
144    /// \frac{\partial P_{iJ}}{\partial F_{KL}^2} = -P_{iL}F_{KJ}^{2-T} - \mathcal{C}_{iJmL}F_{mK}^1
145    /// ```
146    /// ```math
147    /// \frac{\partial R_{IJ}}{\partial F_{kL}} = -F_{IL}^{2-T}P_{kJ} - F_{mI}^1\mathcal{C}_{mJkL}
148    /// ```
149    /// ```math
150    /// \frac{\partial R_{IJ}}{\partial F_{KL}^2} = \mathcal{C}_{IJKL}^2 + F_{IM}^1P_{ML}{F_{KJ}^{2-T}} - \frac{\partial R_{IJ}}{\partial F_{mL}}\,F_{mK}^1
151    /// ```
152    fn internal_variables_tangents(
153        &self,
154        deformation_gradient: &DeformationGradient,
155        deformation_gradient_2: &DeformationGradient2,
156    ) -> Result<
157        (
158            TensorRank4<3, 2, 0, 1, 0>,
159            TensorRank4<3, 1, 0, 2, 0>,
160            FirstPiolaKirchhoffTangentStiffness2,
161        ),
162        ConstitutiveError,
163    > {
164        //
165        // If hyperelastic, tangent_1 should equal a transpose of tangent_2.
166        // Could add a method to utilize that for minimize() in hyperelastic/multiplicative.
167        //
168        let deformation_gradient_2_inverse = deformation_gradient_2.inverse();
169        let deformation_gradient_2_inverse_transpose = deformation_gradient_2_inverse.transpose();
170        let deformation_gradient_1 = deformation_gradient * &deformation_gradient_2_inverse;
171        let deformation_gradient_1_transpose = deformation_gradient_1.transpose();
172        let first_piola_kirchhoff_stress =
173            self.first_piola_kirchhoff_stress_foo(deformation_gradient, deformation_gradient_2)?;
174        let tangent_0 = self.first_piola_kirchhoff_tangent_stiffness_foo(
175            deformation_gradient,
176            deformation_gradient_2,
177        )?;
178        let tangent_1 = TensorRank4::dyad_il_kj(
179            &(deformation_gradient_2_inverse_transpose * -1.0),
180            &first_piola_kirchhoff_stress,
181        ) - &deformation_gradient_1_transpose * &tangent_0;
182        let tangent_2 = TensorRank4::dyad_il_jk(
183            &first_piola_kirchhoff_stress,
184            &(&deformation_gradient_2_inverse * -1.0),
185        ) - tangent_0
186            .contract_third_index_with_first_index_of(&deformation_gradient_1);
187        let tangent_3 = FirstPiolaKirchhoffTangentStiffness2::from(
188            self.1
189                .first_piola_kirchhoff_tangent_stiffness(deformation_gradient_2.into())?,
190        ) - tangent_1
191            .contract_third_index_with_first_index_of(&deformation_gradient_1)
192            + TensorRank4::dyad_il_jk(
193                &(deformation_gradient_1_transpose * first_piola_kirchhoff_stress),
194                &deformation_gradient_2_inverse,
195            );
196        Ok((tangent_1, tangent_2, tangent_3))
197    }
198    fn internal_variables_constraints(&self) -> (&[usize], usize) {
199        (&[10, 11, 14], 9)
200    }
201}
202
203impl<C1, C2> Elastic for Multiplicative<C1, C2>
204where
205    C1: Elastic,
206    C2: Elastic,
207{
208    /// Calculates and returns the Cauchy stress.
209    ///
210    /// ```math
211    /// \boldsymbol{\sigma}(\mathbf{F}) = \frac{1}{J_2}\,\boldsymbol{\sigma}_1(\mathbf{F}_1)
212    /// ```
213    fn cauchy_stress(
214        &self,
215        deformation_gradient: &DeformationGradient,
216    ) -> Result<CauchyStress, ConstitutiveError> {
217        let (deformation_gradient_1, deformation_gradient_2) =
218            self.deformation_gradients(deformation_gradient)?;
219        Ok(self.0.cauchy_stress(&deformation_gradient_1)? / deformation_gradient_2.determinant())
220    }
221    /// Dummy method that will panic.
222    fn cauchy_tangent_stiffness(
223        &self,
224        _: &DeformationGradient,
225    ) -> Result<CauchyTangentStiffness, ConstitutiveError> {
226        panic!()
227    }
228    /// Calculates and returns the first Piola-Kirchhoff stress.
229    ///
230    /// ```math
231    /// \mathbf{P}(\mathbf{F}) = \mathbf{P}_1(\mathbf{F}_1)\cdot\mathbf{F}_2^{-T}
232    /// ```
233    fn first_piola_kirchhoff_stress(
234        &self,
235        deformation_gradient: &DeformationGradient,
236    ) -> Result<FirstPiolaKirchhoffStress, ConstitutiveError> {
237        let (deformation_gradient_1, deformation_gradient_2) =
238            self.deformation_gradients(deformation_gradient)?;
239        let deformation_gradient_2_inverse_transpose: TensorRank2<3, 0, 0> =
240            deformation_gradient_2.inverse_transpose().into();
241        Ok(self
242            .0
243            .first_piola_kirchhoff_stress(&deformation_gradient_1)?
244            * deformation_gradient_2_inverse_transpose)
245    }
246    /// Dummy method that will panic.
247    fn first_piola_kirchhoff_tangent_stiffness(
248        &self,
249        _: &DeformationGradient,
250    ) -> Result<FirstPiolaKirchhoffTangentStiffness, ConstitutiveError> {
251        panic!()
252    }
253    /// Calculates and returns the second Piola-Kirchhoff stress.
254    ///
255    /// ```math
256    /// \mathbf{S}(\mathbf{F}) = \mathbf{F}_2^{-1}\cdot\mathbf{S}_1(\mathbf{F}_1)\cdot\mathbf{F}_2^{-T}
257    /// ```
258    fn second_piola_kirchhoff_stress(
259        &self,
260        deformation_gradient: &DeformationGradient,
261    ) -> Result<SecondPiolaKirchhoffStress, ConstitutiveError> {
262        let (deformation_gradient_1, deformation_gradient_2) =
263            self.deformation_gradients(deformation_gradient)?;
264        let deformation_gradient_2_inverse: TensorRank2<3, 0, 0> =
265            deformation_gradient_2.inverse().into();
266        Ok(&deformation_gradient_2_inverse
267            * self
268                .0
269                .second_piola_kirchhoff_stress(&deformation_gradient_1)?
270            * deformation_gradient_2_inverse.transpose())
271    }
272    /// Dummy method that will panic.
273    fn second_piola_kirchhoff_tangent_stiffness(
274        &self,
275        _: &DeformationGradient,
276    ) -> Result<SecondPiolaKirchhoffTangentStiffness, ConstitutiveError> {
277        panic!()
278    }
279}
280
281impl<C1, C2> MultiplicativeTrait for Multiplicative<C1, C2>
282where
283    C1: Elastic,
284    C2: Elastic,
285{
286    fn deformation_gradients(
287        &self,
288        deformation_gradient: &DeformationGradient,
289    ) -> Result<(DeformationGradient, DeformationGradient), ConstitutiveError> {
290        if deformation_gradient.is_identity() {
291            Ok((IDENTITY_10, IDENTITY_10))
292        } else {
293            match GradientDescent::default().root(
294                |deformation_gradient_2: &DeformationGradient| {
295                    let deformation_gradient_1: DeformationGradient =
296                        (deformation_gradient * deformation_gradient_2.inverse()).into();
297                    let deformation_gradient_2_inverse_transpose: TensorRank2<3, 0, 0> =
298                        deformation_gradient_2.inverse_transpose().into();
299                    let right_hand_side: FirstPiolaKirchhoffStress = (deformation_gradient_1
300                        .transpose()
301                        * self
302                            .0
303                            .first_piola_kirchhoff_stress(&deformation_gradient_1)?
304                        * deformation_gradient_2_inverse_transpose)
305                        .into();
306                    Ok(self
307                        .1
308                        .first_piola_kirchhoff_stress(deformation_gradient_2)?
309                        - right_hand_side)
310                },
311                IDENTITY_10,
312                EqualityConstraint::None,
313            ) {
314                Ok(deformation_gradient_2) => {
315                    let deformation_gradient_1 =
316                        (deformation_gradient * deformation_gradient_2.inverse()).into();
317                    Ok((deformation_gradient_1, deformation_gradient_2))
318                }
319                Err(error) => Err(ConstitutiveError::Upstream(
320                    format!("{error}"),
321                    format!("{self:?}"),
322                )),
323            }
324        }
325    }
326}
conspire/constitutive/hybrid/elastic/multiplicative/mod.rs

conspire/constitutive/hybrid/elastic/multiplicative/
mod.rs
