1#[cfg(test)]
2mod test;
3
4use crate::{
5 constitutive::{
6 ConstitutiveError,
7 solid::{
8 elastic::Elastic, elastic_hyperviscous::ElasticHyperviscous,
9 hyperelastic::Hyperelastic, hyperviscoelastic::Hyperviscoelastic,
10 viscoelastic::Viscoelastic,
11 },
12 },
13 fem::block::element::{
14 ElementNodalCoordinates, ElementNodalReferenceCoordinates, ElementNodalVelocities,
15 FiniteElementError, StandardGradientOperators,
16 solid::{
17 ElementNodalForcesSolid, ElementNodalStiffnessesSolid, SolidFiniteElement,
18 elastic::ElasticFiniteElement, elastic_hyperviscous::ElasticHyperviscousFiniteElement,
19 hyperelastic::HyperelasticFiniteElement,
20 hyperviscoelastic::HyperviscoelasticFiniteElement,
21 viscoelastic::ViscoelasticFiniteElement,
22 },
23 surface::{
24 SurfaceElement, SurfaceFiniteElement, SurfaceFiniteElementMethods,
25 SurfaceFiniteElementMethodsExtra,
26 },
27 },
28 math::{IDENTITY, Scalar, Tensor},
29 mechanics::{
30 DeformationGradient, DeformationGradientList, DeformationGradientRate,
31 DeformationGradientRateList, FirstPiolaKirchhoffRateTangentStiffnesses,
32 FirstPiolaKirchhoffStresses, FirstPiolaKirchhoffTangentStiffnesses,
33 },
34};
35
36#[cfg(test)]
37use crate::fem::block::element::ShapeFunctionsAtIntegrationPoints;
38
39const G: usize = 1;
40const M: usize = 2;
41const N: usize = 3;
42const P: usize = G;
43
44#[cfg(test)]
45const Q: usize = N;
46
47pub type Triangle = SurfaceElement<G, N, P>;
48
49impl SurfaceFiniteElement<G, N, P> for Triangle {
50 fn new(
51 reference_nodal_coordinates: ElementNodalReferenceCoordinates<N>,
52 thickness: Scalar,
53 ) -> Self {
54 let integration_weights = Self::bases(&reference_nodal_coordinates)
55 .iter()
56 .map(|reference_basis| {
57 reference_basis[0].cross(&reference_basis[1]).norm()
58 * Self::integration_weight()
59 * thickness
60 })
61 .collect();
62 let reference_dual_bases = Self::dual_bases(&reference_nodal_coordinates);
63 let gradient_vectors = Self::standard_gradient_operators()
64 .iter()
65 .zip(reference_dual_bases.iter())
66 .map(|(standard_gradient_operator, reference_dual_basis)| {
67 standard_gradient_operator
68 .iter()
69 .map(|standard_gradient_operator_a| {
70 standard_gradient_operator_a
71 .iter()
72 .zip(reference_dual_basis.iter())
73 .map(|(standard_gradient_operator_a_m, reference_dual_basis_m)| {
74 reference_dual_basis_m * standard_gradient_operator_a_m
75 })
76 .sum()
77 })
78 .collect()
79 })
80 .collect();
81 let reference_normals = reference_dual_bases
82 .iter()
83 .map(|reference_dual_basis| {
84 reference_dual_basis[0]
85 .cross(&reference_dual_basis[1])
86 .normalized()
87 })
88 .collect();
89 Self {
90 gradient_vectors,
91 integration_weights,
92 reference_normals,
93 }
94 }
95}
96
97impl Triangle {
98 const fn integration_weight() -> Scalar {
99 1.0 / 2.0
100 }
101 #[cfg(test)]
102 const fn shape_functions_at_integration_points() -> ShapeFunctionsAtIntegrationPoints<G, Q> {
103 ShapeFunctionsAtIntegrationPoints::<G, Q>::const_from([[1.0 / 3.0; Q]])
104 }
105 const fn standard_gradient_operators() -> StandardGradientOperators<M, N, P> {
106 StandardGradientOperators::<M, N, P>::const_from([[[-1.0, -1.0], [1.0, 0.0], [0.0, 1.0]]])
107 }
108}
109
110impl SurfaceFiniteElementMethodsExtra<M, N, P> for Triangle {
111 fn standard_gradient_operators() -> StandardGradientOperators<M, N, P> {
112 Self::standard_gradient_operators()
113 }
114}
115
116impl SolidFiniteElement<G, N> for Triangle {
117 fn deformation_gradients(
118 &self,
119 nodal_coordinates: &ElementNodalCoordinates<N>,
120 ) -> DeformationGradientList<G> {
121 self.gradient_vectors()
122 .iter()
123 .zip(
124 Self::normals(nodal_coordinates)
125 .iter()
126 .zip(self.reference_normals().iter()),
127 )
128 .map(|(gradient_vectors, normal_and_reference_normal)| {
129 nodal_coordinates
130 .iter()
131 .zip(gradient_vectors.iter())
132 .map(|(nodal_coordinate, gradient_vector)| {
133 (nodal_coordinate, gradient_vector).into()
134 })
135 .sum::<DeformationGradient>()
136 + DeformationGradient::from(normal_and_reference_normal)
137 })
138 .collect()
139 }
140 fn deformation_gradient_rates(
141 &self,
142 nodal_coordinates: &ElementNodalCoordinates<N>,
143 nodal_velocities: &ElementNodalVelocities<N>,
144 ) -> DeformationGradientRateList<G> {
145 self.gradient_vectors()
146 .iter()
147 .zip(
148 Self::normal_rates(nodal_coordinates, nodal_velocities)
149 .iter()
150 .zip(self.reference_normals().iter()),
151 )
152 .map(|(gradient_vectors, normal_rate_and_reference_normal)| {
153 nodal_velocities
154 .iter()
155 .zip(gradient_vectors.iter())
156 .map(|(nodal_velocity, gradient_vector)| {
157 (nodal_velocity, gradient_vector).into()
158 })
159 .sum::<DeformationGradientRate>()
160 + DeformationGradientRate::from(normal_rate_and_reference_normal)
161 })
162 .collect()
163 }
164}
165
166impl<C> ElasticFiniteElement<C, G, N> for Triangle
167where
168 C: Elastic,
169{
170 fn nodal_forces(
171 &self,
172 constitutive_model: &C,
173 nodal_coordinates: &ElementNodalCoordinates<N>,
174 ) -> Result<ElementNodalForcesSolid<N>, FiniteElementError> {
175 match self
176 .deformation_gradients(nodal_coordinates)
177 .iter()
178 .map(|deformation_gradient| {
179 constitutive_model.first_piola_kirchhoff_stress(deformation_gradient)
180 })
181 .collect::<Result<FirstPiolaKirchhoffStresses<G>, _>>()
182 {
183 Ok(first_piola_kirchhoff_stresses) => Ok(first_piola_kirchhoff_stresses
184 .iter()
185 .zip(
186 self.gradient_vectors()
187 .iter()
188 .zip(self.integration_weights().iter()),
189 )
190 .map(
191 |(first_piola_kirchhoff_stress, (gradient_vectors, integration_weight))| {
192 gradient_vectors
193 .iter()
194 .map(|gradient_vector| {
195 (first_piola_kirchhoff_stress * gradient_vector)
196 * integration_weight
197 })
198 .collect()
199 },
200 )
201 .sum()),
202 Err(error) => Err(FiniteElementError::Upstream(
203 format!("{error}"),
204 format!("{self:?}"),
205 )),
206 }
207 }
208 fn nodal_stiffnesses(
209 &self,
210 constitutive_model: &C,
211 nodal_coordinates: &ElementNodalCoordinates<N>,
212 ) -> Result<ElementNodalStiffnessesSolid<N>, FiniteElementError> {
213 match self.deformation_gradients(nodal_coordinates).iter()
214 .map(|deformation_gradient| {
215 constitutive_model.first_piola_kirchhoff_tangent_stiffness(deformation_gradient)
216 })
217 .collect::<Result<FirstPiolaKirchhoffTangentStiffnesses<G>, _>>()
218 {
219 Ok(first_piola_kirchhoff_tangent_stiffnesses) => Ok(first_piola_kirchhoff_tangent_stiffnesses
220 .iter()
221 .zip(
222 self.gradient_vectors()
223 .iter()
224 .zip(self.integration_weights().iter()
225 .zip(self.reference_normals().iter()
226 .zip(Self::normal_gradients(nodal_coordinates).iter())
227 )
228 ),
229 )
230 .map(
231 |(
232 first_piola_kirchhoff_tangent_stiffness,
233 (gradient_vectors, (integration_weight, (reference_normal, normal_gradients))),
234 )| {
235 gradient_vectors.iter()
236 .map(|gradient_vector_a|
237 gradient_vectors.iter()
238 .zip(normal_gradients.iter())
239 .map(|(gradient_vector_b, normal_gradient_b)|
240 first_piola_kirchhoff_tangent_stiffness.iter()
241 .map(|first_piola_kirchhoff_tangent_stiffness_m|
242 IDENTITY.iter()
243 .zip(normal_gradient_b.iter())
244 .map(|(identity_n, normal_gradient_b_n)|
245 first_piola_kirchhoff_tangent_stiffness_m.iter()
246 .zip(gradient_vector_a.iter())
247 .map(|(first_piola_kirchhoff_tangent_stiffness_mj, gradient_vector_a_j)|
248 first_piola_kirchhoff_tangent_stiffness_mj.iter()
249 .zip(identity_n.iter()
250 .zip(normal_gradient_b_n.iter()))
251 .map(|(first_piola_kirchhoff_tangent_stiffness_mjk, (identity_nk, normal_gradient_b_n_k))|
252 first_piola_kirchhoff_tangent_stiffness_mjk.iter()
253 .zip(gradient_vector_b.iter()
254 .zip(reference_normal.iter()))
255 .map(|(first_piola_kirchhoff_tangent_stiffness_mjkl, (gradient_vector_b_l, reference_normal_l))|
256 first_piola_kirchhoff_tangent_stiffness_mjkl * gradient_vector_a_j * (
257 identity_nk * gradient_vector_b_l + normal_gradient_b_n_k * reference_normal_l
258 ) * integration_weight
259 ).sum::<Scalar>()
260 ).sum::<Scalar>()
261 ).sum::<Scalar>()
262 ).collect()
263 ).collect()
264 ).collect()
265 ).collect()
266 }
267 )
268 .sum()),
269 Err(error) => Err(FiniteElementError::Upstream(
270 format!("{error}"),
271 format!("{self:?}"),
272 )),
273 }
274 }
275}
276
277impl<C> HyperelasticFiniteElement<C, G, N> for Triangle
278where
279 C: Hyperelastic,
280{
281 fn helmholtz_free_energy(
282 &self,
283 constitutive_model: &C,
284 nodal_coordinates: &ElementNodalCoordinates<N>,
285 ) -> Result<Scalar, FiniteElementError> {
286 match self
287 .deformation_gradients(nodal_coordinates)
288 .iter()
289 .zip(self.integration_weights().iter())
290 .map(|(deformation_gradient, integration_weight)| {
291 Ok::<_, ConstitutiveError>(
292 constitutive_model.helmholtz_free_energy_density(deformation_gradient)?
293 * integration_weight,
294 )
295 })
296 .sum()
297 {
298 Ok(helmholtz_free_energy) => Ok(helmholtz_free_energy),
299 Err(error) => Err(FiniteElementError::Upstream(
300 format!("{error}"),
301 format!("{self:?}"),
302 )),
303 }
304 }
305}
306
307impl<C> ViscoelasticFiniteElement<C, G, N> for Triangle
308where
309 C: Viscoelastic,
310{
311 fn nodal_forces(
312 &self,
313 constitutive_model: &C,
314 nodal_coordinates: &ElementNodalCoordinates<N>,
315 nodal_velocities: &ElementNodalVelocities<N>,
316 ) -> Result<ElementNodalForcesSolid<N>, FiniteElementError> {
317 match self
318 .deformation_gradients(nodal_coordinates)
319 .iter()
320 .zip(
321 self.deformation_gradient_rates(nodal_coordinates, nodal_velocities)
322 .iter(),
323 )
324 .map(|(deformation_gradient, deformation_gradient_rate)| {
325 constitutive_model
326 .first_piola_kirchhoff_stress(deformation_gradient, deformation_gradient_rate)
327 })
328 .collect::<Result<FirstPiolaKirchhoffStresses<G>, _>>()
329 {
330 Ok(first_piola_kirchhoff_stresses) => Ok(first_piola_kirchhoff_stresses
331 .iter()
332 .zip(
333 self.gradient_vectors()
334 .iter()
335 .zip(self.integration_weights().iter()),
336 )
337 .map(
338 |(first_piola_kirchhoff_stress, (gradient_vectors, integration_weight))| {
339 gradient_vectors
340 .iter()
341 .map(|gradient_vector| {
342 (first_piola_kirchhoff_stress * gradient_vector)
343 * integration_weight
344 })
345 .collect()
346 },
347 )
348 .sum()),
349 Err(error) => Err(FiniteElementError::Upstream(
350 format!("{error}"),
351 format!("{self:?}"),
352 )),
353 }
354 }
355 fn nodal_stiffnesses(
356 &self,
357 constitutive_model: &C,
358 nodal_coordinates: &ElementNodalCoordinates<N>,
359 nodal_velocities: &ElementNodalVelocities<N>,
360 ) -> Result<ElementNodalStiffnessesSolid<N>, FiniteElementError> {
361 match self.deformation_gradients(nodal_coordinates).iter().zip(self.deformation_gradient_rates(nodal_coordinates, nodal_velocities).iter())
362 .map(|(deformation_gradient, deformation_gradient_rate)| {
363 constitutive_model.first_piola_kirchhoff_rate_tangent_stiffness(deformation_gradient, deformation_gradient_rate)
364 })
365 .collect::<Result<FirstPiolaKirchhoffRateTangentStiffnesses<G>, _>>()
366 {
367 Ok(first_piola_kirchhoff_rate_tangent_stiffnesses) => Ok(first_piola_kirchhoff_rate_tangent_stiffnesses
368 .iter()
369 .zip(
370 self.gradient_vectors()
371 .iter()
372 .zip(self.integration_weights().iter()
373 .zip(self.reference_normals().iter()
374 .zip(Self::normal_gradients(nodal_coordinates).iter())
375 )
376 ),
377 )
378 .map(
379 |(
380 first_piola_kirchoff_rate_tangent_stiffness_mjkl,
381 (gradient_vectors, (integration_weight, (reference_normal, normal_gradients))),
382 )| {
383 gradient_vectors.iter()
384 .map(|gradient_vector_a|
385 gradient_vectors.iter()
386 .zip(normal_gradients.iter())
387 .map(|(gradient_vector_b, normal_gradient_b)|
388 first_piola_kirchoff_rate_tangent_stiffness_mjkl.iter()
389 .map(|first_piola_kirchhoff_rate_tangent_stiffness_m|
390 IDENTITY.iter()
391 .zip(normal_gradient_b.iter())
392 .map(|(identity_n, normal_gradient_b_n)|
393 first_piola_kirchhoff_rate_tangent_stiffness_m.iter()
394 .zip(gradient_vector_a.iter())
395 .map(|(first_piola_kirchhoff_rate_tangent_stiffness_mj, gradient_vector_a_j)|
396 first_piola_kirchhoff_rate_tangent_stiffness_mj.iter()
397 .zip(identity_n.iter()
398 .zip(normal_gradient_b_n.iter()))
399 .map(|(first_piola_kirchhoff_rate_tangent_stiffness_mjk, (identity_nk, normal_gradient_b_n_k))|
400 first_piola_kirchhoff_rate_tangent_stiffness_mjk.iter()
401 .zip(gradient_vector_b.iter()
402 .zip(reference_normal.iter()))
403 .map(|(first_piola_kirchoff_rate_tangent_stiffness_mjkl, (gradient_vector_b_l, reference_normal_l))|
404 first_piola_kirchoff_rate_tangent_stiffness_mjkl * gradient_vector_a_j * (
405 identity_nk * gradient_vector_b_l + normal_gradient_b_n_k * reference_normal_l
406 ) * integration_weight
407 ).sum::<Scalar>()
408 ).sum::<Scalar>()
409 ).sum::<Scalar>()
410 ).collect()
411 ).collect()
412 ).collect()
413 ).collect()
414 }
415 )
416 .sum()),
417 Err(error) => Err(FiniteElementError::Upstream(
418 format!("{error}"),
419 format!("{self:?}"),
420 )),
421 }
422 }
423}
424
425impl<C> ElasticHyperviscousFiniteElement<C, G, N> for Triangle
426where
427 C: ElasticHyperviscous,
428{
429 fn viscous_dissipation(
430 &self,
431 constitutive_model: &C,
432 nodal_coordinates: &ElementNodalCoordinates<N>,
433 nodal_velocities: &ElementNodalVelocities<N>,
434 ) -> Result<Scalar, FiniteElementError> {
435 match self
436 .deformation_gradients(nodal_coordinates)
437 .iter()
438 .zip(
439 self.deformation_gradient_rates(nodal_coordinates, nodal_velocities)
440 .iter()
441 .zip(self.integration_weights().iter()),
442 )
443 .map(
444 |(deformation_gradient, (deformation_gradient_rate, integration_weight))| {
445 Ok::<_, ConstitutiveError>(
446 constitutive_model
447 .viscous_dissipation(deformation_gradient, deformation_gradient_rate)?
448 * integration_weight,
449 )
450 },
451 )
452 .sum()
453 {
454 Ok(helmholtz_free_energy) => Ok(helmholtz_free_energy),
455 Err(error) => Err(FiniteElementError::Upstream(
456 format!("{error}"),
457 format!("{self:?}"),
458 )),
459 }
460 }
461 fn dissipation_potential(
462 &self,
463 constitutive_model: &C,
464 nodal_coordinates: &ElementNodalCoordinates<N>,
465 nodal_velocities: &ElementNodalVelocities<N>,
466 ) -> Result<Scalar, FiniteElementError> {
467 match self
468 .deformation_gradients(nodal_coordinates)
469 .iter()
470 .zip(
471 self.deformation_gradient_rates(nodal_coordinates, nodal_velocities)
472 .iter()
473 .zip(self.integration_weights().iter()),
474 )
475 .map(
476 |(deformation_gradient, (deformation_gradient_rate, integration_weight))| {
477 Ok::<_, ConstitutiveError>(
478 constitutive_model.dissipation_potential(
479 deformation_gradient,
480 deformation_gradient_rate,
481 )? * integration_weight,
482 )
483 },
484 )
485 .sum()
486 {
487 Ok(helmholtz_free_energy) => Ok(helmholtz_free_energy),
488 Err(error) => Err(FiniteElementError::Upstream(
489 format!("{error}"),
490 format!("{self:?}"),
491 )),
492 }
493 }
494}
495
496impl<C> HyperviscoelasticFiniteElement<C, G, N> for Triangle
497where
498 C: Hyperviscoelastic,
499{
500 fn helmholtz_free_energy(
501 &self,
502 constitutive_model: &C,
503 nodal_coordinates: &ElementNodalCoordinates<N>,
504 ) -> Result<Scalar, FiniteElementError> {
505 match self
506 .deformation_gradients(nodal_coordinates)
507 .iter()
508 .zip(self.integration_weights().iter())
509 .map(|(deformation_gradient, integration_weight)| {
510 Ok::<_, ConstitutiveError>(
511 constitutive_model.helmholtz_free_energy_density(deformation_gradient)?
512 * integration_weight,
513 )
514 })
515 .sum()
516 {
517 Ok(helmholtz_free_energy) => Ok(helmholtz_free_energy),
518 Err(error) => Err(FiniteElementError::Upstream(
519 format!("{error}"),
520 format!("{self:?}"),
521 )),
522 }
523 }
524}