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

1#[cfg(test)]
2mod test;
3
4use crate::{
5    ABS_TOL,
6    constitutive::{
7        ConstitutiveError,
8        hybrid::{Hybrid, Multiplicative, MultiplicativeTrait},
9        solid::{Solid, elastic::Elastic},
10    },
11    math::{IDENTITY_10, Rank2, Tensor, TensorRank2, ZERO_10},
12    mechanics::{
13        CauchyStress, CauchyTangentStiffness, DeformationGradient, FirstPiolaKirchhoffStress,
14        FirstPiolaKirchhoffTangentStiffness, Scalar, SecondPiolaKirchhoffStress,
15        SecondPiolaKirchhoffTangentStiffness,
16    },
17};
18
19impl<C1, C2> Solid for Multiplicative<C1, C2> {
20    /// Dummy method that will panic.
21    fn bulk_modulus(&self) -> &Scalar {
22        panic!()
23    }
24    /// Dummy method that will panic.
25    fn shear_modulus(&self) -> &Scalar {
26        panic!()
27    }
28}
29
30impl<C1, C2> Elastic for Multiplicative<C1, C2>
31where
32    C1: Elastic,
33    C2: Elastic,
34{
35    /// Calculates and returns the Cauchy stress.
36    ///
37    /// ```math
38    /// \boldsymbol{\sigma}(\mathbf{F}) = \frac{1}{J_2}\,\boldsymbol{\sigma}_1(\mathbf{F}_1)
39    /// ```
40    fn cauchy_stress(
41        &self,
42        deformation_gradient: &DeformationGradient,
43    ) -> Result<CauchyStress, ConstitutiveError> {
44        let (deformation_gradient_1, deformation_gradient_2) =
45            self.deformation_gradients(deformation_gradient)?;
46        Ok(self
47            .constitutive_model_1()
48            .cauchy_stress(&deformation_gradient_1)?
49            / deformation_gradient_2.determinant())
50    }
51    /// Dummy method that will panic.
52    fn cauchy_tangent_stiffness(
53        &self,
54        _: &DeformationGradient,
55    ) -> Result<CauchyTangentStiffness, ConstitutiveError> {
56        panic!()
57    }
58    /// Calculates and returns the first Piola-Kirchhoff stress.
59    ///
60    /// ```math
61    /// \mathbf{P}(\mathbf{F}) = \mathbf{P}_1(\mathbf{F}_1)\cdot\mathbf{F}_2^{-T}
62    /// ```
63    fn first_piola_kirchhoff_stress(
64        &self,
65        deformation_gradient: &DeformationGradient,
66    ) -> Result<FirstPiolaKirchhoffStress, ConstitutiveError> {
67        let (deformation_gradient_1, deformation_gradient_2) =
68            self.deformation_gradients(deformation_gradient)?;
69        let deformation_gradient_2_inverse_transpose: TensorRank2<3, 0, 0> =
70            deformation_gradient_2.inverse_transpose().into();
71        Ok(self
72            .constitutive_model_1()
73            .first_piola_kirchhoff_stress(&deformation_gradient_1)?
74            * deformation_gradient_2_inverse_transpose)
75    }
76    /// Dummy method that will panic.
77    fn first_piola_kirchhoff_tangent_stiffness(
78        &self,
79        _: &DeformationGradient,
80    ) -> Result<FirstPiolaKirchhoffTangentStiffness, ConstitutiveError> {
81        panic!()
82    }
83    /// Calculates and returns the second Piola-Kirchhoff stress.
84    ///
85    /// ```math
86    /// \mathbf{S}(\mathbf{F}) = \mathbf{F}_2^{-1}\cdot\mathbf{S}_1(\mathbf{F}_1)\cdot\mathbf{F}_2^{-T}
87    /// ```
88    fn second_piola_kirchhoff_stress(
89        &self,
90        deformation_gradient: &DeformationGradient,
91    ) -> Result<SecondPiolaKirchhoffStress, ConstitutiveError> {
92        let (deformation_gradient_1, deformation_gradient_2) =
93            self.deformation_gradients(deformation_gradient)?;
94        let deformation_gradient_2_inverse: TensorRank2<3, 0, 0> =
95            deformation_gradient_2.inverse().into();
96        Ok(&deformation_gradient_2_inverse
97            * self
98                .constitutive_model_1()
99                .second_piola_kirchhoff_stress(&deformation_gradient_1)?
100            * deformation_gradient_2_inverse.transpose())
101    }
102    /// Dummy method that will panic.
103    fn second_piola_kirchhoff_tangent_stiffness(
104        &self,
105        _: &DeformationGradient,
106    ) -> Result<SecondPiolaKirchhoffTangentStiffness, ConstitutiveError> {
107        panic!()
108    }
109}
110
111impl<C1, C2> MultiplicativeTrait for Multiplicative<C1, C2>
112where
113    C1: Elastic,
114    C2: Elastic,
115{
116    fn deformation_gradients(
117        &self,
118        deformation_gradient: &DeformationGradient,
119    ) -> Result<(DeformationGradient, DeformationGradient), ConstitutiveError> {
120        if deformation_gradient.is_identity() {
121            Ok((IDENTITY_10, IDENTITY_10))
122        } else {
123            let mut deformation_gradient_1 = IDENTITY_10;
124            let mut deformation_gradient_2 = deformation_gradient.clone();
125            let mut deformation_gradient_2_old = IDENTITY_10;
126            let mut deformation_gradient_2_inverse_transpose: TensorRank2<3, 0, 0>;
127            let mut residual: FirstPiolaKirchhoffStress;
128            let mut residual_increment: FirstPiolaKirchhoffStress;
129            let mut residual_norm = 1.0;
130            let mut residual_old = ZERO_10;
131            let mut right_hand_side: FirstPiolaKirchhoffStress;
132            let mut steps: u8 = 0;
133            let steps_maximum: u8 = 50;
134            let mut step_size: Scalar;
135            while steps < steps_maximum {
136                deformation_gradient_1 =
137                    (deformation_gradient * deformation_gradient_2.inverse()).into();
138                deformation_gradient_2_inverse_transpose =
139                    deformation_gradient_2.inverse_transpose().into();
140                right_hand_side = (deformation_gradient_1.transpose()
141                    * self
142                        .constitutive_model_1()
143                        .first_piola_kirchhoff_stress(&deformation_gradient_1)?
144                    * deformation_gradient_2_inverse_transpose)
145                    .into();
146                residual = self
147                    .constitutive_model_2()
148                    .first_piola_kirchhoff_stress(&deformation_gradient_2)?
149                    - right_hand_side;
150                residual_norm = residual.norm();
151                if residual_norm >= ABS_TOL {
152                    residual_increment = residual_old - &residual;
153                    step_size = (deformation_gradient_2_old - &deformation_gradient_2)
154                        .full_contraction(&residual_increment)
155                        .abs()
156                        / residual_increment.norm_squared();
157                    deformation_gradient_2_old = deformation_gradient_2.clone();
158                    residual_old = residual.clone();
159                    deformation_gradient_2 -= residual * step_size;
160                } else {
161                    break;
162                }
163                steps += 1;
164            }
165            if residual_norm >= ABS_TOL && steps == steps_maximum {
166                panic!("MAX STEPS REACHED")
167            } else {
168                Ok((deformation_gradient_1, deformation_gradient_2))
169            }
170        }
171    }
172}
conspire/constitutive/hybrid/elastic/multiplicative/mod.rs

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