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},
11 units::Volume,
12};
13
14const G: usize = 18;
15const N: usize = 15;
16const P: usize = N;
17
18pub type Wedge = QuadraticElement<G, N>;
19
20impl FiniteElement<G, M, N, P> for Wedge {
21 fn integration_points() -> ParametricCoordinates<G, M> {
22 const ONE_SIXTH: Scalar = 1.0 / 12.0;
23 const TWO_THIRDS: Scalar = 2.0 / 3.0;
24 [
25 [ONE_SIXTH, ONE_SIXTH, FRAC_SQRT_3_5],
26 [ONE_SIXTH, TWO_THIRDS, FRAC_SQRT_3_5],
27 [TWO_THIRDS, ONE_SIXTH, FRAC_SQRT_3_5],
28 [ONE_SIXTH, ONE_SIXTH, -FRAC_SQRT_3_5],
29 [ONE_SIXTH, TWO_THIRDS, -FRAC_SQRT_3_5],
30 [TWO_THIRDS, ONE_SIXTH, -FRAC_SQRT_3_5],
31 [ONE_SIXTH, ONE_SIXTH, 0.0],
32 [ONE_SIXTH, TWO_THIRDS, 0.0],
33 [TWO_THIRDS, ONE_SIXTH, 0.0],
34 [0.5, 0.5, FRAC_SQRT_3_5],
35 [0.5, 0.0, FRAC_SQRT_3_5],
36 [0.0, 0.5, FRAC_SQRT_3_5],
37 [0.5, 0.5, -FRAC_SQRT_3_5],
38 [0.5, 0.0, -FRAC_SQRT_3_5],
39 [0.0, 0.5, -FRAC_SQRT_3_5],
40 [0.5, 0.5, 0.0],
41 [0.5, 0.0, 0.0],
42 [0.0, 0.5, 0.0],
43 ]
44 .into()
45 }
46 fn integration_weights(&self) -> &IntegrationWeights<G, Volume> {
47 &self.integration_weights
48 }
49 fn parametric_reference() -> ParametricReference<M, N> {
50 [
51 [0.0, 0.0, -1.0],
52 [1.0, 0.0, -1.0],
53 [0.0, 1.0, -1.0],
54 [0.0, 0.0, 1.0],
55 [1.0, 0.0, 1.0],
56 [0.0, 1.0, 1.0],
57 [0.5, 0.0, -1.0],
58 [0.5, 0.5, -1.0],
59 [0.0, 0.5, -1.0],
60 [0.0, 0.0, 0.0],
61 [1.0, 0.0, 0.0],
62 [0.0, 1.0, 0.0],
63 [0.5, 0.0, 1.0],
64 [0.5, 0.5, 1.0],
65 [0.0, 0.5, 1.0],
66 ]
67 .into()
68 }
69 fn parametric_weights() -> ScalarList<G> {
70 const ONE_TWELFTH: Scalar = 1.0 / 12.0;
71 const TWO_FIFTEENTHS: Scalar = 2.0 / 15.0;
72 const ONE_ONE_HUNDRED_EIGHTH: Scalar = 1.0 / 108.0;
73 const TWO_ONE_HUNDRED_THIRTY_FIFTHS: Scalar = 2.0 / 135.0;
74 [
75 ONE_TWELFTH,
76 ONE_TWELFTH,
77 ONE_TWELFTH,
78 ONE_TWELFTH,
79 ONE_TWELFTH,
80 ONE_TWELFTH,
81 TWO_FIFTEENTHS,
82 TWO_FIFTEENTHS,
83 TWO_FIFTEENTHS,
84 ONE_ONE_HUNDRED_EIGHTH,
85 ONE_ONE_HUNDRED_EIGHTH,
86 ONE_ONE_HUNDRED_EIGHTH,
87 ONE_ONE_HUNDRED_EIGHTH,
88 ONE_ONE_HUNDRED_EIGHTH,
89 ONE_ONE_HUNDRED_EIGHTH,
90 TWO_ONE_HUNDRED_THIRTY_FIFTHS,
91 TWO_ONE_HUNDRED_THIRTY_FIFTHS,
92 TWO_ONE_HUNDRED_THIRTY_FIFTHS,
93 ]
94 .into()
95 }
96 fn shape_functions(parametric_coordinate: ParametricCoordinate<M>) -> ShapeFunctions<N> {
97 let [xi_1, xi_2, xi_3] = parametric_coordinate.into();
98 let xi_0 = 1.0 - xi_1 - xi_2;
99 [
100 -0.5 * xi_0 * (1.0 - xi_3) * (2.0 * xi_1 + 2.0 * xi_2 + xi_3),
101 0.5 * xi_1 * (1.0 - xi_3) * (2.0 * xi_1 - 2.0 - xi_3),
102 0.5 * xi_2 * (1.0 - xi_3) * (2.0 * xi_2 - 2.0 - xi_3),
103 -0.5 * xi_0 * (1.0 + xi_3) * (2.0 * xi_1 + 2.0 * xi_2 - xi_3),
104 0.5 * xi_1 * (1.0 + xi_3) * (2.0 * xi_1 - 2.0 + xi_3),
105 0.5 * xi_2 * (1.0 + xi_3) * (2.0 * xi_2 - 2.0 + xi_3),
106 2.0 * xi_0 * xi_1 * (1.0 - xi_3),
107 2.0 * xi_1 * xi_2 * (1.0 - xi_3),
108 2.0 * xi_0 * xi_2 * (1.0 - xi_3),
109 xi_0 * (1.0 - xi_3 * xi_3),
110 xi_1 * (1.0 - xi_3 * xi_3),
111 xi_2 * (1.0 - xi_3 * xi_3),
112 2.0 * xi_0 * xi_1 * (1.0 + xi_3),
113 2.0 * xi_1 * xi_2 * (1.0 + xi_3),
114 2.0 * xi_0 * xi_2 * (1.0 + xi_3),
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 xi_0 = 1.0 - xi_1 - xi_2;
123 [
124 [
125 0.5 * (1.0 - xi_3) * (4.0 * xi_1 + 4.0 * xi_2 + xi_3 - 2.0),
126 0.5 * (1.0 - xi_3) * (4.0 * xi_1 + 4.0 * xi_2 + xi_3 - 2.0),
127 xi_0 * (xi_1 + xi_2 + xi_3 - 0.5),
128 ],
129 [
130 0.5 * (1.0 - xi_3) * (4.0 * xi_1 - 2.0 - xi_3),
131 0.0,
132 xi_1 * (xi_3 - xi_1 + 0.5),
133 ],
134 [
135 0.0,
136 0.5 * (1.0 - xi_3) * (4.0 * xi_2 - xi_3 - 2.0),
137 xi_2 * (xi_3 - xi_2 + 0.5),
138 ],
139 [
140 0.5 * (1.0 + xi_3) * (4.0 * xi_1 + 4.0 * xi_2 - xi_3 - 2.0),
141 0.5 * (1.0 + xi_3) * (4.0 * xi_1 + 4.0 * xi_2 - xi_3 - 2.0),
142 xi_0 * (-xi_1 - xi_2 + xi_3 + 0.5),
143 ],
144 [
145 0.5 * (1.0 + xi_3) * (4.0 * xi_1 - 2.0 + xi_3),
146 0.0,
147 xi_1 * (xi_3 + xi_1 - 0.5),
148 ],
149 [
150 0.0,
151 0.5 * (1.0 + xi_3) * (4.0 * xi_2 + xi_3 - 2.0),
152 xi_2 * (xi_3 + xi_2 - 0.5),
153 ],
154 [
155 2.0 * (1.0 - xi_3) * (1.0 - 2.0 * xi_1 - xi_2),
156 -2.0 * xi_1 * (1.0 - xi_3),
157 -2.0 * xi_0 * xi_1,
158 ],
159 [
160 2.0 * xi_2 * (1.0 - xi_3),
161 2.0 * xi_1 * (1.0 - xi_3),
162 -2.0 * xi_1 * xi_2,
163 ],
164 [
165 -2.0 * xi_2 * (1.0 - xi_3),
166 2.0 * (1.0 - xi_3) * (xi_0 - xi_2),
167 -2.0 * xi_0 * xi_2,
168 ],
169 [xi_3 * xi_3 - 1.0, xi_3 * xi_3 - 1.0, -2.0 * xi_0 * xi_3],
170 [1.0 - xi_3 * xi_3, 0.0, -2.0 * xi_1 * xi_3],
171 [0.0, 1.0 - xi_3 * xi_3, -2.0 * xi_2 * xi_3],
172 [
173 2.0 * (1.0 + xi_3) * (1.0 - 2.0 * xi_1 - xi_2),
174 -2.0 * xi_1 * (1.0 + xi_3),
175 2.0 * xi_0 * xi_1,
176 ],
177 [
178 2.0 * xi_2 * (1.0 + xi_3),
179 2.0 * xi_1 * (1.0 + xi_3),
180 2.0 * xi_1 * xi_2,
181 ],
182 [
183 -2.0 * xi_2 * (1.0 + xi_3),
184 2.0 * (1.0 + xi_3) * (1.0 - xi_1 - 2.0 * xi_2),
185 2.0 * xi_0 * xi_2,
186 ],
187 ]
188 .into()
189 }
190}
191
192impl QuadraticFiniteElement<G, N> for Wedge {}