49#ifndef __INTREPID2_HGRAD_HEX_C2_FEM_DEF_HPP__
50#define __INTREPID2_HGRAD_HEX_C2_FEM_DEF_HPP__
57 template<
bool serendipity>
58 template<EOperator opType>
59 template<
typename OutputViewType,
60 typename inputViewType>
61 KOKKOS_INLINE_FUNCTION
65 const inputViewType input ) {
67 case OPERATOR_VALUE : {
68 const auto x = input(0);
69 const auto y = input(1);
70 const auto z = input(2);
73 output.access( 0) = 0.125*(-1. + x)*x*(-1. + y)*y*(-1. + z)*z;
74 output.access( 1) = 0.125*x*(1.+ x)*(-1. + y)*y*(-1. + z)*z;
75 output.access( 2) = 0.125*x*(1.+ x)*y*(1.+ y)*(-1. + z)*z;
76 output.access( 3) = 0.125*(-1. + x)*x*y*(1.+ y)*(-1. + z)*z;
77 output.access( 4) = 0.125*(-1. + x)*x*(-1. + y)*y*z*(1.+ z);
78 output.access( 5) = 0.125*x*(1.+ x)*(-1. + y)*y*z*(1.+ z);
79 output.access( 6) = 0.125*x*(1.+ x)*y*(1.+ y)*z*(1.+ z);
80 output.access( 7) = 0.125*(-1. + x)*x*y*(1.+ y)*z*(1.+ z);
81 output.access( 8) = 0.25*(1. - x)*(1. + x)*(-1. + y)*y*(-1. + z)*z;
82 output.access( 9) = 0.25*x*(1.+ x)*(1. - y)*(1. + y)*(-1. + z)*z;
83 output.access(10) = 0.25*(1. - x)*(1. + x)*y*(1.+ y)*(-1. + z)*z;
84 output.access(11) = 0.25*(-1. + x)*x*(1. - y)*(1. + y)*(-1. + z)*z;
85 output.access(12) = 0.25*(-1. + x)*x*(-1. + y)*y*(1. - z)*(1. + z);
86 output.access(13) = 0.25*x*(1.+ x)*(-1. + y)*y*(1. - z)*(1. + z);
87 output.access(14) = 0.25*x*(1.+ x)*y*(1.+ y)*(1. - z)*(1. + z);
88 output.access(15) = 0.25*(-1. + x)*x*y*(1.+ y)*(1. - z)*(1. + z);
89 output.access(16) = 0.25*(1. - x)*(1. + x)*(-1. + y)*y*z*(1.+ z);
90 output.access(17) = 0.25*x*(1.+ x)*(1. - y)*(1. + y)*z*(1.+ z);
91 output.access(18) = 0.25*(1. - x)*(1. + x)*y*(1.+ y)*z*(1.+ z);
92 output.access(19) = 0.25*(-1. + x)*x*(1. - y)*(1. + y)*z*(1.+ z);
93 if constexpr (!serendipity) {
94 output.access(20) = (1. - x)*(1. + x)*(1. - y)*(1. + y)*(1. - z)*(1. + z);
95 output.access(21) = 0.5*(1. - x)*(1. + x)*(1. - y)*(1. + y)*(-1. + z)*z;
96 output.access(22) = 0.5*(1. - x)*(1. + x)*(1. - y)*(1. + y)*z*(1.+ z);
97 output.access(23) = 0.5*(-1. + x)*x*(1. - y)*(1. + y)*(1. - z)*(1. + z);
98 output.access(24) = 0.5*x*(1.+ x)*(1. - y)*(1. + y)*(1. - z)*(1. + z);
99 output.access(25) = 0.5*(1. - x)*(1. + x)*(-1. + y)*y*(1. - z)*(1. + z);
100 output.access(26) = 0.5*(1. - x)*(1. + x)*y*(1.+ y)*(1. - z)*(1. + z);
106 const auto x = input(0);
107 const auto y = input(1);
108 const auto z = input(2);
111 output.access(0, 0) = (-0.125 + 0.25*x)*(-1. + y)*y*(-1. + z)*z;
112 output.access(0, 1) = (-1. + x)*x*(-0.125 + 0.25*y)*(-1. + z)*z;
113 output.access(0, 2) = (-1. + x)*x*(-1. + y)*y*(-0.125 + 0.25*z);
115 output.access(1, 0) = (0.125 + 0.25*x)*(-1. + y)*y*(-1. + z)*z;
116 output.access(1, 1) = x*(1. + x)*(-0.125 + 0.25*y)*(-1. + z)*z;
117 output.access(1, 2) = x*(1. + x)*(-1. + y)*y*(-0.125 + 0.25*z);
119 output.access(2, 0) = (0.125 + 0.25*x)*y*(1. + y)*(-1. + z)*z;
120 output.access(2, 1) = x*(1. + x)*(0.125 + 0.25*y)*(-1. + z)*z;
121 output.access(2, 2) = x*(1. + x)*y*(1. + y)*(-0.125 + 0.25*z);
123 output.access(3, 0) = (-0.125 + 0.25*x)*y*(1. + y)*(-1. + z)*z;
124 output.access(3, 1) = (-1. + x)*x*(0.125 + 0.25*y)*(-1. + z)*z;
125 output.access(3, 2) = (-1. + x)*x*y*(1. + y)*(-0.125 + 0.25*z);
127 output.access(4, 0) = (-0.125 + 0.25*x)*(-1. + y)*y*z*(1. + z);
128 output.access(4, 1) = (-1. + x)*x*(-0.125 + 0.25*y)*z*(1. + z);
129 output.access(4, 2) = (-1. + x)*x*(-1. + y)*y*(0.125 + 0.25*z);
131 output.access(5, 0) = (0.125 + 0.25*x)*(-1. + y)*y*z*(1. + z);
132 output.access(5, 1) = x*(1. + x)*(-0.125 + 0.25*y)*z*(1. + z);
133 output.access(5, 2) = x*(1. + x)*(-1. + y)*y*(0.125 + 0.25*z);
135 output.access(6, 0) = (0.125 + 0.25*x)*y*(1. + y)*z*(1. + z);
136 output.access(6, 1) = x*(1. + x)*(0.125 + 0.25*y)*z*(1. + z);
137 output.access(6, 2) = x*(1. + x)*y*(1. + y)*(0.125 + 0.25*z);
139 output.access(7, 0) = (-0.125 + 0.25*x)*y*(1. + y)*z*(1. + z);
140 output.access(7, 1) = (-1. + x)*x*(0.125 + 0.25*y)*z*(1. + z);
141 output.access(7, 2) = (-1. + x)*x*y*(1. + y)*(0.125 + 0.25*z);
143 output.access(8, 0) = -0.5*x*(-1. + y)*y*(-1. + z)*z;
144 output.access(8, 1) = (1. - x)*(1. + x)*(-0.25 + 0.5*y)*(-1. + z)*z;
145 output.access(8, 2) = (1. - x)*(1. + x)*(-1. + y)*y*(-0.25 + 0.5*z);
147 output.access(9, 0) = (0.25 + 0.5*x)*(1. - y)*(1. + y)*(-1. + z)*z;
148 output.access(9, 1) = x*(1. + x)*(-0.5*y)*(-1. + z)*z;
149 output.access(9, 2) = x*(1. + x)*(1. - y)*(1. + y)*(-0.25 + 0.5*z);
151 output.access(10, 0) = -0.5*x*y*(1. + y)*(-1. + z)*z;
152 output.access(10, 1) = (1. - x)*(1. + x)*(0.25 + 0.5*y)*(-1. + z)*z;
153 output.access(10, 2) = (1. - x)*(1. + x)*y*(1. + y)*(-0.25 + 0.5*z);
155 output.access(11, 0) = (-0.25 + 0.5*x)*(1. - y)*(1. + y)*(-1. + z)*z;
156 output.access(11, 1) = (-1. + x)*x*(-0.5*y)*(-1. + z)*z;
157 output.access(11, 2) = (-1. + x)*x*(1. - y)*(1. + y)*(-0.25 + 0.5*z);
159 output.access(12, 0) = (-0.25 + 0.5*x)*(-1. + y)*y*(1. - z)*(1. + z);
160 output.access(12, 1) = (-1. + x)*x*(-0.25 + 0.5*y)*(1. - z)*(1. + z);
161 output.access(12, 2) = (-1. + x)*x*(-1. + y)*y*(-0.5*z);
163 output.access(13, 0) = (0.25 + 0.5*x)*(-1. + y)*y*(1. - z)*(1. + z);
164 output.access(13, 1) = x*(1. + x)*(-0.25 + 0.5*y)*(1. - z)*(1. + z);
165 output.access(13, 2) = x*(1. + x)*(-1. + y)*y*(-0.5*z);
167 output.access(14, 0) = (0.25 + 0.5*x)*y*(1. + y)*(1. - z)*(1. + z);
168 output.access(14, 1) = x*(1. + x)*(0.25 + 0.5*y)*(1. - z)*(1. + z);
169 output.access(14, 2) = x*(1. + x)*y*(1. + y)*(-0.5*z);
171 output.access(15, 0) = (-0.25 + 0.5*x)*y*(1. + y)*(1. - z)*(1. + z);
172 output.access(15, 1) = (-1. + x)*x*(0.25 + 0.5*y)*(1. - z)*(1. + z);
173 output.access(15, 2) = (-1. + x)*x*y*(1. + y)*(-0.5*z);
175 output.access(16, 0) = -0.5*x*(-1. + y)*y*z*(1. + z);
176 output.access(16, 1) = (1. - x)*(1. + x)*(-0.25 + 0.5*y)*z*(1. + z);
177 output.access(16, 2) = (1. - x)*(1. + x)*(-1. + y)*y*(0.25 + 0.5*z);
179 output.access(17, 0) = (0.25 + 0.5*x)*(1. - y)*(1. + y)*z*(1. + z);
180 output.access(17, 1) = x*(1. + x)*(-0.5*y)*z*(1. + z);
181 output.access(17, 2) = x*(1. + x)*(1. - y)*(1. + y)*(0.25 + 0.5*z);
183 output.access(18, 0) = -0.5*x*y*(1. + y)*z*(1. + z);
184 output.access(18, 1) = (1. - x)*(1. + x)*(0.25 + 0.5*y)*z*(1. + z);
185 output.access(18, 2) = (1. - x)*(1. + x)*y*(1. + y)*(0.25 + 0.5*z);
187 output.access(19, 0) = (-0.25 + 0.5*x)*(1. - y)*(1. + y)*z*(1. + z);
188 output.access(19, 1) = (-1. + x)*x*(-0.5*y)*z*(1. + z);
189 output.access(19, 2) = (-1. + x)*x*(1. - y)*(1. + y)*(0.25 + 0.5*z);
191 if constexpr (!serendipity) {
192 output.access(20, 0) = -2.*x*(1. - y)*(1. + y)*(1. - z)*(1. + z);
193 output.access(20, 1) = (1. - x)*(1. + x)*(-2.*y)*(1. - z)*(1. + z);
194 output.access(20, 2) = (1. - x)*(1. + x)*(1. - y)*(1. + y)*(-2.*z);
196 output.access(21, 0) = -x*(1. - y)*(1. + y)*(-1. + z)*z;
197 output.access(21, 1) = (1. - x)*(1. + x)*(-y)*(-1. + z)*z;
198 output.access(21, 2) = (1. - x)*(1. + x)*(1. - y)*(1. + y)*(-0.5 + z);
200 output.access(22, 0) = -x*(1. - y)*(1. + y)*z*(1. + z);
201 output.access(22, 1) = (1. - x)*(1. + x)*(-y)*z*(1. + z);
202 output.access(22, 2) = (1. - x)*(1. + x)*(1. - y)*(1. + y)*(0.5 + z);
204 output.access(23, 0) = (-0.5 + x)*(1. - y)*(1. + y)*(1. - z)*(1. + z);
205 output.access(23, 1) = (-1. + x)*x*(-y)*(1. - z)*(1. + z);
206 output.access(23, 2) = (-1. + x)*x*(1. - y)*(1. + y)*(-z);
208 output.access(24, 0) = (0.5 + x)*(1. - y)*(1. + y)*(1. - z)*(1. + z);
209 output.access(24, 1) = x*(1. + x)*(-y)*(1. - z)*(1. + z);
210 output.access(24, 2) = x*(1. + x)*(1. - y)*(1. + y)*(-z);
212 output.access(25, 0) = -x*(-1. + y)*y*(1. - z)*(1. + z);
213 output.access(25, 1) = (1. - x)*(1. + x)*(-0.5 + y)*(1. - z)*(1. + z);
214 output.access(25, 2) = (1. - x)*(1. + x)*(-1. + y)*y*(-z);
216 output.access(26, 0) = -x*y*(1. + y)*(1. - z)*(1. + z);
217 output.access(26, 1) = (1. - x)*(1. + x)*(0.5 + y)*(1. - z)*(1. + z);
218 output.access(26, 2) = (1. - x)*(1. + x)*y*(1. + y)*(-z);
223 const auto x = input(0);
224 const auto y = input(1);
225 const auto z = input(2);
228 output.access(0, 0) = 0.25*(-1. + y)*y*(-1. + z)*z;
229 output.access(0, 1) = (-0.125 + y*(0.25 - 0.25*z) + x*(0.25 + y*(-0.5 + 0.5*z) - 0.25*z) + 0.125*z)*z;
230 output.access(0, 2) = y*(-0.125 + x*(0.25 + y*(-0.25 + 0.5*z) - 0.5*z) + y*(0.125 - 0.25*z) + 0.25*z);
231 output.access(0, 3) = 0.25*(-1. + x)*x*(-1. + z)*z;
232 output.access(0, 4) = x*(-0.125 + y*(0.25 - 0.5*z) + x*(0.125 + y*(-0.25 + 0.5*z) - 0.25*z) + 0.25*z);
233 output.access(0, 5) = 0.25*(-1. + x)*x*(-1. + y)*y;
235 output.access(1, 0) = 0.25*(-1. + y)*y*(-1. + z)*z;
236 output.access(1, 1) = (0.125 + x*(0.25 + y*(-0.5 + 0.5*z) - 0.25*z) + y*(-0.25 + 0.25*z) - 0.125*z)*z;
237 output.access(1, 2) = y*(0.125 + x*(0.25 + y*(-0.25 + 0.5*z) - 0.5*z) + y*(-0.125 + 0.25*z) - 0.25*z);
238 output.access(1, 3) = 0.25*x*(1 + x)*(-1. + z)*z;
239 output.access(1, 4) = x*(1. + x)*(0.125 + y*(-0.25 + 0.5*z) - 0.25*z);
240 output.access(1, 5) = 0.25*x*(1 + x)*(-1. + y)*y;
242 output.access(2, 0) = 0.25*y*(1 + y)*(-1. + z)*z;
243 output.access(2, 1) = (0.125 + x*(0.25 + 0.5*y) + 0.25*y)*(-1. + z)*z;
244 output.access(2, 2) = y*(1. + y)*(-0.125 + x*(-0.25 + 0.5*z) + 0.25*z);
245 output.access(2, 3) = 0.25*x*(1 + x)*(-1. + z)*z;
246 output.access(2, 4) = x*(1. + x)*(-0.125 + y*(-0.25 + 0.5*z) + 0.25*z);
247 output.access(2, 5) = 0.25*x*(1 + x)*y*(1 + y);
249 output.access(3, 0) = 0.25*y*(1 + y)*(-1. + z)*z;
250 output.access(3, 1) = (0.125 + y*(0.25 - 0.25*z) + x*(-0.25 + y*(-0.5 + 0.5*z) + 0.25*z) - 0.125*z)*z;
251 output.access(3, 2) = y*(1. + y)*(0.125 + x*(-0.25 + 0.5*z) - 0.25*z);
252 output.access(3, 3) = 0.25*(-1. + x)*x*(-1. + z)*z;
253 output.access(3, 4) = x*(0.125 + y*(0.25 - 0.5*z) + x*(-0.125 + y*(-0.25 + 0.5*z) + 0.25*z) - 0.25*z);
254 output.access(3, 5) = 0.25*(-1. + x)*x*y*(1 + y);
256 output.access(4, 0) = 0.25*(-1. + y)*y*z*(1 + z);
257 output.access(4, 1) = (0.125 + x*(-0.25 + 0.5*y) - 0.25*y)*z*(1. + z);
258 output.access(4, 2) = y*(0.125 + x*(-0.25 + y*(0.25 + 0.5*z) - 0.5*z) + y*(-0.125 - 0.25*z) + 0.25*z);
259 output.access(4, 3) = 0.25*(-1. + x)*x*z*(1 + z);
260 output.access(4, 4) = x*(0.125 + y*(-0.25 - 0.5*z) + x*(-0.125 + y*(0.25 + 0.5*z) - 0.25*z) + 0.25*z);
261 output.access(4, 5) = 0.25*(-1. + x)*x*(-1. + y)*y;
263 output.access(5, 0) = 0.25*(-1. + y)*y*z*(1 + z);
264 output.access(5, 1) = (-0.125 + x*(-0.25 + 0.5*y) + 0.25*y)*z*(1. + z);
265 output.access(5, 2) = (-1. + y)*y*(0.125 + x*(0.25 + 0.5*z) + 0.25*z);
266 output.access(5, 3) = 0.25*x*(1 + x)*z*(1 + z);
267 output.access(5, 4) = x*(1. + x)*(-0.125 + y*(0.25 + 0.5*z) - 0.25*z);
268 output.access(5, 5) = 0.25*x*(1 + x)*(-1. + y)*y;
270 output.access(6, 0) = 0.25*y*(1 + y)*z*(1 + z);
271 output.access(6, 1) = (0.125 + x*(0.25 + 0.5*y) + 0.25*y)*z*(1. + z);
272 output.access(6, 2) = y*(1. + y)*(0.125 + x*(0.25 + 0.5*z) + 0.25*z);
273 output.access(6, 3) = 0.25*x*(1 + x)*z*(1 + z);
274 output.access(6, 4) = x*(1. + x)*(0.125 + y*(0.25 + 0.5*z) + 0.25*z);
275 output.access(6, 5) = 0.25*x*(1 + x)*y*(1 + y);
277 output.access(7, 0) = 0.25*y*(1 + y)*z*(1 + z);
278 output.access(7, 1) = (-0.125 + x*(0.25 + 0.5*y) - 0.25*y)*z*(1. + z);
279 output.access(7, 2) = y*(1. + y)*(-0.125 + x*(0.25 + 0.5*z) - 0.25*z);
280 output.access(7, 3) = 0.25*(-1. + x)*x*z*(1 + z);
281 output.access(7, 4) = (-1. + x)*x*(0.125 + y*(0.25 + 0.5*z) + 0.25*z);
282 output.access(7, 5) = 0.25*(-1. + x)*x*y*(1 + y);
284 output.access(8, 0) = -0.5*(-1. + y)*y*(-1. + z)*z;
285 output.access(8, 1) = (0. + x*(-0.5 + y))*z + (x*(0.5 - y) )*(z*z);
286 output.access(8, 2) = (y*y)*(x*(0.5 - z) ) + y*(x*(-0.5 + z));
287 output.access(8, 3) = 0.5*(1. - x)*(1. + x)*(-1. + z)*z;
288 output.access(8, 4) = 0.25 + (x*x)*(-0.25 + y*(0.5 - z) + 0.5*z) - 0.5*z + y*(-0.5 + z);
289 output.access(8, 5) = 0.5*(1. - x)*(1. + x)*(-1. + y)*y;
291 output.access(9, 0) = 0.5*(1. - y)*(1. + y)*(-1. + z)*z;
292 output.access(9, 1) = (0.5*y + x*(y))*z + (x*(-y) - 0.5*y)*(z*z);
293 output.access(9, 2) = -0.25 + (y*y)*(0.25 - 0.5*z) + 0.5*z + x*(-0.5 + (y*y)*(0.5 - z) + z);
294 output.access(9, 3) = -0.5*x*(1 + x)*(-1. + z)*z;
295 output.access(9, 4) = x*(y*(0.5 - z) ) + (x*x)*(y*(0.5 - z) );
296 output.access(9, 5) = 0.5*x*(1 + x)*(1. - y)*(1. + y);
298 output.access(10, 0) = -0.5*y*(1 + y)*(-1. + z)*z;
299 output.access(10, 1) = (0. + x*(0.5 + y))*z + (x*(-0.5 - y) )*(z*z);
300 output.access(10, 2) = y*(x*(0.5 - z) ) + (y*y)*(x*(0.5 - z) );
301 output.access(10, 3) = 0.5*(1. - x)*(1. + x)*(-1. + z)*z;
302 output.access(10, 4) = -0.25 + (x*x)*(0.25 + y*(0.5 - z) - 0.5*z) + 0.5*z + y*(-0.5 + z);
303 output.access(10, 5) = 0.5*(1. - x)*(1. + x)*y*(1 + y);
305 output.access(11, 0) = 0.5*(1. - y)*(1. + y)*(-1. + z)*z;
306 output.access(11, 1) = (-0.5*y + x*(y))*z + (x*(-y) + 0.5*y)*(z*z);
307 output.access(11, 2) = 0.25 + (y*y)*(-0.25 + 0.5*z) - 0.5*z + x*(-0.5 + (y*y)*(0.5 - z) + z);
308 output.access(11, 3) = -0.5*(-1. + x)*x*(-1. + z)*z;
309 output.access(11, 4) = (x*x)*(y*(0.5 - z) ) + x*(y*(-0.5 + z));
310 output.access(11, 5) = 0.5*(-1. + x)*x*(1. - y)*(1. + y);
312 output.access(12, 0) = 0.5*(-1. + y)*y*(1. - z)*(1. + z);
313 output.access(12, 1) = 0.25 - 0.25*(z*z) + y*(-0.5 + 0.5*(z*z)) + x*(-0.5 + 0.5*(z*z) + y*(1. - (z*z)));
314 output.access(12, 2) = (y*y)*(x*(-z) + 0.5*z) + y*(-0.5*z + x*(z));
315 output.access(12, 3) = 0.5*(-1. + x)*x*(1. - z)*(1. + z);
316 output.access(12, 4) = (x*x)*(y*(-z) + 0.5*z) + x*(-0.5*z + y*(z));
317 output.access(12, 5) = -0.5*(-1. + x)*x*(-1. + y)*y;
319 output.access(13, 0) = 0.5*(-1. + y)*y*(1. - z)*(1. + z);
320 output.access(13, 1) = -0.25 + 0.25*(z*z) + y*(0.5 - 0.5*(z*z)) + x*(-0.5 + 0.5*(z*z) + y*(1. - (z*z)));
321 output.access(13, 2) = (y*y)*(x*(-z) - 0.5*z) + y*(0.5*z + x*(z));
322 output.access(13, 3) = 0.5*x*(1 + x)*(1. - z)*(1. + z);
323 output.access(13, 4) = x*(y*(-z) + 0.5*z) + (x*x)*(y*(-z) + 0.5*z);
324 output.access(13, 5) = -0.5*x*(1 + x)*(-1. + y)*y;
326 output.access(14, 0) = 0.5*y*(1 + y)*(1. - z)*(1. + z);
327 output.access(14, 1) = 0.25 - 0.25*(z*z) + y*(0.5 - 0.5*(z*z)) + x*(0.5 - 0.5*(z*z) + y*(1. - (z*z)));
328 output.access(14, 2) = y*(x*(-z) - 0.5*z) + (y*y)*(x*(-z) - 0.5*z);
329 output.access(14, 3) = 0.5*x*(1 + x)*(1. - z)*(1. + z);
330 output.access(14, 4) = x*(y*(-z) - 0.5*z) + (x*x)*(y*(-z) - 0.5*z);
331 output.access(14, 5) = -0.5*x*(1 + x)*y*(1 + y);
333 output.access(15, 0) = 0.5*y*(1 + y)*(1. - z)*(1. + z);
334 output.access(15, 1) = -0.25 + 0.25*(z*z) + y*(-0.5 + 0.5*(z*z)) + x*(0.5 - 0.5*(z*z) + y*(1. - (z*z)));
335 output.access(15, 2) = y*(x*(-z) + 0.5*z) + (y*y)*(x*(-z) + 0.5*z);
336 output.access(15, 3) = 0.5*(-1. + x)*x*(1. - z)*(1. + z);
337 output.access(15, 4) = (x*x)*(y*(-z) - 0.5*z) + x*(0.5*z + y*(z));
338 output.access(15, 5) = -0.5*(-1. + x)*x*y*(1 + y);
340 output.access(16, 0) = -0.5*(-1. + y)*y*z*(1 + z);
341 output.access(16, 1) = (x*(0.5 - y) )*z + (x*(0.5 - y) )*(z*z);
342 output.access(16, 2) = (y*y)*(x*(-0.5 - z) ) + y*(x*(0.5 + z));
343 output.access(16, 3) = 0.5*(1. - x)*(1. + x)*z*(1 + z);
344 output.access(16, 4) = -0.25 + (x*x)*(0.25 + y*(-0.5 - z) + 0.5*z) - 0.5*z + y*(0.5 + z);
345 output.access(16, 5) = 0.5*(1. - x)*(1. + x)*(-1. + y)*y;
347 output.access(17, 0) = 0.5*(1. - y)*(1. + y)*z*(1 + z);
348 output.access(17, 1) = (x*(-y) - 0.5*y)*z + (x*(-y) - 0.5*y)*(z*z);
349 output.access(17, 2) = 0.25 + (y*y)*(-0.25 - 0.5*z) + 0.5*z + x*(0.5 + (y*y)*(-0.5 - z) + z);
350 output.access(17, 3) = -0.5*x*(1 + x)*z*(1 + z);
351 output.access(17, 4) = x*(y*(-0.5 - z) ) + (x*x)*(y*(-0.5 - z) );
352 output.access(17, 5) = 0.5*x*(1 + x)*(1. - y)*(1. + y);
354 output.access(18, 0) = -0.5*y*(1 + y)*z*(1 + z);
355 output.access(18, 1) = (x*(-0.5 - y) )*z + (x*(-0.5 - y) )*(z*z);
356 output.access(18, 2) = y*(x*(-0.5 - z) ) + (y*y)*(x*(-0.5 - z) );
357 output.access(18, 3) = 0.5*(1. - x)*(1. + x)*z*(1 + z);
358 output.access(18, 4) = 0.25 + (x*x)*(-0.25 + y*(-0.5 - z) - 0.5*z) + 0.5*z + y*(0.5 + z);
359 output.access(18, 5) = 0.5*(1. - x)*(1. + x)*y*(1 + y);
361 output.access(19, 0) = 0.5*(1. - y)*(1. + y)*z*(1 + z);
362 output.access(19, 1) = (x*(-y) + 0.5*y)*z + (x*(-y) + 0.5*y)*(z*z);
363 output.access(19, 2) = -0.25 + (y*y)*(0.25 + 0.5*z) - 0.5*z + x*(0.5 + (y*y)*(-0.5 - z) + z);
364 output.access(19, 3) = -0.5*(-1. + x)*x*z*(1 + z);
365 output.access(19, 4) = (x*x)*(y*(-0.5 - z) ) + x*(y*(0.5 + z));
366 output.access(19, 5) = 0.5*(-1. + x)*x*(1. - y)*(1. + y);
368 if constexpr (!serendipity) {
369 output.access(20, 0) = -2.*(1. - y)*(1. + y)*(1. - z)*(1. + z);
370 output.access(20, 1) = -4.*x*y*(-1. + z*z);
371 output.access(20, 2) = x*((y*y)*(-4.*z) + 4.*z);
372 output.access(20, 3) = -2.*(1. - x)*(1. + x)*(1. - z)*(1. + z);
373 output.access(20, 4) = (x*x)*(y*(-4.*z) ) + y*(4.*z);
374 output.access(20, 5) = -2.*(1. - x)*(1. + x)*(1. - y)*(1. + y);
376 output.access(21, 0) = -(1. - y)*(1. + y)*(-1. + z)*z;
377 output.access(21, 1) = (x*(-2.*y) )*z + (0. + x*(2.*y))*(z*z);
378 output.access(21, 2) = x*(1. - 2.*z + (y*y)*(-1. + 2.*z));
379 output.access(21, 3) = -(1. - x)*(1. + x)*(-1. + z)*z;
380 output.access(21, 4) = y*(1. - 2.*z) + (x*x)*(y*(-1. + 2.*z));
381 output.access(21, 5) = (1. - x)*(1. + x)*(1. - y)*(1. + y);
383 output.access(22, 0) = -(1. - y)*(1. + y)*z*(1 + z);
384 output.access(22, 1) = (0. + x*(2.*y))*z + (0. + x*(2.*y))*(z*z);
385 output.access(22, 2) = x*(-1. - 2.*z + (y*y)*(1. + 2.*z));
386 output.access(22, 3) = -(1. - x)*(1. + x)*z*(1 + z);
387 output.access(22, 4) = y*(-1. - 2.*z) + (x*x)*(y*(1. + 2.*z));
388 output.access(22, 5) = (1. - x)*(1. + x)*(1. - y)*(1. + y);
390 output.access(23, 0) = (1. - y)*(1. + y)*(1. - z)*(1. + z);
391 output.access(23, 1) = (-1. + 2.*x)*y*(-1. + z*z);
392 output.access(23, 2) = (-1. + 2.*x)*(-1. + y*y)*z;
393 output.access(23, 3) =-(-1. + x)*x*(1. - z)*(1. + z);
394 output.access(23, 4) = 2.*(-1. + x)*x*y*z;
395 output.access(23, 5) =-(-1. + x)*x*(1. - y)*(1. + y);
397 output.access(24, 0) = (1. - y)*(1. + y)*(1. - z)*(1. + z);
398 output.access(24, 1) = (1. + 2.*x)*y*(-1. + z*z);
399 output.access(24, 2) = (1. + 2.*x)*(-1. + y*y)*z;
400 output.access(24, 3) = x*(1. + x)*(-1. + z)*(1. + z);
401 output.access(24, 4) = 2.*x*(1. + x)*y*z;
402 output.access(24, 5) = x*(1. + x)*(-1. + y)*(1. + y);
404 output.access(25, 0) = -(-1. + y)*y*(1. - z)*(1. + z);
405 output.access(25, 1) = x*(-1. + 2.*y)*(-1. + z*z);
406 output.access(25, 2) = 2.*x*(-1. + y)*y*z;
407 output.access(25, 3) = (1. - x)*(1. + x)*(1. - z)*(1. + z);
408 output.access(25, 4) = (-1. + x*x)*(-1. + 2.*y)*z;
409 output.access(25, 5) =-(1. - x)*(1. + x)*(-1. + y)*y;
411 output.access(26, 0) = y*(1. + y)*(-1. + z)*(1. + z);
412 output.access(26, 1) = x*(1. + 2.*y)*(-1. + z*z);
413 output.access(26, 2) = 2.*x*y*(1. + y)*z;
414 output.access(26, 3) = (-1. + x)*(1. + x)*(-1. + z)*(1. + z);
415 output.access(26, 4) = (-1. + x*x)*(1. + 2.*y)*z;
416 output.access(26, 5) = (-1. + x)*(1. + x)*y*(1. + y);
421 const auto x = input(0);
422 const auto y = input(1);
423 const auto z = input(2);
425 output.access(0, 0) = 0.;
426 output.access(0, 1) = ((-1.+ 2.*y)*(-1.+ z)*z)/4.;
427 output.access(0, 2) = ((-1.+ y)*y*(-1.+ 2.*z))/4.;
428 output.access(0, 3) = ((-1.+ 2.*x)*(-1.+ z)*z)/4.;
429 output.access(0, 4) = ((-1.+ 2.*x)*(-1.+ 2.*y)*(-1.+ 2.*z))/8.;
430 output.access(0, 5) = ((-1.+ 2.*x)*(-1.+ y)*y)/4.;
431 output.access(0, 6) = 0.;
432 output.access(0, 7) = ((-1.+ x)*x*(-1.+ 2.*z))/4.;
433 output.access(0, 8) = ((-1.+ x)*x*(-1.+ 2.*y))/4.;
434 output.access(0, 9) = 0.;
436 output.access(1, 0) = 0.;
437 output.access(1, 1) = ((-1.+ 2.*y)*(-1.+ z)*z)/4.;
438 output.access(1, 2) = ((-1.+ y)*y*(-1.+ 2.*z))/4.;
439 output.access(1, 3) = ((1.+ 2.*x)*(-1.+ z)*z)/4.;
440 output.access(1, 4) = ((1.+ 2.*x)*(-1.+ 2.*y)*(-1.+ 2.*z))/8.;
441 output.access(1, 5) = ((1.+ 2.*x)*(-1.+ y)*y)/4.;
442 output.access(1, 6) = 0.;
443 output.access(1, 7) = (x*(1.+ x)*(-1.+ 2.*z))/4.;
444 output.access(1, 8) = (x*(1.+ x)*(-1.+ 2.*y))/4.;
445 output.access(1, 9) = 0.;
447 output.access(2, 0) = 0.;
448 output.access(2, 1) = ((1.+ 2.*y)*(-1.+ z)*z)/4.;
449 output.access(2, 2) = (y*(1.+ y)*(-1.+ 2.*z))/4.;
450 output.access(2, 3) = ((1.+ 2.*x)*(-1.+ z)*z)/4.;
451 output.access(2, 4) = ((1.+ 2.*x)*(1.+ 2.*y)*(-1.+ 2.*z))/8.;
452 output.access(2, 5) = ((1.+ 2.*x)*y*(1.+ y))/4.;
453 output.access(2, 6) = 0.;
454 output.access(2, 7) = (x*(1.+ x)*(-1.+ 2.*z))/4.;
455 output.access(2, 8) = (x*(1.+ x)*(1.+ 2.*y))/4.;
456 output.access(2, 9) = 0.;
458 output.access(3, 0) = 0.;
459 output.access(3, 1) = ((1.+ 2.*y)*(-1.+ z)*z)/4.;
460 output.access(3, 2) = (y*(1.+ y)*(-1.+ 2.*z))/4.;
461 output.access(3, 3) = ((-1.+ 2.*x)*(-1.+ z)*z)/4.;
462 output.access(3, 4) = ((-1.+ 2.*x)*(1.+ 2.*y)*(-1.+ 2.*z))/8.;
463 output.access(3, 5) = ((-1.+ 2.*x)*y*(1.+ y))/4.;
464 output.access(3, 6) = 0.;
465 output.access(3, 7) = ((-1.+ x)*x*(-1.+ 2.*z))/4.;
466 output.access(3, 8) = ((-1.+ x)*x*(1.+ 2.*y))/4.;
467 output.access(3, 9) = 0.;
469 output.access(4, 0) = 0.;
470 output.access(4, 1) = ((-1.+ 2.*y)*z*(1.+ z))/4.;
471 output.access(4, 2) = ((-1.+ y)*y*(1.+ 2.*z))/4.;
472 output.access(4, 3) = ((-1.+ 2.*x)*z*(1.+ z))/4.;
473 output.access(4, 4) = ((-1.+ 2.*x)*(-1.+ 2.*y)*(1.+ 2.*z))/8.;
474 output.access(4, 5) = ((-1.+ 2.*x)*(-1.+ y)*y)/4.;
475 output.access(4, 6) = 0.;
476 output.access(4, 7) = ((-1.+ x)*x*(1.+ 2.*z))/4.;
477 output.access(4, 8) = ((-1.+ x)*x*(-1.+ 2.*y))/4.;
478 output.access(4, 9) = 0.;
480 output.access(5, 0) = 0.;
481 output.access(5, 1) = ((-1.+ 2.*y)*z*(1.+ z))/4.;
482 output.access(5, 2) = ((-1.+ y)*y*(1.+ 2.*z))/4.;
483 output.access(5, 3) = ((1.+ 2.*x)*z*(1.+ z))/4.;
484 output.access(5, 4) = ((1.+ 2.*x)*(-1.+ 2.*y)*(1.+ 2.*z))/8.;
485 output.access(5, 5) = ((1.+ 2.*x)*(-1.+ y)*y)/4.;
486 output.access(5, 6) = 0.;
487 output.access(5, 7) = (x*(1.+ x)*(1.+ 2.*z))/4.;
488 output.access(5, 8) = (x*(1.+ x)*(-1.+ 2.*y))/4.;
489 output.access(5, 9) = 0.;
491 output.access(6, 0) = 0.;
492 output.access(6, 1) = ((1.+ 2.*y)*z*(1.+ z))/4.;
493 output.access(6, 2) = (y*(1.+ y)*(1.+ 2.*z))/4.;
494 output.access(6, 3) = ((1.+ 2.*x)*z*(1.+ z))/4.;
495 output.access(6, 4) = ((1.+ 2.*x)*(1.+ 2.*y)*(1.+ 2.*z))/8.;
496 output.access(6, 5) = ((1.+ 2.*x)*y*(1.+ y))/4.;
497 output.access(6, 6) = 0.;
498 output.access(6, 7) = (x*(1.+ x)*(1.+ 2.*z))/4.;
499 output.access(6, 8) = (x*(1.+ x)*(1.+ 2.*y))/4.;
500 output.access(6, 9) = 0.;
502 output.access(7, 0) = 0.;
503 output.access(7, 1) = ((1.+ 2.*y)*z*(1.+ z))/4.;
504 output.access(7, 2) = (y*(1.+ y)*(1.+ 2.*z))/4.;
505 output.access(7, 3) = ((-1.+ 2.*x)*z*(1.+ z))/4.;
506 output.access(7, 4) = ((-1.+ 2.*x)*(1.+ 2.*y)*(1.+ 2.*z))/8.;
507 output.access(7, 5) = ((-1.+ 2.*x)*y*(1.+ y))/4.;
508 output.access(7, 6) = 0.;
509 output.access(7, 7) = ((-1.+ x)*x*(1.+ 2.*z))/4.;
510 output.access(7, 8) = ((-1.+ x)*x*(1.+ 2.*y))/4.;
511 output.access(7, 9) = 0.;
513 output.access(8, 0) = 0.;
514 output.access(8, 1) = -((-1.+ 2.*y)*(-1.+ z)*z)/2.;
515 output.access(8, 2) = -((-1.+ y)*y*(-1.+ 2.*z))/2.;
516 output.access(8, 3) = -(x*(-1.+ z)*z);
517 output.access(8, 4) = -(x*(-1.+ 2.*y)*(-1.+ 2.*z))/2.;
518 output.access(8, 5) = -(x*(-1.+ y)*y);
519 output.access(8, 6) = 0.;
520 output.access(8, 7) = -((-1.+ (x*x))*(-1.+ 2.*z))/2.;
521 output.access(8, 8) = -((-1.+ (x*x))*(-1.+ 2.*y))/2.;
522 output.access(8, 9) = 0.;
524 output.access(9, 0) = 0.;
525 output.access(9, 1) = -(y*(-1.+ z)*z);
526 output.access(9, 2) = -((-1.+ (y*y))*(-1.+ 2.*z))/2.;
527 output.access(9, 3) = -((1.+ 2.*x)*(-1.+ z)*z)/2.;
528 output.access(9, 4) = -((1.+ 2.*x)*y*(-1.+ 2.*z))/2.;
529 output.access(9, 5) = -((1.+ 2.*x)*(-1.+ (y*y)))/2.;
530 output.access(9, 6) = 0.;
531 output.access(9, 7) = -(x*(1.+ x)*(-1.+ 2.*z))/2.;
532 output.access(9, 8) = -(x*(1.+ x)*y);
533 output.access(9, 9) = 0.;
535 output.access(10, 0) = 0.;
536 output.access(10, 1) = -((1.+ 2.*y)*(-1.+ z)*z)/2.;
537 output.access(10, 2) = -(y*(1.+ y)*(-1.+ 2.*z))/2.;
538 output.access(10, 3) = -(x*(-1.+ z)*z);
539 output.access(10, 4) = -(x*(1.+ 2.*y)*(-1.+ 2.*z))/2.;
540 output.access(10, 5) = -(x*y*(1.+ y));
541 output.access(10, 6) = 0.;
542 output.access(10, 7) = -((-1.+ (x*x))*(-1.+ 2.*z))/2.;
543 output.access(10, 8) = -((-1.+ (x*x))*(1.+ 2.*y))/2.;
544 output.access(10, 9) = 0.;
546 output.access(11, 0) = 0.;
547 output.access(11, 1) = -(y*(-1.+ z)*z);
548 output.access(11, 2) = -((-1.+ (y*y))*(-1.+ 2.*z))/2.;
549 output.access(11, 3) = -((-1.+ 2.*x)*(-1.+ z)*z)/2.;
550 output.access(11, 4) = -((-1.+ 2.*x)*y*(-1.+ 2.*z))/2.;
551 output.access(11, 5) = -((-1.+ 2.*x)*(-1.+ (y*y)))/2.;
552 output.access(11, 6) = 0.;
553 output.access(11, 7) = -((-1.+ x)*x*(-1.+ 2.*z))/2.;
554 output.access(11, 8) = -((-1.+ x)*x*y);
555 output.access(11, 9) = 0.;
557 output.access(12, 0) = 0.;
558 output.access(12, 1) = -((-1.+ 2.*y)*(-1.+ (z*z)))/2.;
559 output.access(12, 2) = -((-1.+ y)*y*z);
560 output.access(12, 3) = -((-1.+ 2.*x)*(-1.+ (z*z)))/2.;
561 output.access(12, 4) = -((-1.+ 2.*x)*(-1.+ 2.*y)*z)/2.;
562 output.access(12, 5) = -((-1.+ 2.*x)*(-1.+ y)*y)/2.;
563 output.access(12, 6) = 0.;
564 output.access(12, 7) = -((-1.+ x)*x*z);
565 output.access(12, 8) = -((-1.+ x)*x*(-1.+ 2.*y))/2.;
566 output.access(12, 9) = 0.;
568 output.access(13, 0) = 0.;
569 output.access(13, 1) = -((-1.+ 2.*y)*(-1.+ (z*z)))/2.;
570 output.access(13, 2) = -((-1.+ y)*y*z);
571 output.access(13, 3) = -((1.+ 2.*x)*(-1.+ (z*z)))/2.;
572 output.access(13, 4) = -((1.+ 2.*x)*(-1.+ 2.*y)*z)/2.;
573 output.access(13, 5) = -((1.+ 2.*x)*(-1.+ y)*y)/2.;
574 output.access(13, 6) = 0.;
575 output.access(13, 7) = -(x*(1.+ x)*z);
576 output.access(13, 8) = -(x*(1.+ x)*(-1.+ 2.*y))/2.;
577 output.access(13, 9) = 0.;
579 output.access(14, 0) = 0.;
580 output.access(14, 1) = -((1.+ 2.*y)*(-1.+ (z*z)))/2.;
581 output.access(14, 2) = -(y*(1.+ y)*z);
582 output.access(14, 3) = -((1.+ 2.*x)*(-1.+ (z*z)))/2.;
583 output.access(14, 4) = -((1.+ 2.*x)*(1.+ 2.*y)*z)/2.;
584 output.access(14, 5) = -((1.+ 2.*x)*y*(1.+ y))/2.;
585 output.access(14, 6) = 0.;
586 output.access(14, 7) = -(x*(1.+ x)*z);
587 output.access(14, 8) = -(x*(1.+ x)*(1.+ 2.*y))/2.;
588 output.access(14, 9) = 0.;
590 output.access(15, 0) = 0.;
591 output.access(15, 1) = -((1.+ 2.*y)*(-1.+ (z*z)))/2.;
592 output.access(15, 2) = -(y*(1.+ y)*z);
593 output.access(15, 3) = -((-1.+ 2.*x)*(-1.+ (z*z)))/2.;
594 output.access(15, 4) = -((-1.+ 2.*x)*(1.+ 2.*y)*z)/2.;
595 output.access(15, 5) = -((-1.+ 2.*x)*y*(1.+ y))/2.;
596 output.access(15, 6) = 0.;
597 output.access(15, 7) = -((-1.+ x)*x*z);
598 output.access(15, 8) = -((-1.+ x)*x*(1.+ 2.*y))/2.;
599 output.access(15, 9) = 0.;
601 output.access(16, 0) = 0.;
602 output.access(16, 1) = -((-1.+ 2.*y)*z*(1.+ z))/2.;
603 output.access(16, 2) = -((-1.+ y)*y*(1.+ 2.*z))/2.;
604 output.access(16, 3) = -(x*z*(1.+ z));
605 output.access(16, 4) = -(x*(-1.+ 2.*y)*(1.+ 2.*z))/2.;
606 output.access(16, 5) = -(x*(-1.+ y)*y);
607 output.access(16, 6) = 0.;
608 output.access(16, 7) = -((-1.+ (x*x))*(1.+ 2.*z))/2.;
609 output.access(16, 8) = -((-1.+ (x*x))*(-1.+ 2.*y))/2.;
610 output.access(16, 9) = 0.;
612 output.access(17, 0) = 0.;
613 output.access(17, 1) = -(y*z*(1.+ z));
614 output.access(17, 2) = -((-1.+ (y*y))*(1.+ 2.*z))/2.;
615 output.access(17, 3) = -((1.+ 2.*x)*z*(1.+ z))/2.;
616 output.access(17, 4) = -((1.+ 2.*x)*y*(1.+ 2.*z))/2.;
617 output.access(17, 5) = -((1.+ 2.*x)*(-1.+ (y*y)))/2.;
618 output.access(17, 6) = 0.;
619 output.access(17, 7) = -(x*(1.+ x)*(1.+ 2.*z))/2.;
620 output.access(17, 8) = -(x*(1.+ x)*y);
621 output.access(17, 9) = 0.;
623 output.access(18, 0) = 0.;
624 output.access(18, 1) = -((1.+ 2.*y)*z*(1.+ z))/2.;
625 output.access(18, 2) = -(y*(1.+ y)*(1.+ 2.*z))/2.;
626 output.access(18, 3) = -(x*z*(1.+ z));
627 output.access(18, 4) = -(x*(1.+ 2.*y)*(1.+ 2.*z))/2.;
628 output.access(18, 5) = -(x*y*(1.+ y));
629 output.access(18, 6) = 0.;
630 output.access(18, 7) = -((-1.+ (x*x))*(1.+ 2.*z))/2.;
631 output.access(18, 8) = -((-1.+ (x*x))*(1.+ 2.*y))/2.;
632 output.access(18, 9) = 0.;
634 output.access(19, 0) = 0.;
635 output.access(19, 1) = -(y*z*(1.+ z));
636 output.access(19, 2) = -((-1.+ (y*y))*(1.+ 2.*z))/2.;
637 output.access(19, 3) = -((-1.+ 2.*x)*z*(1.+ z))/2.;
638 output.access(19, 4) = -((-1.+ 2.*x)*y*(1.+ 2.*z))/2.;
639 output.access(19, 5) = -((-1.+ 2.*x)*(-1.+ (y*y)))/2.;
640 output.access(19, 6) = 0.;
641 output.access(19, 7) = -((-1.+ x)*x*(1.+ 2.*z))/2.;
642 output.access(19, 8) = -((-1.+ x)*x*y);
643 output.access(19, 9) = 0.;
645 if constexpr(!serendipity) {
646 output.access(20, 0) = 0.;
647 output.access(20, 1) = -4*y*(-1.+ (z*z));
648 output.access(20, 2) = -4*(-1.+ (y*y))*z;
649 output.access(20, 3) = -4*x*(-1.+ (z*z));
650 output.access(20, 4) = -8*x*y*z;
651 output.access(20, 5) = -4*x*(-1.+ (y*y));
652 output.access(20, 6) = 0.;
653 output.access(20, 7) = -4*(-1.+ (x*x))*z;
654 output.access(20, 8) = -4*(-1.+ (x*x))*y;
655 output.access(20, 9) = 0.;
657 output.access(21, 0) = 0.;
658 output.access(21, 1) = 2.*y*(-1.+ z)*z;
659 output.access(21, 2) = (-1.+ (y*y))*(-1.+ 2.*z);
660 output.access(21, 3) = 2.*x*(-1.+ z)*z;
661 output.access(21, 4) = 2.*x*y*(-1.+ 2.*z);
662 output.access(21, 5) = 2.*x*(-1.+ (y*y));
663 output.access(21, 6) = 0.;
664 output.access(21, 7) = (-1.+ (x*x))*(-1.+ 2.*z);
665 output.access(21, 8) = 2.*(-1.+ (x*x))*y;
666 output.access(21, 9) = 0.;
668 output.access(22, 0) = 0.;
669 output.access(22, 1) = 2.*y*z*(1.+ z);
670 output.access(22, 2) = (-1.+ (y*y))*(1.+ 2.*z);
671 output.access(22, 3) = 2.*x*z*(1.+ z);
672 output.access(22, 4) = 2.*x*y*(1.+ 2.*z);
673 output.access(22, 5) = 2.*x*(-1.+ (y*y));
674 output.access(22, 6) = 0.;
675 output.access(22, 7) = (-1.+ (x*x))*(1.+ 2.*z);
676 output.access(22, 8) = 2.*(-1.+ (x*x))*y;
677 output.access(22, 9) = 0.;
679 output.access(23, 0) = 0.;
680 output.access(23, 1) = 2.*y*(-1.+ (z*z));
681 output.access(23, 2) = 2.*(-1.+ (y*y))*z;
682 output.access(23, 3) = (-1.+ 2.*x)*(-1.+ (z*z));
683 output.access(23, 4) = 2.*(-1.+ 2.*x)*y*z;
684 output.access(23, 5) = (-1.+ 2.*x)*(-1.+ (y*y));
685 output.access(23, 6) = 0.;
686 output.access(23, 7) = 2.*(-1.+ x)*x*z;
687 output.access(23, 8) = 2.*(-1.+ x)*x*y;
688 output.access(23, 9) = 0.;
690 output.access(24, 0) = 0.;
691 output.access(24, 1) = 2.*y*(-1.+ (z*z));
692 output.access(24, 2) = 2.*(-1.+ (y*y))*z;
693 output.access(24, 3) = (1.+ 2.*x)*(-1.+ (z*z));
694 output.access(24, 4) = 2.*(1.+ 2.*x)*y*z;
695 output.access(24, 5) = (1.+ 2.*x)*(-1.+ (y*y));
696 output.access(24, 6) = 0.;
697 output.access(24, 7) = 2.*x*(1.+ x)*z;
698 output.access(24, 8) = 2.*x*(1.+ x)*y;
699 output.access(24, 9) = 0.;
701 output.access(25, 0) = 0.;
702 output.access(25, 1) = (-1.+ 2.*y)*(-1.+ (z*z));
703 output.access(25, 2) = 2.*(-1.+ y)*y*z;
704 output.access(25, 3) = 2.*x*(-1.+ (z*z));
705 output.access(25, 4) = 2.*x*(-1.+ 2.*y)*z;
706 output.access(25, 5) = 2.*x*(-1.+ y)*y;
707 output.access(25, 6) = 0.;
708 output.access(25, 7) = 2.*(-1.+ (x*x))*z;
709 output.access(25, 8) = (-1.+ (x*x))*(-1.+ 2.*y);
710 output.access(25, 9) = 0.;
712 output.access(26, 0) = 0.;
713 output.access(26, 1) = (1.+ 2.*y)*(-1.+ (z*z));
714 output.access(26, 2) = 2.*y*(1.+ y)*z;
715 output.access(26, 3) = 2.*x*(-1.+ (z*z));
716 output.access(26, 4) = 2.*x*(1.+ 2.*y)*z;
717 output.access(26, 5) = 2.*x*y*(1.+ y);
718 output.access(26, 6) = 0.;
719 output.access(26, 7) = 2.*(-1.+ (x*x))*z;
720 output.access(26, 8) = (-1.+ (x*x))*(1.+ 2.*y);
721 output.access(26, 9) = 0.;
728 const ordinal_type jend = output.extent(1);
729 const ordinal_type iend = output.extent(0);
731 for (ordinal_type j=0;j<jend;++j)
732 for (ordinal_type i=0;i<iend;++i)
733 output.access(i, j) = 0.0;
735 const auto x = input(0);
736 const auto y = input(1);
737 const auto z = input(2);
738 output.access(0, 3) = ((-1.+ z)*z)/2.;
739 output.access(0, 4) = ((-1.+ 2.*y)*(-1.+ 2.*z))/4.;
740 output.access(0, 5) = ((-1.+ y)*y)/2.;
741 output.access(0, 7) = ((-1.+ 2.*x)*(-1.+ 2.*z))/4.;
742 output.access(0, 8) = ((-1.+ 2.*x)*(-1.+ 2.*y))/4.;
743 output.access(0, 12)= ((-1.+ x)*x)/2.;
745 output.access(1, 3) = ((-1.+ z)*z)/2.;
746 output.access(1, 4) = ((-1.+ 2.*y)*(-1.+ 2.*z))/4.;
747 output.access(1, 5) = ((-1.+ y)*y)/2.;
748 output.access(1, 7) = ((1. + 2.*x)*(-1.+ 2.*z))/4.;
749 output.access(1, 8) = ((1. + 2.*x)*(-1.+ 2.*y))/4.;
750 output.access(1, 12)= (x*(1. + x))/2.;
752 output.access(2, 3) = ((-1.+ z)*z)/2.;
753 output.access(2, 4) = ((1. + 2.*y)*(-1.+ 2.*z))/4.;
754 output.access(2, 5) = (y*(1. + y))/2.;
755 output.access(2, 7) = ((1. + 2.*x)*(-1.+ 2.*z))/4.;
756 output.access(2, 8) = ((1. + 2.*x)*(1. + 2.*y))/4.;
757 output.access(2, 12)= (x*(1. + x))/2.;
759 output.access(3, 3) = ((-1.+ z)*z)/2.;
760 output.access(3, 4) = ((1. + 2.*y)*(-1.+ 2.*z))/4.;
761 output.access(3, 5) = (y*(1. + y))/2.;
762 output.access(3, 7) = ((-1.+ 2.*x)*(-1.+ 2.*z))/4.;
763 output.access(3, 8) = ((-1.+ 2.*x)*(1. + 2.*y))/4.;
764 output.access(3, 12)= ((-1.+ x)*x)/2.;
766 output.access(4, 3) = (z*(1. + z))/2.;
767 output.access(4, 4) = ((-1.+ 2.*y)*(1. + 2.*z))/4.;
768 output.access(4, 5) = ((-1.+ y)*y)/2.;
769 output.access(4, 7) = ((-1.+ 2.*x)*(1. + 2.*z))/4.;
770 output.access(4, 8) = ((-1.+ 2.*x)*(-1.+ 2.*y))/4.;
771 output.access(4, 12)= ((-1.+ x)*x)/2.;
773 output.access(5, 3) = (z*(1. + z))/2.;
774 output.access(5, 4) = ((-1.+ 2.*y)*(1. + 2.*z))/4.;
775 output.access(5, 5) = ((-1.+ y)*y)/2.;
776 output.access(5, 7) = ((1. + 2.*x)*(1. + 2.*z))/4.;
777 output.access(5, 8) = ((1. + 2.*x)*(-1.+ 2.*y))/4.;
778 output.access(5, 12)= (x*(1. + x))/2.;
780 output.access(6, 3) = (z*(1. + z))/2.;
781 output.access(6, 4) = ((1. + 2.*y)*(1. + 2.*z))/4.;
782 output.access(6, 5) = (y*(1. + y))/2.;
783 output.access(6, 7) = ((1. + 2.*x)*(1. + 2.*z))/4.;
784 output.access(6, 8) = ((1. + 2.*x)*(1. + 2.*y))/4.;
785 output.access(6, 12)= (x*(1. + x))/2.;
787 output.access(7, 3) = (z*(1. + z))/2.;
788 output.access(7, 4) = ((1. + 2.*y)*(1. + 2.*z))/4.;
789 output.access(7, 5) = (y*(1. + y))/2.;
790 output.access(7, 7) = ((-1.+ 2.*x)*(1. + 2.*z))/4.;
791 output.access(7, 8) = ((-1.+ 2.*x)*(1. + 2.*y))/4.;
792 output.access(7, 12)= ((-1.+ x)*x)/2.;
794 output.access(8, 3) = -((-1.+ z)*z);
795 output.access(8, 4) = -0.5 + y + z - 2.*y*z;
796 output.access(8, 5) = -((-1.+ y)*y);
797 output.access(8, 7) = x - 2.*x*z;
798 output.access(8, 8) = x - 2.*x*y;
799 output.access(8, 12)= 1. - x*x;
801 output.access(9, 3) = -((-1.+ z)*z);
802 output.access(9, 4) = y - 2.*y*z;
803 output.access(9, 5) = 1 - y*y;
804 output.access(9, 7) = 0.5 + x - z - 2.*x*z;
805 output.access(9, 8) = -((1. + 2.*x)*y);
806 output.access(9, 12)= -(x*(1. + x));
808 output.access(10, 3) = -((-1.+ z)*z);
809 output.access(10, 4) = 0.5 + y - z - 2.*y*z;
810 output.access(10, 5) = -(y*(1. + y));
811 output.access(10, 7) = x - 2.*x*z;
812 output.access(10, 8) = -(x*(1. + 2.*y));
813 output.access(10, 12)= 1. - x*x;
815 output.access(11, 3) = -((-1.+ z)*z);
816 output.access(11, 4) = y - 2.*y*z;
817 output.access(11, 5) = 1. - y*y;
818 output.access(11, 7) = -0.5 + x + z - 2.*x*z;
819 output.access(11, 8) = y - 2.*x*y;
820 output.access(11, 12)= -((-1.+ x)*x);
822 output.access(12, 3) = 1. - z*z;
823 output.access(12, 4) = z - 2.*y*z;
824 output.access(12, 5) = -((-1.+ y)*y);
825 output.access(12, 7) = z - 2.*x*z;
826 output.access(12, 8) = -0.5 + x + y - 2.*x*y;
827 output.access(12, 12)= -((-1.+ x)*x);
829 output.access(13, 3) = 1. - z*z;
830 output.access(13, 4) = z - 2.*y*z;
831 output.access(13, 5) = -((-1.+ y)*y);
832 output.access(13, 7) = -((1. + 2.*x)*z);
833 output.access(13, 8) = 0.5 + x - y - 2.*x*y;
834 output.access(13, 12)= -(x*(1. + x));
836 output.access(14, 3) = 1. - z*z;
837 output.access(14, 4) = -((1. + 2.*y)*z);
838 output.access(14, 5) = -(y*(1. + y));
839 output.access(14, 7) = -((1. + 2.*x)*z);
840 output.access(14, 8) = -((1. + 2.*x)*(1. + 2.*y))/2.;
841 output.access(14, 12)= -(x*(1. + x));
843 output.access(15, 3) = 1. - z*z;
844 output.access(15, 4) = -((1. + 2.*y)*z);
845 output.access(15, 5) = -(y*(1. + y));
846 output.access(15, 7) = z - 2.*x*z;
847 output.access(15, 8) = 0.5 + y - x*(1. + 2.*y);
848 output.access(15, 12)= -((-1.+ x)*x);
850 output.access(16, 3) = -(z*(1. + z));
851 output.access(16, 4) = 0.5 + z - y*(1. + 2.*z);
852 output.access(16, 5) = -((-1.+ y)*y);
853 output.access(16, 7) = -(x*(1. + 2.*z));
854 output.access(16, 8) = x - 2.*x*y;
855 output.access(16, 12)= 1. - x*x;
857 output.access(17, 3) = -(z*(1. + z));
858 output.access(17, 4) = -(y*(1. + 2.*z));
859 output.access(17, 5) = 1. - y*y;
860 output.access(17, 7) = -((1. + 2.*x)*(1. + 2.*z))/2.;
861 output.access(17, 8) = -((1. + 2.*x)*y);
862 output.access(17, 12)= -(x*(1. + x));
864 output.access(18, 3) = -(z*(1. + z));
865 output.access(18, 4) = -((1. + 2.*y)*(1. + 2.*z))/2.;
866 output.access(18, 5) = -(y*(1. + y));
867 output.access(18, 7) = -(x*(1. + 2.*z));
868 output.access(18, 8) = -(x*(1. + 2.*y));
869 output.access(18, 12)= 1. - x*x;
871 output.access(19, 3) = -(z*(1. + z));
872 output.access(19, 4) = -(y*(1. + 2.*z));
873 output.access(19, 5) = 1. - y*y;
874 output.access(19, 7) = 0.5 + z - x*(1. + 2.*z);
875 output.access(19, 8) = y - 2.*x*y;
876 output.access(19, 12)= -((-1.+ x)*x);
878 if constexpr (!serendipity) {
879 output.access(20, 3) = 4. - 4.*z*z;
880 output.access(20, 4) = -8.*y*z;
881 output.access(20, 5) = 4. - 4.*y*y;
882 output.access(20, 7) = -8.*x*z;
883 output.access(20, 8) = -8.*x*y;
884 output.access(20, 12)= 4. - 4.*x*x;
886 output.access(21, 3) = 2.*(-1.+ z)*z;
887 output.access(21, 4) = 2.*y*(-1.+ 2.*z);
888 output.access(21, 5) = 2.*(-1.+ y*y);
889 output.access(21, 7) = 2.*x*(-1.+ 2.*z);
890 output.access(21, 8) = 4.*x*y;
891 output.access(21, 12)= 2.*(-1.+ x*x);
893 output.access(22, 3) = 2.*z*(1. + z);
894 output.access(22, 4) = 2.*(y + 2.*y*z);
895 output.access(22, 5) = 2.*(-1.+ y*y);
896 output.access(22, 7) = 2.*(x + 2.*x*z);
897 output.access(22, 8) = 4.*x*y;
898 output.access(22, 12)= 2.*(-1.+ x*x);
900 output.access(23, 3) = 2.*(-1.+ z*z);
901 output.access(23, 4) = 4.*y*z;
902 output.access(23, 5) = 2.*(-1.+ y*y);
903 output.access(23, 7) = 2.*(-1.+ 2.*x)*z;
904 output.access(23, 8) = 2.*(-1.+ 2.*x)*y;
905 output.access(23, 12)= 2.*(-1.+ x)*x;
907 output.access(24, 3) = 2.*(-1.+ z*z);
908 output.access(24, 4) = 4.*y*z;
909 output.access(24, 5) = 2.*(-1.+ y*y);
910 output.access(24, 7) = 2.*(z + 2.*x*z);
911 output.access(24, 8) = 2.*(y + 2.*x*y);
912 output.access(24, 12)= 2.*x*(1. + x);
914 output.access(25, 3) = 2.*(-1.+ z*z);
915 output.access(25, 4) = 2.*(-1.+ 2.*y)*z;
916 output.access(25, 5) = 2.*(-1.+ y)*y;
917 output.access(25, 7) = 4.*x*z;
918 output.access(25, 8) = 2.*x*(-1.+ 2.*y);
919 output.access(25, 12)= 2.*(-1.+ x*x);
921 output.access(26, 3) = 2.*(-1.+ z*z);
922 output.access(26, 4) = 2.*(z + 2.*y*z);
923 output.access(26, 5) = 2.*y*(1. + y);
924 output.access(26, 7) = 4.*x*z;
925 output.access(26, 8) = 2.*(x + 2.*x*y);
926 output.access(26, 12)= 2.*(-1.+ x*x);
930 case OPERATOR_MAX : {
931 const ordinal_type jend = output.extent(1);
932 const ordinal_type iend = output.extent(0);
934 for (ordinal_type j=0;j<jend;++j)
935 for (ordinal_type i=0;i<iend;++i)
936 output.access(i, j) = 0.0;
940 INTREPID2_TEST_FOR_ABORT( opType != OPERATOR_VALUE &&
941 opType != OPERATOR_GRAD &&
942 opType != OPERATOR_CURL &&
943 opType != OPERATOR_D1 &&
944 opType != OPERATOR_D2 &&
945 opType != OPERATOR_D3 &&
946 opType != OPERATOR_D4 &&
947 opType != OPERATOR_MAX,
948 ">>> ERROR: (Intrepid2::Basis_HGRAD_HEX_DEG2_FEM::Serial::getValues) operator is not supported");
954 template<
bool serendipity>
955 template<
typename DT,
956 typename outputValueValueType,
class ...outputValueProperties,
957 typename inputPointValueType,
class ...inputPointProperties>
960 getValues( Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
961 const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
962 const EOperator operatorType ) {
963 typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
964 typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
965 typedef typename ExecSpace<typename inputPointViewType::execution_space,typename DT::execution_space>::ExecSpaceType ExecSpaceType;
968 const auto loopSize = inputPoints.extent(0);
969 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
971 switch (operatorType) {
973 case OPERATOR_VALUE: {
974 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_VALUE> FunctorType;
975 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
980 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_GRAD> FunctorType;
981 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
985 INTREPID2_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_CURL), std::invalid_argument,
986 ">>> ERROR (Basis_HGRAD_HEX_DEG2_FEM): CURL is invalid operator for rank-0 (scalar) functions in 3D");
990 INTREPID2_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
991 ">>> ERROR (Basis_HGRAD_HEX_DEG2_FEM): DIV is invalid operator for rank-0 (scalar) functions in 3D");
995 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_D2> FunctorType;
996 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
1000 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_D3> FunctorType;
1001 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
1005 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_D4> FunctorType;
1006 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
1011 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
1012 ">>> ERROR (Basis_HGRAD_HEX_DEG2_FEM): operator not supported");
1018 case OPERATOR_D10: {
1019 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_MAX> FunctorType;
1020 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
1024 INTREPID2_TEST_FOR_EXCEPTION( !( Intrepid2::isValidOperator(operatorType) ), std::invalid_argument,
1025 ">>> ERROR (Basis_HGRAD_HEX_DEG2_FEM): Invalid operator type");
1034 template<
bool serendipity,
typename DT,
typename OT,
typename PT>
1039 this ->
basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >() );
1047 const ordinal_type tagSize = 4;
1048 const ordinal_type posScDim = 0;
1049 const ordinal_type posScOrd = 1;
1050 const ordinal_type posDfOrd = 2;
1053 ordinal_type tags[108] = { 0, 0, 0, 1,
1085 OrdinalTypeArray1DHost tagView(&tags[0], serendipity ? 80 : 108);
1098 Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
1101 dofCoords(0,0) = -1.0; dofCoords(0,1) = -1.0; dofCoords(0,2) = -1.0;
1102 dofCoords(1,0) = 1.0; dofCoords(1,1) = -1.0; dofCoords(1,2) = -1.0;
1103 dofCoords(2,0) = 1.0; dofCoords(2,1) = 1.0; dofCoords(2,2) = -1.0;
1104 dofCoords(3,0) = -1.0; dofCoords(3,1) = 1.0; dofCoords(3,2) = -1.0;
1105 dofCoords(4,0) = -1.0; dofCoords(4,1) = -1.0; dofCoords(4,2) = 1.0;
1106 dofCoords(5,0) = 1.0; dofCoords(5,1) = -1.0; dofCoords(5,2) = 1.0;
1107 dofCoords(6,0) = 1.0; dofCoords(6,1) = 1.0; dofCoords(6,2) = 1.0;
1108 dofCoords(7,0) = -1.0; dofCoords(7,1) = 1.0; dofCoords(7,2) = 1.0;
1110 dofCoords(8,0) = 0.0; dofCoords(8,1) = -1.0; dofCoords(8,2) = -1.0;
1111 dofCoords(9,0) = 1.0; dofCoords(9,1) = 0.0; dofCoords(9,2) = -1.0;
1112 dofCoords(10,0) = 0.0; dofCoords(10,1) = 1.0; dofCoords(10,2) = -1.0;
1113 dofCoords(11,0) = -1.0; dofCoords(11,1) = 0.0; dofCoords(11,2) = -1.0;
1114 dofCoords(12,0) = -1.0; dofCoords(12,1) = -1.0; dofCoords(12,2) = 0.0;
1115 dofCoords(13,0) = 1.0; dofCoords(13,1) = -1.0; dofCoords(13,2) = 0.0;
1116 dofCoords(14,0) = 1.0; dofCoords(14,1) = 1.0; dofCoords(14,2) = 0.0;
1117 dofCoords(15,0) = -1.0; dofCoords(15,1) = 1.0; dofCoords(15,2) = 0.0;
1118 dofCoords(16,0) = 0.0; dofCoords(16,1) = -1.0; dofCoords(16,2) = 1.0;
1119 dofCoords(17,0) = 1.0; dofCoords(17,1) = 0.0; dofCoords(17,2) = 1.0;
1120 dofCoords(18,0) = 0.0; dofCoords(18,1) = 1.0; dofCoords(18,2) = 1.0;
1121 dofCoords(19,0) = -1.0; dofCoords(19,1) = 0.0; dofCoords(19,2) = 1.0;
1123 if constexpr (!serendipity) {
1124 dofCoords(20,0) = 0.0; dofCoords(20,1) = 0.0; dofCoords(20,2) = 0.0;
1126 dofCoords(21,0) = 0.0; dofCoords(21,1) = 0.0; dofCoords(21,2) = -1.0;
1127 dofCoords(22,0) = 0.0; dofCoords(22,1) = 0.0; dofCoords(22,2) = 1.0;
1128 dofCoords(23,0) = -1.0; dofCoords(23,1) = 0.0; dofCoords(23,2) = 0.0;
1129 dofCoords(24,0) = 1.0; dofCoords(24,1) = 0.0; dofCoords(24,2) = 0.0;
1130 dofCoords(25,0) = 0.0; dofCoords(25,1) = -1.0; dofCoords(25,2) = 0.0;
1131 dofCoords(26,0) = 0.0; dofCoords(26,1) = 1.0; dofCoords(26,2) = 0.0;
1134 this->
dofCoords_ = Kokkos::create_mirror_view(
typename DT::memory_space(), dofCoords);
1135 Kokkos::deep_copy(this->
dofCoords_, dofCoords);
Basis_HGRAD_HEX_DEG2_FEM()
Constructor.
ECoordinates basisCoordinates_
ordinal_type basisDegree_
void setOrdinalTagData(OrdinalTypeView3D &tagToOrdinal, OrdinalTypeView2D &ordinalToTag, const OrdinalTypeView1D tags, const ordinal_type basisCard, const ordinal_type tagSize, const ordinal_type posScDim, const ordinal_type posScOrd, const ordinal_type posDfOrd)
OrdinalTypeArray2DHost ordinalToTag_
Kokkos::DynRankView< scalarType, DeviceType > dofCoords_
ordinal_type basisCardinality_
OrdinalTypeArray3DHost tagToOrdinal_
shards::CellTopology basisCellTopology_
EFunctionSpace functionSpace_
See Intrepid2::Basis_HGRAD_HEX_DEG2_FEM.