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