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
#[cfg(test)]
pub mod test;

pub mod triangle;

use super::*;

pub trait CompositeSurfaceElement<
    'a,
    C,
    const G: usize,
    const M: usize,
    const N: usize,
    const O: usize,
    const P: usize,
    const Q: usize,
> where
    C: Constitutive<'a>,
    Self: CompositeElement<'a, C, G, M, N, O, P, Q>,
{
    fn calculate_bases<const I: usize>(nodal_coordinates: &Coordinates<I, O>) -> Bases<I, P> {
        Self::calculate_standard_gradient_operators()
            .iter()
            .map(|standard_gradient_operator| {
                standard_gradient_operator
                    .iter()
                    .zip(nodal_coordinates.iter())
                    .map(|(standard_gradient_operator_a, nodal_coordinate_a)| {
                        standard_gradient_operator_a
                            .iter()
                            .map(|standard_gradient_operator_a_m| {
                                nodal_coordinate_a * standard_gradient_operator_a_m
                            })
                            .collect()
                    })
                    .sum()
            })
            .collect()
    }
    fn calculate_deformation_gradients_composite_surface_element(
        &self,
        nodal_coordinates: &NodalCoordinates<O>,
    ) -> DeformationGradients<G> {
        let normals = Self::calculate_normals(nodal_coordinates);
        self.get_projected_gradient_vectors()
            .iter()
            .zip(self.get_scaled_reference_normals().iter())
            .map(|(projected_gradient_vectors, scaled_reference_normals)| {
                nodal_coordinates
                    .iter()
                    .zip(projected_gradient_vectors.iter())
                    .map(|(nodal_coordinate, projected_gradient_vector)| {
                        DeformationGradient::dyad(nodal_coordinate, projected_gradient_vector)
                    })
                    .sum::<DeformationGradient>()
                    + scaled_reference_normals
                        .iter()
                        .zip(normals.iter())
                        .map(|(scaled_reference_normal, normal)| {
                            DeformationGradient::dyad(normal, scaled_reference_normal)
                        })
                        .sum::<DeformationGradient>()
            })
            .collect()
    }
    fn calculate_deformation_gradient_rates_composite_surface_element(
        &self,
        nodal_coordinates: &NodalCoordinates<O>,
        nodal_velocities: &NodalVelocities<O>,
    ) -> DeformationGradientRates<G> {
        let normal_rates = Self::calculate_normal_rates(nodal_coordinates, nodal_velocities);
        self.get_projected_gradient_vectors()
            .iter()
            .zip(self.get_scaled_reference_normals().iter())
            .map(|(projected_gradient_vectors, scaled_reference_normals)| {
                nodal_velocities
                    .iter()
                    .zip(projected_gradient_vectors.iter())
                    .map(|(nodal_velocity, projected_gradient_vector)| {
                        DeformationGradientRate::dyad(nodal_velocity, projected_gradient_vector)
                    })
                    .sum::<DeformationGradientRate>()
                    + scaled_reference_normals
                        .iter()
                        .zip(normal_rates.iter())
                        .map(|(scaled_reference_normal, normal_rate)| {
                            DeformationGradientRate::dyad(normal_rate, scaled_reference_normal)
                        })
                        .sum::<DeformationGradientRate>()
            })
            .collect()
    }
    fn calculate_dual_bases<const I: usize>(nodal_coordinates: &Coordinates<I, O>) -> Bases<I, P> {
        Self::calculate_bases(nodal_coordinates)
            .iter()
            .map(|basis_vectors| {
                basis_vectors
                    .iter()
                    .map(|basis_vectors_m| {
                        basis_vectors
                            .iter()
                            .map(|basis_vectors_n| basis_vectors_m * basis_vectors_n)
                            .collect()
                    })
                    .collect::<TensorRank2<M, I, I>>()
                    .inverse()
                    .iter()
                    .map(|metric_tensor_m| {
                        metric_tensor_m
                            .iter()
                            .zip(basis_vectors.iter())
                            .map(|(metric_tensor_mn, basis_vectors_n)| {
                                basis_vectors_n * metric_tensor_mn
                            })
                            .sum()
                    })
                    .collect()
            })
            .collect()
    }
    fn calculate_normals(nodal_coordinates: &NodalCoordinates<O>) -> Normals<P> {
        Self::calculate_bases(nodal_coordinates)
            .iter()
            .map(|basis_vectors| basis_vectors[0].cross(&basis_vectors[1]).normalized())
            .collect()
    }
    fn calculate_normal_gradients(nodal_coordinates: &Coordinates<1, O>) -> NormalGradientss<P, O> {
        let levi_civita_symbol = LEVI_CIVITA;
        let mut normalization: Scalar = 0.0;
        let mut normal_vector = tensor_rank_1_zero();
        Self::calculate_standard_gradient_operators().iter()
        .zip(Self::calculate_bases(nodal_coordinates).iter())
        .map(|(standard_gradient_operator, basis_vectors)|{
            normalization = basis_vectors[0].cross(&basis_vectors[1]).norm();
            normal_vector = basis_vectors[0].cross(&basis_vectors[1])/normalization;
            standard_gradient_operator.iter()
            .map(|standard_gradient_operator_a|
                levi_civita_symbol.iter()
                .map(|levi_civita_symbol_m|
                    IDENTITY.iter()
                    .zip(normal_vector.iter())
                    .map(|(identity_i, normal_vector_i)|
                        levi_civita_symbol_m.iter()
                        .zip(basis_vectors[0].iter()
                        .zip(basis_vectors[1].iter()))
                        .map(|(levi_civita_symbol_mn, (basis_vector_0_n, basis_vector_1_n))|
                            levi_civita_symbol_mn.iter()
                            .zip(identity_i.iter()
                            .zip(normal_vector.iter()))
                            .map(|(levi_civita_symbol_mno, (identity_io, normal_vector_o))|
                                levi_civita_symbol_mno * (identity_io - normal_vector_i * normal_vector_o)
                            ).sum::<Scalar>() * (
                                standard_gradient_operator_a[0] * basis_vector_1_n
                              - standard_gradient_operator_a[1] * basis_vector_0_n
                            )
                        ).sum::<Scalar>() / normalization
                    ).collect()
                ).collect()
            ).collect()
        }).collect()
    }
    fn calculate_normal_rates(
        nodal_coordinates: &NodalCoordinates<O>,
        nodal_velocities: &NodalVelocities<O>,
    ) -> NormalRates<P> {
        let levi_civita_symbol = LEVI_CIVITA;
        let mut normalization: Scalar = 0.0;
        let mut normal_vector = tensor_rank_1_zero();
        Self::calculate_standard_gradient_operators().iter()
        .zip(Self::calculate_bases(nodal_coordinates).iter())
        .map(|(standard_gradient_operator, basis_vectors)|{
            normalization = basis_vectors[0].cross(&basis_vectors[1]).norm();
            normal_vector = basis_vectors[0].cross(&basis_vectors[1])/normalization;
            IDENTITY.iter()
            .zip(normal_vector.iter())
            .map(|(identity_i, normal_vector_i)|
                nodal_velocities.iter()
                .zip(standard_gradient_operator.iter())
                .map(|(nodal_velocity_a, standard_gradient_operator_a)|
                    levi_civita_symbol.iter()
                    .zip(nodal_velocity_a.iter())
                    .map(|(levi_civita_symbol_m, nodal_velocity_a_m)|
                        levi_civita_symbol_m.iter()
                        .zip(basis_vectors[0].iter()
                        .zip(basis_vectors[1].iter()))
                        .map(|(levi_civita_symbol_mn, (basis_vector_0_n, basis_vector_1_n))|
                            levi_civita_symbol_mn.iter()
                            .zip(identity_i.iter()
                            .zip(normal_vector.iter()))
                            .map(|(levi_civita_symbol_mno, (identity_io, normal_vector_o))|
                                levi_civita_symbol_mno * (identity_io - normal_vector_i * normal_vector_o)
                            ).sum::<Scalar>() * (
                                standard_gradient_operator_a[0] * basis_vector_1_n
                              - standard_gradient_operator_a[1] * basis_vector_0_n
                            )
                        ).sum::<Scalar>() * nodal_velocity_a_m
                    ).sum::<Scalar>()
                ).sum::<Scalar>() / normalization
            ).collect()
        }).collect()
    }
    fn calculate_normal_tangents(nodal_coordinates: &Coordinates<1, O>) -> NormalTangentss<P, O> {
        let levi_civita_symbol = LEVI_CIVITA;
        let mut normalization: Scalar = 0.0;
        let mut normal_vector = tensor_rank_1_zero();
        Self::calculate_standard_gradient_operators().iter()
        .zip(Self::calculate_bases(nodal_coordinates).iter()
        .zip(Self::calculate_normal_gradients(nodal_coordinates).iter()))
        .map(|(standard_gradient_operator, (basis_vectors, normal_gradients))|{
            normalization = basis_vectors[0].cross(&basis_vectors[1]).norm();
            normal_vector = basis_vectors[0].cross(&basis_vectors[1])/normalization;
            normal_gradients.iter()
            .zip(standard_gradient_operator.iter())
            .map(|(normal_gradient_a, standard_gradient_operator_a)|
                normal_gradients.iter()
                .zip(standard_gradient_operator.iter())
                .map(|(normal_gradient_b, standard_gradient_operator_b)|
                    normal_gradient_a.iter()
                    .zip(levi_civita_symbol.iter())
                    .map(|(normal_gradient_a_m, levi_civita_symbol_m)|
                        normal_gradient_b.iter()
                        .zip(levi_civita_symbol.iter()
                        .zip(levi_civita_symbol_m.iter()))
                        .map(|(normal_gradient_b_n, (levi_civita_symbol_n, levi_civita_symbol_mn))|
                            normal_gradient_a_m.iter()
                            .zip(normal_gradient_b_n.iter()
                            .zip(normal_vector.iter()
                            .zip(IDENTITY.iter())))
                            .map(|(normal_gradient_a_m_i, (normal_gradient_b_n_i, (normal_vector_i, identity_i)))|
                                (levi_civita_symbol_m.iter()
                                .zip(levi_civita_symbol_n.iter()
                                .zip(basis_vectors[0].iter()
                                .zip(basis_vectors[1].iter())))
                                .map(|(levi_civita_symbol_mr, (levi_civita_symbol_nr, (basis_vector_0_r, basis_vector_1_r)))|
                                    levi_civita_symbol_mr.iter()
                                    .zip(normal_vector.iter()
                                    .zip(normal_gradient_b_n.iter()))
                                    .map(|(levi_civita_symbol_mrs, (normal_vector_s, normal_gradient_b_n_s))|
                                        levi_civita_symbol_mrs * (
                                            normal_gradient_b_n_i * normal_vector_s
                                        + normal_gradient_b_n_s * normal_vector_i
                                        )
                                    ).sum::<Scalar>() * (
                                        standard_gradient_operator_a[1] * basis_vector_0_r
                                    - standard_gradient_operator_a[0] * basis_vector_1_r
                                    ) +
                                    levi_civita_symbol_nr.iter()
                                    .zip(normal_vector.iter())
                                    .map(|(levi_civita_symbol_nrs, normal_vector_s)|
                                        levi_civita_symbol_nrs * normal_vector_s * normal_gradient_a_m_i
                                    ).sum::<Scalar>() * (
                                        standard_gradient_operator_b[1] * basis_vector_0_r
                                    - standard_gradient_operator_b[0] * basis_vector_1_r
                                    )
                                ).sum::<Scalar>() +
                                levi_civita_symbol_mn * (
                                    identity_i - &normal_vector * normal_vector_i
                                ) * (
                                    standard_gradient_operator_a[0] * standard_gradient_operator_b[1]
                                - standard_gradient_operator_a[1] * standard_gradient_operator_b[0]
                                )) / normalization
                            ).collect()
                        ).collect()
                    ).collect()
                ).collect()
            ).collect()
        }).collect()
    }
    fn calculate_objects(
        &self,
        normal_gradients: &NormalGradientss<P, O>,
    ) -> TensorRank3List2D<3, 1, 1, 0, O, G> {
        self.get_scaled_reference_normals()
            .iter()
            .map(|scaled_reference_normals| {
                normal_gradients
                    .iter()
                    .zip(scaled_reference_normals.iter())
                    .map(|(normal_gradient, scaled_reference_normal)| {
                        normal_gradient
                            .iter()
                            .map(|normal_gradient_a| {
                                normal_gradient_a
                                    .iter()
                                    .map(|normal_gradient_a_m| {
                                        TensorRank2::dyad(
                                            normal_gradient_a_m,
                                            scaled_reference_normal,
                                        )
                                    })
                                    .collect::<TensorRank3<3, 1, 1, 0>>()
                            })
                            .collect::<TensorRank3List<3, 1, 1, 0, O>>()
                    })
                    .sum::<TensorRank3List<3, 1, 1, 0, O>>()
            })
            .collect()
    }
    fn calculate_projected_gradient_vectors_composite_surface_element(
        reference_nodal_coordinates: &ReferenceNodalCoordinates<O>,
    ) -> ProjectedGradientVectors<G, N> {
        let reference_dual_bases = Self::calculate_dual_bases(reference_nodal_coordinates);
        let reference_jacobians_subelements =
            Self::calculate_reference_jacobians_subelements(reference_nodal_coordinates);
        let inverse_projection_matrix =
            Self::calculate_inverse_projection_matrix(&reference_jacobians_subelements);
        Self::calculate_shape_functions_at_integration_points()
            .iter()
            .map(|shape_functions_at_integration_point| {
                Self::calculate_standard_gradient_operators_transposed()
                    .iter()
                    .map(|standard_gradient_operators_a| {
                        Self::calculate_shape_function_integrals()
                            .iter()
                            .zip(
                                standard_gradient_operators_a.iter().zip(
                                    reference_dual_bases
                                        .iter()
                                        .zip(reference_jacobians_subelements.iter()),
                                ),
                            )
                            .map(
                                |(
                                    shape_function_integral,
                                    (
                                        standard_gradient_operator,
                                        (
                                            reference_dual_basis_vectors,
                                            reference_jacobian_subelement,
                                        ),
                                    ),
                                )| {
                                    reference_dual_basis_vectors
                                        .iter()
                                        .zip(standard_gradient_operator.iter())
                                        .map(
                                            |(
                                                reference_dual_basis_vector,
                                                standard_gradient_operator_mu,
                                            )| {
                                                reference_dual_basis_vector
                                                    * standard_gradient_operator_mu
                                            },
                                        )
                                        .sum::<Vector<0>>()
                                        * reference_jacobian_subelement
                                        * (shape_functions_at_integration_point
                                            * (&inverse_projection_matrix
                                                * shape_function_integral))
                                },
                            )
                            .sum()
                    })
                    .collect()
            })
            .collect()
    }
    fn calculate_reference_normals(
        reference_nodal_coordinates: &ReferenceNodalCoordinates<O>,
    ) -> ReferenceNormals<P> {
        Self::calculate_dual_bases(reference_nodal_coordinates)
            .iter()
            .map(|reference_dual_basis_vectors| {
                reference_dual_basis_vectors[0]
                    .cross(&reference_dual_basis_vectors[1])
                    .normalized()
            })
            .collect()
    }
    fn calculate_scaled_reference_normals(
        reference_nodal_coordinates: &ReferenceNodalCoordinates<O>,
    ) -> ScaledReferenceNormals<G, P> {
        let reference_jacobians_subelements =
            Self::calculate_reference_jacobians_subelements(reference_nodal_coordinates);
        let inverse_projection_matrix =
            Self::calculate_inverse_projection_matrix(&reference_jacobians_subelements);
        let reference_normals = Self::calculate_reference_normals(reference_nodal_coordinates);
        let shape_function_integrals = Self::calculate_shape_function_integrals();
        Self::calculate_shape_functions_at_integration_points()
            .iter()
            .map(|shape_function| {
                reference_normals
                    .iter()
                    .zip(
                        shape_function_integrals
                            .iter()
                            .zip(reference_jacobians_subelements.iter()),
                    )
                    .map(
                        |(
                            reference_normal,
                            (shape_function_integral, reference_jacobian_subelement),
                        )| {
                            reference_normal
                                * ((shape_function
                                    * (&inverse_projection_matrix * shape_function_integral))
                                    * reference_jacobian_subelement)
                        },
                    )
                    .collect()
            })
            .collect()
    }
    fn get_scaled_reference_normals(&self) -> &ScaledReferenceNormals<G, P>;
}