1use crate::{
2 constitutive::{
3 ConstitutiveError,
4 solid::{elastic::Elastic, hyperelastic::Hyperelastic},
5 },
6 fem::block::element::{
7 Element, FiniteElement, FiniteElementError, ParametricCoordinate, ParametricCoordinates,
8 ParametricReference, ShapeFunctions, ShapeFunctionsGradients, basic_from,
9 surface::SurfaceElement,
10 },
11 math::{Scalar, ScalarList, Tensor, TensorArray, TensorRank1List, TensorRank2List2D},
12 mechanics::{
13 DeformationGradient, DeformationGradientList, FirstPiolaKirchhoffStressList,
14 FirstPiolaKirchhoffTangentStiffnessList,
15 },
16};
17
18const M: usize = 2;
19
20pub type Quadrilateral = Element<2, 4, 4, 1>;
21pub type Triangle = Element<2, 1, 3, 1>;
22
23pub type PlanarElementNodalCoordinates<const N: usize> = TensorRank1List<M, 1, N>;
24pub type PlanarElementNodalReferenceCoordinates<const N: usize> = TensorRank1List<M, 0, N>;
25pub type PlanarElementNodalForcesSolid<const N: usize> = TensorRank1List<M, 1, N>;
26pub type PlanarElementNodalStiffnessesSolid<const N: usize> = TensorRank2List2D<M, 1, 1, N, N>;
27
28impl<const G: usize, const N: usize, const O: usize, const P: usize> FiniteElement<G, M, N, P>
29 for Element<2, G, N, O>
30where
31 SurfaceElement<G, N, O>: FiniteElement<G, M, N, P>,
32{
33 fn integration_points() -> ParametricCoordinates<G, M> {
34 SurfaceElement::<G, N, O>::integration_points()
35 }
36 fn integration_weights(&self) -> &ScalarList<G> {
37 &self.integration_weights
38 }
39 fn parametric_reference() -> ParametricReference<M, N> {
40 SurfaceElement::<G, N, O>::parametric_reference()
41 }
42 fn parametric_weights() -> ScalarList<G> {
43 SurfaceElement::<G, N, O>::parametric_weights()
44 }
45 fn shape_functions(parametric_coordinate: ParametricCoordinate<M>) -> ShapeFunctions<P> {
46 SurfaceElement::<G, N, O>::shape_functions(parametric_coordinate)
47 }
48 fn shape_functions_gradients(
49 parametric_coordinate: ParametricCoordinate<M>,
50 ) -> ShapeFunctionsGradients<M, P> {
51 SurfaceElement::<G, N, O>::shape_functions_gradients(parametric_coordinate)
52 }
53}
54
55impl<const G: usize, const N: usize, const O: usize> From<PlanarElementNodalReferenceCoordinates<N>>
56 for Element<2, G, N, O>
57where
58 Self: FiniteElement<G, M, N, N>,
59{
60 fn from(reference_nodal_coordinates: PlanarElementNodalReferenceCoordinates<N>) -> Self {
61 basic_from(reference_nodal_coordinates)
62 }
63}
64
65pub trait PlanarSolidFiniteElement<const G: usize, const N: usize, const P: usize>
66where
67 Self: FiniteElement<G, M, N, P>,
68{
69 fn deformation_gradients(
70 &self,
71 nodal_coordinates: &PlanarElementNodalCoordinates<N>,
72 ) -> DeformationGradientList<G>;
73}
74
75impl<const G: usize, const N: usize, const O: usize, const P: usize>
76 PlanarSolidFiniteElement<G, N, P> for Element<2, G, N, O>
77where
78 Self: FiniteElement<G, M, N, P>,
79{
80 fn deformation_gradients(
81 &self,
82 nodal_coordinates: &PlanarElementNodalCoordinates<N>,
83 ) -> DeformationGradientList<G> {
84 self.gradient_vectors()
85 .iter()
86 .map(|gradient_vectors| {
87 let mut deformation_gradient = DeformationGradient::zero();
88 deformation_gradient[2][2] = 1.0;
89 nodal_coordinates.iter().zip(gradient_vectors).for_each(
90 |(nodal_coordinate, gradient_vector)| {
91 (0..M).for_each(|i| {
92 (0..M).for_each(|j| {
93 deformation_gradient[i][j] +=
94 nodal_coordinate[i] * gradient_vector[j]
95 })
96 })
97 },
98 );
99 deformation_gradient
100 })
101 .collect()
102 }
103}
104
105pub trait PlanarElasticFiniteElement<C, const G: usize, const N: usize, const P: usize>
106where
107 C: Elastic,
108 Self: PlanarSolidFiniteElement<G, N, P>,
109{
110 fn nodal_forces(
111 &self,
112 constitutive_model: &C,
113 nodal_coordinates: &PlanarElementNodalCoordinates<N>,
114 ) -> Result<PlanarElementNodalForcesSolid<N>, FiniteElementError>;
115 fn nodal_stiffnesses(
116 &self,
117 constitutive_model: &C,
118 nodal_coordinates: &PlanarElementNodalCoordinates<N>,
119 ) -> Result<PlanarElementNodalStiffnessesSolid<N>, FiniteElementError>;
120}
121
122impl<C, const G: usize, const N: usize, const O: usize, const P: usize>
123 PlanarElasticFiniteElement<C, G, N, P> for Element<2, G, N, O>
124where
125 C: Elastic,
126 Self: PlanarSolidFiniteElement<G, N, P>,
127{
128 fn nodal_forces(
129 &self,
130 constitutive_model: &C,
131 nodal_coordinates: &PlanarElementNodalCoordinates<N>,
132 ) -> Result<PlanarElementNodalForcesSolid<N>, FiniteElementError> {
133 match self
134 .deformation_gradients(nodal_coordinates)
135 .iter()
136 .map(|deformation_gradient| {
137 constitutive_model.first_piola_kirchhoff_stress(deformation_gradient)
138 })
139 .collect::<Result<FirstPiolaKirchhoffStressList<G>, _>>()
140 {
141 Ok(first_piola_kirchhoff_stresses) => Ok(first_piola_kirchhoff_stresses
142 .iter()
143 .zip(
144 self.gradient_vectors()
145 .iter()
146 .zip(self.integration_weights()),
147 )
148 .map(
149 |(first_piola_kirchhoff_stress, (gradient_vectors, integration_weight))| {
150 gradient_vectors
151 .iter()
152 .map(|gradient_vector| {
153 (0..M)
154 .map(|i| {
155 (0..M)
156 .map(|j| {
157 first_piola_kirchhoff_stress[i][j]
158 * gradient_vector[j]
159 })
160 .sum::<Scalar>()
161 * integration_weight
162 })
163 .collect()
164 })
165 .collect()
166 },
167 )
168 .sum()),
169 Err(error) => Err(FiniteElementError::Upstream(
170 format!("{error}"),
171 format!("{self:?}"),
172 )),
173 }
174 }
175 fn nodal_stiffnesses(
176 &self,
177 constitutive_model: &C,
178 nodal_coordinates: &PlanarElementNodalCoordinates<N>,
179 ) -> Result<PlanarElementNodalStiffnessesSolid<N>, FiniteElementError> {
180 match self
181 .deformation_gradients(nodal_coordinates)
182 .iter()
183 .map(|deformation_gradient| {
184 constitutive_model.first_piola_kirchhoff_tangent_stiffness(deformation_gradient)
185 })
186 .collect::<Result<FirstPiolaKirchhoffTangentStiffnessList<G>, _>>()
187 {
188 Ok(first_piola_kirchhoff_tangent_stiffnesses) => {
189 Ok(first_piola_kirchhoff_tangent_stiffnesses
190 .iter()
191 .zip(
192 self.gradient_vectors()
193 .iter()
194 .zip(self.integration_weights()),
195 )
196 .map(
197 |(
198 first_piola_kirchhoff_tangent_stiffness,
199 (gradient_vectors, integration_weight),
200 )| {
201 gradient_vectors
202 .iter()
203 .map(|gradient_vector_a| {
204 gradient_vectors
205 .iter()
206 .map(|gradient_vector_b| {
207 (0..M)
208 .map(|i| {
209 (0..M)
210 .map(|k| {
211 (0..M)
212 .map(|j| {
213 (0..M)
214 .map(|l| {
215 first_piola_kirchhoff_tangent_stiffness
216 [i][j][k][l]
217 * gradient_vector_a[j]
218 * gradient_vector_b[l]
219 })
220 .sum::<Scalar>()
221 })
222 .sum::<Scalar>()
223 * integration_weight
224 })
225 .collect()
226 })
227 .collect()
228 })
229 .collect()
230 })
231 .collect()
232 },
233 )
234 .sum())
235 }
236 Err(error) => Err(FiniteElementError::Upstream(
237 format!("{error}"),
238 format!("{self:?}"),
239 )),
240 }
241 }
242}
243
244pub trait PlanarHyperelasticFiniteElement<C, const G: usize, const N: usize, const P: usize>
245where
246 C: Hyperelastic,
247 Self: PlanarElasticFiniteElement<C, G, N, P>,
248{
249 fn helmholtz_free_energy(
250 &self,
251 constitutive_model: &C,
252 nodal_coordinates: &PlanarElementNodalCoordinates<N>,
253 ) -> Result<Scalar, FiniteElementError>;
254}
255
256impl<C, const G: usize, const N: usize, const O: usize, const P: usize>
257 PlanarHyperelasticFiniteElement<C, G, N, P> for Element<2, G, N, O>
258where
259 C: Hyperelastic,
260 Self: PlanarElasticFiniteElement<C, G, N, P>,
261{
262 fn helmholtz_free_energy(
263 &self,
264 constitutive_model: &C,
265 nodal_coordinates: &PlanarElementNodalCoordinates<N>,
266 ) -> Result<Scalar, FiniteElementError> {
267 match self
268 .deformation_gradients(nodal_coordinates)
269 .iter()
270 .zip(self.integration_weights())
271 .map(|(deformation_gradient, integration_weight)| {
272 Ok::<_, ConstitutiveError>(
273 constitutive_model.helmholtz_free_energy_density(deformation_gradient)?
274 * integration_weight,
275 )
276 })
277 .sum()
278 {
279 Ok(helmholtz_free_energy) => Ok(helmholtz_free_energy),
280 Err(error) => Err(FiniteElementError::Upstream(
281 format!("{error}"),
282 format!("{self:?}"),
283 )),
284 }
285 }
286}