1#[cfg(test)]
2mod test;
3
4use crate::{
5 fem::block::element::{
6 FRAC_SQRT_3_5, FiniteElement, ParametricCoordinate, ParametricCoordinates,
7 ParametricReference, ShapeFunctions, ShapeFunctionsGradients,
8 quadratic::{M, QuadraticElement, QuadraticFiniteElement},
9 },
10 math::{Scalar, ScalarList, TensorRank1},
11};
12
13const G: usize = 27;
14const N: usize = 13;
15const P: usize = N;
16
17pub type Pyramid = QuadraticElement<G, N>;
18
19impl FiniteElement<G, M, N, P> for Pyramid {
20 fn integration_points() -> ParametricCoordinates<G, M> {
21 const X: [f64; 3] = [0.294997790111502, 0.652996233961648, 0.927005975926850];
22 let u1_2d = [
23 -FRAC_SQRT_3_5,
24 0.0,
25 FRAC_SQRT_3_5,
26 -FRAC_SQRT_3_5,
27 0.0,
28 FRAC_SQRT_3_5,
29 -FRAC_SQRT_3_5,
30 0.0,
31 FRAC_SQRT_3_5,
32 ];
33 let u2_2d = [
34 -FRAC_SQRT_3_5,
35 -FRAC_SQRT_3_5,
36 -FRAC_SQRT_3_5,
37 0.0,
38 0.0,
39 0.0,
40 FRAC_SQRT_3_5,
41 FRAC_SQRT_3_5,
42 FRAC_SQRT_3_5,
43 ];
44 X.into_iter()
45 .flat_map(|x| {
46 u1_2d
47 .into_iter()
48 .zip(u2_2d)
49 .map(move |(u1, u2)| TensorRank1::from([x * u1, x * u2, 1.0 - x]))
50 })
51 .collect()
52 }
53 fn integration_weights(&self) -> &ScalarList<G> {
54 &self.integration_weights
55 }
56 fn parametric_reference() -> ParametricReference<M, N> {
57 [
58 [-1.0, -1.0, 0.0],
59 [1.0, -1.0, 0.0],
60 [1.0, 1.0, 0.0],
61 [-1.0, 1.0, 0.0],
62 [0.0, 0.0, 1.0],
63 [0.0, -1.0, 0.0],
64 [1.0, 0.0, 0.0],
65 [0.0, 1.0, 0.0],
66 [-1.0, 0.0, 0.0],
67 [-0.5, -0.5, 0.5],
68 [0.5, -0.5, 0.5],
69 [0.5, 0.5, 0.5],
70 [-0.5, 0.5, 0.5],
71 ]
72 .into()
73 }
74 fn parametric_weights() -> ScalarList<G> {
75 const B: [f64; 3] = [0.029950703008581, 0.146246269259866, 0.157136361064887];
76 const W1: f64 = 5.0 / 9.0;
77 const W2: f64 = 8.0 / 9.0;
78 let w_2d = [
79 W1 * W1,
80 W2 * W1,
81 W1 * W1,
82 W1 * W2,
83 W2 * W2,
84 W1 * W2,
85 W1 * W1,
86 W2 * W1,
87 W1 * W1,
88 ];
89 B.into_iter()
90 .flat_map(|b| w_2d.into_iter().map(move |w| w * b))
91 .collect()
92 }
93 fn shape_functions(parametric_coordinate: ParametricCoordinate<M>) -> ShapeFunctions<N> {
94 let [xi_1, xi_2, xi_3] = parametric_coordinate.into();
95 let bottom = bottom(xi_3);
96 [
97 0.25 * (-xi_1 - xi_2 - 1.0)
98 * ((1.0 - xi_1) * (1.0 - xi_2) - xi_3 + xi_1 * xi_2 * xi_3 / bottom),
99 0.25 * (xi_1 - xi_2 - 1.0)
100 * ((1.0 + xi_1) * (1.0 - xi_2) - xi_3 - xi_1 * xi_2 * xi_3 / bottom),
101 0.25 * (xi_1 + xi_2 - 1.0)
102 * ((1.0 + xi_1) * (1.0 + xi_2) - xi_3 + xi_1 * xi_2 * xi_3 / bottom),
103 0.25 * (-xi_1 + xi_2 - 1.0)
104 * ((1.0 - xi_1) * (1.0 + xi_2) - xi_3 - xi_1 * xi_2 * xi_3 / bottom),
105 xi_3 * (2.0 * xi_3 - 1.0),
106 0.5 * (1.0 + xi_1 - xi_3) * (1.0 - xi_1 - xi_3) * (1.0 - xi_2 - xi_3) / bottom,
107 0.5 * (1.0 + xi_2 - xi_3) * (1.0 - xi_2 - xi_3) * (1.0 + xi_1 - xi_3) / bottom,
108 0.5 * (1.0 + xi_1 - xi_3) * (1.0 - xi_1 - xi_3) * (1.0 + xi_2 - xi_3) / bottom,
109 0.5 * (1.0 + xi_2 - xi_3) * (1.0 - xi_2 - xi_3) * (1.0 - xi_1 - xi_3) / bottom,
110 xi_3 * (1.0 - xi_1 - xi_3) * (1.0 - xi_2 - xi_3) / bottom,
111 xi_3 * (1.0 + xi_1 - xi_3) * (1.0 - xi_2 - xi_3) / bottom,
112 xi_3 * (1.0 + xi_1 - xi_3) * (1.0 + xi_2 - xi_3) / bottom,
113 xi_3 * (1.0 - xi_1 - xi_3) * (1.0 + xi_2 - xi_3) / bottom,
114 ]
115 .into()
116 }
117 fn shape_functions_gradients(
118 parametric_coordinate: ParametricCoordinate<M>,
119 ) -> ShapeFunctionsGradients<M, N> {
120 let [xi_1, xi_2, xi_3] = parametric_coordinate.into();
121 let bottom = bottom(xi_3);
122 let bottom_squared = bottom * bottom;
123 [
124 [
125 0.25 * ((-xi_1 - xi_2 - 1.0) * (-1.0 + xi_2 + xi_2 * xi_3 / bottom)
126 - ((1.0 - xi_1) * (1.0 - xi_2) - xi_3 + xi_1 * xi_2 * xi_3 / bottom)),
127 0.25 * ((-xi_1 - xi_2 - 1.0) * (-1.0 + xi_1 + xi_1 * xi_3 / bottom)
128 - ((1.0 - xi_1) * (1.0 - xi_2) - xi_3 + xi_1 * xi_2 * xi_3 / bottom)),
129 0.25 * (-xi_1 - xi_2 - 1.0)
130 * (-1.0 + xi_1 * xi_2 / bottom + xi_1 * xi_2 * xi_3 / bottom_squared),
131 ],
132 [
133 0.25 * ((xi_1 - xi_2 - 1.0) * (1.0 - xi_2 - xi_2 * xi_3 / bottom)
134 + ((1.0 + xi_1) * (1.0 - xi_2) - xi_3 - xi_1 * xi_2 * xi_3 / bottom)),
135 0.25 * ((xi_1 - xi_2 - 1.0) * (-1.0 - xi_1 - xi_1 * xi_3 / bottom)
136 - ((1.0 + xi_1) * (1.0 - xi_2) - xi_3 - xi_1 * xi_2 * xi_3 / bottom)),
137 0.25 * (xi_1 - xi_2 - 1.0)
138 * (-1.0 - xi_1 * xi_2 / bottom - xi_1 * xi_2 * xi_3 / bottom_squared),
139 ],
140 [
141 0.25 * ((xi_1 + xi_2 - 1.0) * (1.0 + xi_2 + xi_2 * xi_3 / bottom)
142 + ((1.0 + xi_1) * (1.0 + xi_2) - xi_3 + xi_1 * xi_2 * xi_3 / bottom)),
143 0.25 * ((xi_1 + xi_2 - 1.0) * (1.0 + xi_1 + xi_1 * xi_3 / bottom)
144 + ((1.0 + xi_1) * (1.0 + xi_2) - xi_3 + xi_1 * xi_2 * xi_3 / bottom)),
145 0.25 * (xi_1 + xi_2 - 1.0)
146 * (-1.0 + xi_1 * xi_2 / bottom + xi_1 * xi_2 * xi_3 / bottom_squared),
147 ],
148 [
149 0.25 * ((-xi_1 + xi_2 - 1.0) * (-1.0 - xi_2 - xi_2 * xi_3 / bottom)
150 - ((1.0 - xi_1) * (1.0 + xi_2) - xi_3 - xi_1 * xi_2 * xi_3 / bottom)),
151 0.25 * ((-xi_1 + xi_2 - 1.0) * (1.0 - xi_1 - xi_1 * xi_3 / bottom)
152 + ((1.0 - xi_1) * (1.0 + xi_2) - xi_3 - xi_1 * xi_2 * xi_3 / bottom)),
153 0.25 * (-xi_1 + xi_2 - 1.0)
154 * (-1.0 - xi_1 * xi_2 / bottom - xi_1 * xi_2 * xi_3 / bottom_squared),
155 ],
156 [0.0, 0.0, 4.0 * xi_3 - 1.0],
157 [
158 -xi_1 * (1.0 - xi_2 - xi_3) / bottom,
159 -0.5 * (1.0 - xi_1 - xi_3) * (1.0 + xi_1 - xi_3) / bottom,
160 0.5 * (xi_1 * xi_1 * xi_2 / bottom_squared + xi_2) - 1.0 + xi_3,
161 ],
162 [
163 0.5 * (1.0 - xi_2 - xi_3) * (1.0 + xi_2 - xi_3) / bottom,
164 -xi_2 * (1.0 + xi_1 - xi_3) / bottom,
165 -0.5 * (xi_1 * xi_2 * xi_2 / bottom_squared + xi_1) - 1.0 + xi_3,
166 ],
167 [
168 -xi_1 * (1.0 + xi_2 - xi_3) / bottom,
169 0.5 * (1.0 - xi_1 - xi_3) * (1.0 + xi_1 - xi_3) / bottom,
170 -0.5 * (xi_1 * xi_1 * xi_2 / bottom_squared + xi_2) - 1.0 + xi_3,
171 ],
172 [
173 -0.5 * (1.0 - xi_2 - xi_3) * (1.0 + xi_2 - xi_3) / bottom,
174 -xi_2 * (1.0 - xi_1 - xi_3) / bottom,
175 0.5 * (xi_1 * xi_2 * xi_2 / bottom_squared + xi_1) - 1.0 + xi_3,
176 ],
177 [
178 -(1.0 - xi_2 - xi_3) * xi_3 / bottom,
179 -(1.0 - xi_1 - xi_3) * xi_3 / bottom,
180 xi_1 * xi_2 / bottom_squared + 1.0 - xi_1 - xi_2 - 2.0 * xi_3,
181 ],
182 [
183 (1.0 - xi_2 - xi_3) * xi_3 / bottom,
184 -(1.0 + xi_1 - xi_3) * xi_3 / bottom,
185 -xi_1 * xi_2 / bottom_squared + 1.0 + xi_1 - xi_2 - 2.0 * xi_3,
186 ],
187 [
188 (1.0 + xi_2 - xi_3) * xi_3 / bottom,
189 (1.0 + xi_1 - xi_3) * xi_3 / bottom,
190 xi_1 * xi_2 / bottom_squared + 1.0 + xi_1 + xi_2 - 2.0 * xi_3,
191 ],
192 [
193 -(1.0 + xi_2 - xi_3) * xi_3 / bottom,
194 (1.0 - xi_1 - xi_3) * xi_3 / bottom,
195 -xi_1 * xi_2 / bottom_squared + 1.0 - xi_1 + xi_2 - 2.0 * xi_3,
196 ],
197 ]
198 .into()
199 }
200}
201
202fn bottom(xi_3: Scalar) -> Scalar {
203 const SMALL: Scalar = 4e1 * f64::EPSILON;
204 if (1.0 - xi_3).abs() > SMALL {
205 1.0 - xi_3
206 } else {
207 SMALL
208 }
209}
210
211impl QuadraticFiniteElement<G, N> for Pyramid {}