flavio 0.5.0

flavio welcomes you
Documentation
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
#[cfg(test)]
mod test;

pub mod localization;
pub mod surface;
pub mod tetrahedron;

use super::*;

pub trait CompositeElement<
    'a,
    C,
    const G: usize,
    const M: usize,
    const N: usize,
    const O: usize,
    const P: usize,
    const Q: usize,
> where
    C: Constitutive<'a>,
{
    fn calculate_deformation_gradients(
        &self,
        nodal_coordinates: &NodalCoordinates<N>,
    ) -> DeformationGradients<G> {
        self.get_projected_gradient_vectors()
            .iter()
            .map(|projected_gradient_vectors| {
                nodal_coordinates
                    .iter()
                    .zip(projected_gradient_vectors.iter())
                    .map(|(nodal_coordinate, projected_gradient_vector)| {
                        DeformationGradient::dyad(nodal_coordinate, projected_gradient_vector)
                    })
                    .sum()
            })
            .collect()
    }
    fn calculate_deformation_gradient_rates(
        &self,
        _: &NodalCoordinates<N>,
        nodal_velocities: &NodalVelocities<N>,
    ) -> DeformationGradientRates<G> {
        self.get_projected_gradient_vectors()
            .iter()
            .map(|projected_gradient_vectors| {
                nodal_velocities
                    .iter()
                    .zip(projected_gradient_vectors.iter())
                    .map(|(nodal_velocity, projected_gradient_vector)| {
                        DeformationGradientRate::dyad(nodal_velocity, projected_gradient_vector)
                    })
                    .sum()
            })
            .collect()
    }
    fn calculate_inverse_normalized_projection_matrix() -> NormalizedProjectionMatrix<Q>;
    fn calculate_inverse_projection_matrix(
        reference_jacobians_subelements: &Scalars<P>,
    ) -> NormalizedProjectionMatrix<Q> {
        Self::calculate_shape_function_integrals_products()
            .iter()
            .zip(reference_jacobians_subelements.iter())
            .map(
                |(shape_function_integrals_products, reference_jacobian_subelement)| {
                    shape_function_integrals_products * reference_jacobian_subelement
                },
            )
            .sum::<ProjectionMatrix<Q>>()
            .inverse()
    }
    fn calculate_projected_gradient_vectors(
        reference_nodal_coordinates: &ReferenceNodalCoordinates<O>,
    ) -> ProjectedGradientVectors<G, N>;
    fn calculate_reference_jacobians(
        reference_nodal_coordinates: &ReferenceNodalCoordinates<O>,
    ) -> Scalars<G> {
        let vector = Self::calculate_inverse_normalized_projection_matrix()
            * Self::calculate_shape_function_integrals()
                .iter()
                .zip(
                    Self::calculate_reference_jacobians_subelements(reference_nodal_coordinates)
                        .iter(),
                )
                .map(|(shape_function_integral, reference_jacobian_subelement)| {
                    shape_function_integral * reference_jacobian_subelement
                })
                .sum::<TensorRank1<Q, 9>>();
        Self::calculate_shape_functions_at_integration_points()
            .iter()
            .map(|shape_functions_at_integration_point| {
                shape_functions_at_integration_point * &vector
            })
            .collect()
    }
    fn calculate_reference_jacobians_subelements(
        reference_nodal_coordinates: &ReferenceNodalCoordinates<O>,
    ) -> Scalars<P>;
    fn calculate_shape_function_integrals() -> ShapeFunctionIntegrals<P, Q>;
    fn calculate_shape_function_integrals_products() -> ShapeFunctionIntegralsProducts<P, Q>;
    fn calculate_shape_functions_at_integration_points() -> ShapeFunctionsAtIntegrationPoints<G, Q>;
    fn calculate_standard_gradient_operators() -> StandardGradientOperators<M, O, P>;
    fn calculate_standard_gradient_operators_transposed(
    ) -> StandardGradientOperatorsTransposed<M, O, P>;
    fn get_constitutive_models(&self) -> &[C; G];
    fn get_integration_weights(&self) -> &Scalars<G>;
    fn get_projected_gradient_vectors(&self) -> &ProjectedGradientVectors<G, N>;
}

pub trait ElasticCompositeElement<
    'a,
    C,
    const G: usize,
    const M: usize,
    const N: usize,
    const O: usize,
    const P: usize,
    const Q: usize,
> where
    C: Elastic<'a>,
    Self: CompositeElement<'a, C, G, M, N, O, P, Q>,
{
    fn calculate_nodal_forces_composite_element(
        &self,
        nodal_coordinates: &NodalCoordinates<N>,
    ) -> Result<NodalForces<N>, ConstitutiveError> {
        Ok(self
            .get_constitutive_models()
            .iter()
            .zip(
                self.calculate_deformation_gradients(nodal_coordinates)
                    .iter(),
            )
            .map(|(constitutive_model, deformation_gradient)| {
                constitutive_model.calculate_first_piola_kirchoff_stress(deformation_gradient)
            })
            .collect::<Result<FirstPiolaKirchoffStresses<G>, _>>()?
            .iter()
            .zip(
                self.get_projected_gradient_vectors()
                    .iter()
                    .zip(self.get_integration_weights().iter()),
            )
            .map(
                |(
                    first_piola_kirchoff_stress,
                    (projected_gradient_vectors, scaled_composite_jacobian),
                )| {
                    projected_gradient_vectors
                        .iter()
                        .map(|projected_gradient_vector| {
                            (first_piola_kirchoff_stress * projected_gradient_vector)
                                * scaled_composite_jacobian
                        })
                        .collect()
                },
            )
            .sum())
    }
    fn calculate_nodal_stiffnesses_composite_element(
        &self,
        nodal_coordinates: &NodalCoordinates<N>,
    ) -> Result<NodalStiffnesses<N>, ConstitutiveError> {
        Ok(self
            .get_constitutive_models()
            .iter()
            .zip(
                self.calculate_deformation_gradients(nodal_coordinates)
                    .iter(),
            )
            .map(|(constitutive_model, deformation_gradient)| {
                constitutive_model
                    .calculate_first_piola_kirchoff_tangent_stiffness(deformation_gradient)
            })
            .collect::<Result<FirstPiolaKirchoffTangentStiffnesses<G>, _>>()?
            .iter()
            .zip(
                self.get_projected_gradient_vectors()
                    .iter()
                    .zip(self.get_integration_weights().iter()),
            )
            .map(
                |(
                    first_piola_kirchoff_tangent_stiffness,
                    (projected_gradient_vectors, scaled_composite_jacobian),
                )| {
                    projected_gradient_vectors
                        .iter()
                        .map(|projected_gradient_vector_a| {
                            projected_gradient_vectors
                                .iter()
                                .map(|projected_gradient_vector_b| {
                                    first_piola_kirchoff_tangent_stiffness
                                        .contract_second_fourth_indices_with_first_indices_of(
                                            projected_gradient_vector_a,
                                            projected_gradient_vector_b,
                                        )
                                        * scaled_composite_jacobian
                                })
                                .collect()
                        })
                        .collect()
                },
            )
            .sum())
    }
}

pub trait HyperelasticCompositeElement<
    'a,
    C,
    const G: usize,
    const M: usize,
    const N: usize,
    const O: usize,
    const P: usize,
    const Q: usize,
> where
    C: Hyperelastic<'a>,
    Self: ElasticCompositeElement<'a, C, G, M, N, O, P, Q>,
{
    fn calculate_helmholtz_free_energy_composite_element(
        &self,
        nodal_coordinates: &NodalCoordinates<N>,
    ) -> Result<Scalar, ConstitutiveError> {
        self.get_constitutive_models()
            .iter()
            .zip(
                self.calculate_deformation_gradients(nodal_coordinates)
                    .iter()
                    .zip(self.get_integration_weights().iter()),
            )
            .map(
                |(constitutive_model, (deformation_gradient, scaled_composite_jacobian))| {
                    Ok(constitutive_model
                        .calculate_helmholtz_free_energy_density(deformation_gradient)?
                        * scaled_composite_jacobian)
                },
            )
            .sum()
    }
}

pub trait ViscoelasticCompositeElement<
    'a,
    C,
    const G: usize,
    const M: usize,
    const N: usize,
    const O: usize,
    const P: usize,
    const Q: usize,
> where
    C: Viscoelastic<'a>,
    Self: CompositeElement<'a, C, G, M, N, O, P, Q>,
{
    fn calculate_nodal_forces_composite_element(
        &self,
        nodal_coordinates: &NodalCoordinates<N>,
        nodal_velocities: &NodalVelocities<N>,
    ) -> Result<NodalForces<N>, ConstitutiveError> {
        Ok(self
            .get_constitutive_models()
            .iter()
            .zip(
                self.calculate_deformation_gradients(nodal_coordinates)
                    .iter()
                    .zip(
                        self.calculate_deformation_gradient_rates(
                            nodal_coordinates,
                            nodal_velocities,
                        )
                        .iter(),
                    ),
            )
            .map(
                |(constitutive_model, (deformation_gradient, deformation_gradient_rate))| {
                    constitutive_model.calculate_first_piola_kirchoff_stress(
                        deformation_gradient,
                        deformation_gradient_rate,
                    )
                },
            )
            .collect::<Result<FirstPiolaKirchoffStresses<G>, _>>()?
            .iter()
            .zip(
                self.get_projected_gradient_vectors()
                    .iter()
                    .zip(self.get_integration_weights().iter()),
            )
            .map(
                |(
                    first_piola_kirchoff_stress,
                    (projected_gradient_vectors, scaled_composite_jacobian),
                )| {
                    projected_gradient_vectors
                        .iter()
                        .map(|projected_gradient_vector| {
                            (first_piola_kirchoff_stress * projected_gradient_vector)
                                * scaled_composite_jacobian
                        })
                        .collect()
                },
            )
            .sum())
    }
    fn calculate_nodal_stiffnesses_composite_element(
        &self,
        nodal_coordinates: &NodalCoordinates<N>,
        nodal_velocities: &NodalVelocities<N>,
    ) -> Result<NodalStiffnesses<N>, ConstitutiveError> {
        Ok(self
            .get_constitutive_models()
            .iter()
            .zip(
                self.calculate_deformation_gradients(nodal_coordinates)
                    .iter()
                    .zip(
                        self.calculate_deformation_gradient_rates(
                            nodal_coordinates,
                            nodal_velocities,
                        )
                        .iter(),
                    ),
            )
            .map(
                |(constitutive_model, (deformation_gradient, deformation_gradient_rate))| {
                    constitutive_model.calculate_first_piola_kirchoff_rate_tangent_stiffness(
                        deformation_gradient,
                        deformation_gradient_rate,
                    )
                },
            )
            .collect::<Result<FirstPiolaKirchoffRateTangentStiffnesses<G>, _>>()?
            .iter()
            .zip(
                self.get_projected_gradient_vectors().iter().zip(
                    self.get_projected_gradient_vectors()
                        .iter()
                        .zip(self.get_integration_weights().iter()),
                ),
            )
            .map(
                |(
                    first_piola_kirchoff_rate_tangent_stiffness,
                    (
                        projected_gradient_vectors_a,
                        (projected_gradient_vectors_b, scaled_composite_jacobian),
                    ),
                )| {
                    projected_gradient_vectors_a
                        .iter()
                        .map(|projected_gradient_vector_a| {
                            projected_gradient_vectors_b
                                .iter()
                                .map(|projected_gradient_vector_b| {
                                    first_piola_kirchoff_rate_tangent_stiffness
                                        .contract_second_fourth_indices_with_first_indices_of(
                                            projected_gradient_vector_a,
                                            projected_gradient_vector_b,
                                        )
                                        * scaled_composite_jacobian
                                })
                                .collect()
                        })
                        .collect()
                },
            )
            .sum())
    }
}

pub trait ElasticHyperviscousCompositeElement<
    'a,
    C,
    const G: usize,
    const M: usize,
    const N: usize,
    const O: usize,
    const P: usize,
    const Q: usize,
> where
    C: ElasticHyperviscous<'a>,
    Self: ViscoelasticCompositeElement<'a, C, G, M, N, O, P, Q>,
{
    fn calculate_viscous_dissipation_composite_element(
        &self,
        nodal_coordinates: &NodalCoordinates<N>,
        nodal_velocities: &NodalVelocities<N>,
    ) -> Result<Scalar, ConstitutiveError> {
        self.get_constitutive_models()
            .iter()
            .zip(
                self.calculate_deformation_gradients(nodal_coordinates)
                    .iter()
                    .zip(
                        self.calculate_deformation_gradient_rates(
                            nodal_coordinates,
                            nodal_velocities,
                        )
                        .iter()
                        .zip(self.get_integration_weights().iter()),
                    ),
            )
            .map(
                |(
                    constitutive_model,
                    (deformation_gradient, (deformation_gradient_rate, scaled_composite_jacobian)),
                )| {
                    Ok(constitutive_model.calculate_viscous_dissipation(
                        deformation_gradient,
                        deformation_gradient_rate,
                    )? * scaled_composite_jacobian)
                },
            )
            .sum()
    }
    fn calculate_dissipation_potential_composite_element(
        &self,
        nodal_coordinates: &NodalCoordinates<N>,
        nodal_velocities: &NodalVelocities<N>,
    ) -> Result<Scalar, ConstitutiveError> {
        self.get_constitutive_models()
            .iter()
            .zip(
                self.calculate_deformation_gradients(nodal_coordinates)
                    .iter()
                    .zip(
                        self.calculate_deformation_gradient_rates(
                            nodal_coordinates,
                            nodal_velocities,
                        )
                        .iter()
                        .zip(self.get_integration_weights().iter()),
                    ),
            )
            .map(
                |(
                    constitutive_model,
                    (deformation_gradient, (deformation_gradient_rate, scaled_composite_jacobian)),
                )| {
                    Ok(constitutive_model.calculate_dissipation_potential(
                        deformation_gradient,
                        deformation_gradient_rate,
                    )? * scaled_composite_jacobian)
                },
            )
            .sum()
    }
}

pub trait HyperviscoelasticCompositeElement<
    'a,
    C,
    const G: usize,
    const M: usize,
    const N: usize,
    const O: usize,
    const P: usize,
    const Q: usize,
> where
    C: Hyperviscoelastic<'a>,
    Self: ElasticHyperviscousCompositeElement<'a, C, G, M, N, O, P, Q>,
{
    fn calculate_helmholtz_free_energy_composite_element(
        &self,
        nodal_coordinates: &NodalCoordinates<N>,
    ) -> Result<Scalar, ConstitutiveError> {
        self.get_constitutive_models()
            .iter()
            .zip(
                self.calculate_deformation_gradients(nodal_coordinates)
                    .iter()
                    .zip(self.get_integration_weights().iter()),
            )
            .map(
                |(constitutive_model, (deformation_gradient, scaled_composite_jacobian))| {
                    Ok(constitutive_model
                        .calculate_helmholtz_free_energy_density(deformation_gradient)?
                        * scaled_composite_jacobian)
                },
            )
            .sum()
    }
}