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