conspire 0.6.0

The Rust interface to conspire.
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

#[cfg(test)]
mod test;

use crate::{
    fem::block::element::{
        ElementNodalEitherCoordinates, FRAC_SQRT_3_5, FiniteElement, ParametricCoordinate,
        ParametricCoordinates, ParametricReference, ShapeFunctions, ShapeFunctionsGradients,
        quadratic::{M, QuadraticElement, QuadraticFiniteElement},
    },
    math::{Scalar, ScalarList, TensorRank1},
};

const G: usize = 27;
const N: usize = 13;
const P: usize = N;

pub type Pyramid = QuadraticElement<G, N>;

impl FiniteElement<G, M, N, P> for Pyramid {
    fn integration_points() -> ParametricCoordinates<G, M> {
        const X: [f64; 3] = [0.294997790111502, 0.652996233961648, 0.927005975926850];
        let u1_2d = [
            -FRAC_SQRT_3_5,
            0.0,
            FRAC_SQRT_3_5,
            -FRAC_SQRT_3_5,
            0.0,
            FRAC_SQRT_3_5,
            -FRAC_SQRT_3_5,
            0.0,
            FRAC_SQRT_3_5,
        ];
        let u2_2d = [
            -FRAC_SQRT_3_5,
            -FRAC_SQRT_3_5,
            -FRAC_SQRT_3_5,
            0.0,
            0.0,
            0.0,
            FRAC_SQRT_3_5,
            FRAC_SQRT_3_5,
            FRAC_SQRT_3_5,
        ];
        X.into_iter()
            .flat_map(|x| {
                u1_2d
                    .into_iter()
                    .zip(u2_2d)
                    .map(move |(u1, u2)| TensorRank1::from([x * u1, x * u2, 1.0 - x]))
            })
            .collect()
    }
    fn integration_weights(&self) -> &ScalarList<G> {
        &self.integration_weights
    }
    fn parametric_reference() -> ParametricReference<M, N> {
        [
            [-1.0, -1.0, 0.0],
            [1.0, -1.0, 0.0],
            [1.0, 1.0, 0.0],
            [-1.0, 1.0, 0.0],
            [0.0, 0.0, 1.0],
            [0.0, -1.0, 0.0],
            [1.0, 0.0, 0.0],
            [0.0, 1.0, 0.0],
            [-1.0, 0.0, 0.0],
            [-0.5, -0.5, 0.5],
            [0.5, -0.5, 0.5],
            [0.5, 0.5, 0.5],
            [-0.5, 0.5, 0.5],
        ]
        .into()
    }
    fn parametric_weights() -> ScalarList<G> {
        const B: [f64; 3] = [0.029950703008581, 0.146246269259866, 0.157136361064887];
        const W1: f64 = 5.0 / 9.0;
        const W2: f64 = 8.0 / 9.0;
        let w_2d = [
            W1 * W1,
            W2 * W1,
            W1 * W1,
            W1 * W2,
            W2 * W2,
            W1 * W2,
            W1 * W1,
            W2 * W1,
            W1 * W1,
        ];
        B.into_iter()
            .flat_map(|b| w_2d.into_iter().map(move |w| w * b))
            .collect()
    }
    fn scaled_jacobians<const I: usize>(
        _nodal_coordinates: ElementNodalEitherCoordinates<I, N>,
    ) -> ScalarList<P> {
        todo!()
    }
    fn shape_functions(parametric_coordinate: ParametricCoordinate<M>) -> ShapeFunctions<N> {
        let [xi_1, xi_2, xi_3] = parametric_coordinate.into();
        let bottom = bottom(xi_3);
        [
            0.25 * (-xi_1 - xi_2 - 1.0)
                * ((1.0 - xi_1) * (1.0 - xi_2) - xi_3 + xi_1 * xi_2 * xi_3 / bottom),
            0.25 * (xi_1 - xi_2 - 1.0)
                * ((1.0 + xi_1) * (1.0 - xi_2) - xi_3 - xi_1 * xi_2 * xi_3 / bottom),
            0.25 * (xi_1 + xi_2 - 1.0)
                * ((1.0 + xi_1) * (1.0 + xi_2) - xi_3 + xi_1 * xi_2 * xi_3 / bottom),
            0.25 * (-xi_1 + xi_2 - 1.0)
                * ((1.0 - xi_1) * (1.0 + xi_2) - xi_3 - xi_1 * xi_2 * xi_3 / bottom),
            xi_3 * (2.0 * xi_3 - 1.0),
            0.5 * (1.0 + xi_1 - xi_3) * (1.0 - xi_1 - xi_3) * (1.0 - xi_2 - xi_3) / bottom,
            0.5 * (1.0 + xi_2 - xi_3) * (1.0 - xi_2 - xi_3) * (1.0 + xi_1 - xi_3) / bottom,
            0.5 * (1.0 + xi_1 - xi_3) * (1.0 - xi_1 - xi_3) * (1.0 + xi_2 - xi_3) / bottom,
            0.5 * (1.0 + xi_2 - xi_3) * (1.0 - xi_2 - xi_3) * (1.0 - xi_1 - xi_3) / bottom,
            xi_3 * (1.0 - xi_1 - xi_3) * (1.0 - xi_2 - xi_3) / bottom,
            xi_3 * (1.0 + xi_1 - xi_3) * (1.0 - xi_2 - xi_3) / bottom,
            xi_3 * (1.0 + xi_1 - xi_3) * (1.0 + xi_2 - xi_3) / bottom,
            xi_3 * (1.0 - xi_1 - xi_3) * (1.0 + xi_2 - xi_3) / bottom,
        ]
        .into()
    }
    fn shape_functions_gradients(
        parametric_coordinate: ParametricCoordinate<M>,
    ) -> ShapeFunctionsGradients<M, N> {
        let [xi_1, xi_2, xi_3] = parametric_coordinate.into();
        let bottom = bottom(xi_3);
        let bottom_squared = bottom * bottom;
        [
            [
                0.25 * ((-xi_1 - xi_2 - 1.0) * (-1.0 + xi_2 + xi_2 * xi_3 / bottom)
                    - ((1.0 - xi_1) * (1.0 - xi_2) - xi_3 + xi_1 * xi_2 * xi_3 / bottom)),
                0.25 * ((-xi_1 - xi_2 - 1.0) * (-1.0 + xi_1 + xi_1 * xi_3 / bottom)
                    - ((1.0 - xi_1) * (1.0 - xi_2) - xi_3 + xi_1 * xi_2 * xi_3 / bottom)),
                0.25 * (-xi_1 - xi_2 - 1.0)
                    * (-1.0 + xi_1 * xi_2 / bottom + xi_1 * xi_2 * xi_3 / bottom_squared),
            ],
            [
                0.25 * ((xi_1 - xi_2 - 1.0) * (1.0 - xi_2 - xi_2 * xi_3 / bottom)
                    + ((1.0 + xi_1) * (1.0 - xi_2) - xi_3 - xi_1 * xi_2 * xi_3 / bottom)),
                0.25 * ((xi_1 - xi_2 - 1.0) * (-1.0 - xi_1 - xi_1 * xi_3 / bottom)
                    - ((1.0 + xi_1) * (1.0 - xi_2) - xi_3 - xi_1 * xi_2 * xi_3 / bottom)),
                0.25 * (xi_1 - xi_2 - 1.0)
                    * (-1.0 - xi_1 * xi_2 / bottom - xi_1 * xi_2 * xi_3 / bottom_squared),
            ],
            [
                0.25 * ((xi_1 + xi_2 - 1.0) * (1.0 + xi_2 + xi_2 * xi_3 / bottom)
                    + ((1.0 + xi_1) * (1.0 + xi_2) - xi_3 + xi_1 * xi_2 * xi_3 / bottom)),
                0.25 * ((xi_1 + xi_2 - 1.0) * (1.0 + xi_1 + xi_1 * xi_3 / bottom)
                    + ((1.0 + xi_1) * (1.0 + xi_2) - xi_3 + xi_1 * xi_2 * xi_3 / bottom)),
                0.25 * (xi_1 + xi_2 - 1.0)
                    * (-1.0 + xi_1 * xi_2 / bottom + xi_1 * xi_2 * xi_3 / bottom_squared),
            ],
            [
                0.25 * ((-xi_1 + xi_2 - 1.0) * (-1.0 - xi_2 - xi_2 * xi_3 / bottom)
                    - ((1.0 - xi_1) * (1.0 + xi_2) - xi_3 - xi_1 * xi_2 * xi_3 / bottom)),
                0.25 * ((-xi_1 + xi_2 - 1.0) * (1.0 - xi_1 - xi_1 * xi_3 / bottom)
                    + ((1.0 - xi_1) * (1.0 + xi_2) - xi_3 - xi_1 * xi_2 * xi_3 / bottom)),
                0.25 * (-xi_1 + xi_2 - 1.0)
                    * (-1.0 - xi_1 * xi_2 / bottom - xi_1 * xi_2 * xi_3 / bottom_squared),
            ],
            [0.0, 0.0, 4.0 * xi_3 - 1.0],
            [
                -xi_1 * (1.0 - xi_2 - xi_3) / bottom,
                -0.5 * (1.0 - xi_1 - xi_3) * (1.0 + xi_1 - xi_3) / bottom,
                0.5 * (xi_1 * xi_1 * xi_2 / bottom_squared + xi_2) - 1.0 + xi_3,
            ],
            [
                0.5 * (1.0 - xi_2 - xi_3) * (1.0 + xi_2 - xi_3) / bottom,
                -xi_2 * (1.0 + xi_1 - xi_3) / bottom,
                -0.5 * (xi_1 * xi_2 * xi_2 / bottom_squared + xi_1) - 1.0 + xi_3,
            ],
            [
                -xi_1 * (1.0 + xi_2 - xi_3) / bottom,
                0.5 * (1.0 - xi_1 - xi_3) * (1.0 + xi_1 - xi_3) / bottom,
                -0.5 * (xi_1 * xi_1 * xi_2 / bottom_squared + xi_2) - 1.0 + xi_3,
            ],
            [
                -0.5 * (1.0 - xi_2 - xi_3) * (1.0 + xi_2 - xi_3) / bottom,
                -xi_2 * (1.0 - xi_1 - xi_3) / bottom,
                0.5 * (xi_1 * xi_2 * xi_2 / bottom_squared + xi_1) - 1.0 + xi_3,
            ],
            [
                -(1.0 - xi_2 - xi_3) * xi_3 / bottom,
                -(1.0 - xi_1 - xi_3) * xi_3 / bottom,
                xi_1 * xi_2 / bottom_squared + 1.0 - xi_1 - xi_2 - 2.0 * xi_3,
            ],
            [
                (1.0 - xi_2 - xi_3) * xi_3 / bottom,
                -(1.0 + xi_1 - xi_3) * xi_3 / bottom,
                -xi_1 * xi_2 / bottom_squared + 1.0 + xi_1 - xi_2 - 2.0 * xi_3,
            ],
            [
                (1.0 + xi_2 - xi_3) * xi_3 / bottom,
                (1.0 + xi_1 - xi_3) * xi_3 / bottom,
                xi_1 * xi_2 / bottom_squared + 1.0 + xi_1 + xi_2 - 2.0 * xi_3,
            ],
            [
                -(1.0 + xi_2 - xi_3) * xi_3 / bottom,
                (1.0 - xi_1 - xi_3) * xi_3 / bottom,
                -xi_1 * xi_2 / bottom_squared + 1.0 - xi_1 + xi_2 - 2.0 * xi_3,
            ],
        ]
        .into()
    }
}

fn bottom(xi_3: Scalar) -> Scalar {
    const SMALL: Scalar = 4e1 * f64::EPSILON;
    if (1.0 - xi_3).abs() > SMALL {
        1.0 - xi_3
    } else {
        SMALL
    }
}

impl QuadraticFiniteElement<G, N> for Pyramid {}