conspire/domain/fem/block/element/solid/elastic/
mod.rs

1use crate::{
2    constitutive::solid::elastic::Elastic,
3    fem::block::element::{
4        Element, ElementNodalCoordinates, FiniteElement, FiniteElementError, GradientVectors,
5        solid::{ElementNodalForcesSolid, ElementNodalStiffnessesSolid, SolidFiniteElement},
6        surface::{SurfaceElement, SurfaceFiniteElement},
7    },
8    math::{ContractSecondFourthIndicesWithFirstIndicesOf, IDENTITY, Scalar, Tensor},
9    mechanics::{FirstPiolaKirchhoffStressList, FirstPiolaKirchhoffTangentStiffnessList},
10};
11
12pub trait ElasticFiniteElement<C, const G: usize, const M: usize, const N: usize, const P: usize>
13where
14    C: Elastic,
15    Self: SolidFiniteElement<G, M, N, P>,
16{
17    fn nodal_forces(
18        &self,
19        constitutive_model: &C,
20        nodal_coordinates: &ElementNodalCoordinates<N>,
21    ) -> Result<ElementNodalForcesSolid<N>, FiniteElementError>;
22    fn nodal_stiffnesses(
23        &self,
24        constitutive_model: &C,
25        nodal_coordinates: &ElementNodalCoordinates<N>,
26    ) -> Result<ElementNodalStiffnessesSolid<N>, FiniteElementError>;
27}
28
29impl<C, const G: usize, const N: usize, const O: usize, const P: usize>
30    ElasticFiniteElement<C, G, 3, N, P> for Element<G, N, O>
31where
32    C: Elastic,
33    Self: SolidFiniteElement<G, 3, N, P>,
34{
35    fn nodal_forces(
36        &self,
37        constitutive_model: &C,
38        nodal_coordinates: &ElementNodalCoordinates<N>,
39    ) -> Result<ElementNodalForcesSolid<N>, FiniteElementError> {
40        nodal_forces::<_, _, _, _, _, O, _>(
41            self,
42            constitutive_model,
43            self.gradient_vectors(),
44            nodal_coordinates,
45        )
46    }
47    fn nodal_stiffnesses(
48        &self,
49        constitutive_model: &C,
50        nodal_coordinates: &ElementNodalCoordinates<N>,
51    ) -> Result<ElementNodalStiffnessesSolid<N>, FiniteElementError> {
52        match self
53            .deformation_gradients(nodal_coordinates)
54            .iter()
55            .map(|deformation_gradient| {
56                constitutive_model.first_piola_kirchhoff_tangent_stiffness(deformation_gradient)
57            })
58            .collect::<Result<FirstPiolaKirchhoffTangentStiffnessList<G>, _>>()
59        {
60            Ok(first_piola_kirchhoff_tangent_stiffnesses) => {
61                Ok(first_piola_kirchhoff_tangent_stiffnesses
62                    .iter()
63                    .zip(
64                        self.gradient_vectors()
65                            .iter()
66                            .zip(self.integration_weights()),
67                    )
68                    .map(
69                        |(
70                            first_piola_kirchhoff_tangent_stiffness,
71                            (gradient_vectors, integration_weight),
72                        )| {
73                            gradient_vectors
74                                .iter()
75                                .map(|gradient_vector_a| {
76                                    gradient_vectors
77                                        .iter()
78                                        .map(|gradient_vector_b| {
79                                            first_piola_kirchhoff_tangent_stiffness
80                                            .contract_second_fourth_indices_with_first_indices_of(
81                                                gradient_vector_a,
82                                                gradient_vector_b,
83                                            )
84                                            * integration_weight
85                                        })
86                                        .collect()
87                                })
88                                .collect()
89                        },
90                    )
91                    .sum())
92            }
93            Err(error) => Err(FiniteElementError::Upstream(
94                format!("{error}"),
95                format!("{self:?}"),
96            )),
97        }
98    }
99}
100
101impl<C, const G: usize, const N: usize, const O: usize> ElasticFiniteElement<C, G, 2, N, N>
102    for SurfaceElement<G, N, O>
103where
104    C: Elastic,
105    Self: SolidFiniteElement<G, 2, N, N>,
106{
107    fn nodal_forces(
108        &self,
109        constitutive_model: &C,
110        nodal_coordinates: &ElementNodalCoordinates<N>,
111    ) -> Result<ElementNodalForcesSolid<N>, FiniteElementError> {
112        nodal_forces::<_, _, _, _, _, O, _>(
113            self,
114            constitutive_model,
115            self.gradient_vectors(),
116            nodal_coordinates,
117        )
118    }
119    fn nodal_stiffnesses(
120        &self,
121        constitutive_model: &C,
122        nodal_coordinates: &ElementNodalCoordinates<N>,
123    ) -> Result<ElementNodalStiffnessesSolid<N>, FiniteElementError> {
124        match self.deformation_gradients(nodal_coordinates).iter()
125            .map(|deformation_gradient| {
126                constitutive_model.first_piola_kirchhoff_tangent_stiffness(deformation_gradient)
127            })
128            .collect::<Result<FirstPiolaKirchhoffTangentStiffnessList<G>, _>>()
129            {
130            Ok(first_piola_kirchhoff_tangent_stiffnesses) => Ok(first_piola_kirchhoff_tangent_stiffnesses
131            .iter()
132            .zip(
133                self.gradient_vectors()
134                    .iter()
135                    .zip(self.integration_weights().iter()
136                    .zip(self.reference_normals().iter()
137                    .zip(Self::normal_gradients(nodal_coordinates))
138                )
139                ),
140            )
141            .map(
142                |(
143                    first_piola_kirchhoff_tangent_stiffness,
144                    (gradient_vectors, (integration_weight, (reference_normal, normal_gradients))),
145                )| {
146                    gradient_vectors.iter()
147                    .map(|gradient_vector_a|
148                        gradient_vectors.iter()
149                        .zip(normal_gradients.iter())
150                        .map(|(gradient_vector_b, normal_gradient_b)|
151                            first_piola_kirchhoff_tangent_stiffness.iter()
152                            .map(|first_piola_kirchhoff_tangent_stiffness_m|
153                                IDENTITY.iter()
154                                .zip(normal_gradient_b.iter())
155                                .map(|(identity_n, normal_gradient_b_n)|
156                                    first_piola_kirchhoff_tangent_stiffness_m.iter()
157                                    .zip(gradient_vector_a.iter())
158                                    .map(|(first_piola_kirchhoff_tangent_stiffness_mj, gradient_vector_a_j)|
159                                        first_piola_kirchhoff_tangent_stiffness_mj.iter()
160                                        .zip(identity_n.iter()
161                                        .zip(normal_gradient_b_n.iter()))
162                                        .map(|(first_piola_kirchhoff_tangent_stiffness_mjk, (identity_nk, normal_gradient_b_n_k))|
163                                            first_piola_kirchhoff_tangent_stiffness_mjk.iter()
164                                            .zip(gradient_vector_b.iter()
165                                            .zip(reference_normal.iter()))
166                                            .map(|(first_piola_kirchhoff_tangent_stiffness_mjkl, (gradient_vector_b_l, reference_normal_l))|
167                                                first_piola_kirchhoff_tangent_stiffness_mjkl * gradient_vector_a_j * (
168                                                    identity_nk * gradient_vector_b_l + normal_gradient_b_n_k * reference_normal_l
169                                                ) * integration_weight
170                                            ).sum::<Scalar>()
171                                        ).sum::<Scalar>()
172                                    ).sum::<Scalar>()
173                                ).collect()
174                            ).collect()
175                        ).collect()
176                    ).collect()
177                }
178            )
179            .sum()),
180            Err(error) => Err(FiniteElementError::Upstream(
181                format!("{error}"),
182                format!("{self:?}"),
183            )),
184        }
185    }
186}
187
188fn nodal_forces<
189    C,
190    F,
191    const G: usize,
192    const M: usize,
193    const N: usize,
194    const O: usize,
195    const P: usize,
196>(
197    element: &F,
198    constitutive_model: &C,
199    gradient_vectors: &GradientVectors<G, N>,
200    nodal_coordinates: &ElementNodalCoordinates<N>,
201) -> Result<ElementNodalForcesSolid<N>, FiniteElementError>
202where
203    C: Elastic,
204    F: SolidFiniteElement<G, M, N, P>,
205{
206    match element
207        .deformation_gradients(nodal_coordinates)
208        .iter()
209        .map(|deformation_gradient| {
210            constitutive_model.first_piola_kirchhoff_stress(deformation_gradient)
211        })
212        .collect::<Result<FirstPiolaKirchhoffStressList<G>, _>>()
213    {
214        Ok(first_piola_kirchhoff_stresses) => Ok(first_piola_kirchhoff_stresses
215            .iter()
216            .zip(gradient_vectors.iter().zip(element.integration_weights()))
217            .map(
218                |(first_piola_kirchhoff_stress, (gradient_vectors, integration_weight))| {
219                    gradient_vectors
220                        .iter()
221                        .map(|gradient_vector| {
222                            (first_piola_kirchhoff_stress * gradient_vector) * integration_weight
223                        })
224                        .collect()
225                },
226            )
227            .sum()),
228        Err(error) => Err(FiniteElementError::Upstream(
229            format!("{error}"),
230            format!("{element:?}"),
231        )),
232    }
233}
conspire/domain/fem/block/element/solid/elastic/mod.rs

conspire/domain/fem/block/element/solid/elastic/
mod.rs
