Intrepid2
Intrepid2_HGRAD_HEX_C2_FEMDef.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Intrepid2 Package
4//
5// Copyright 2007 NTESS and the Intrepid2 contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
16#ifndef __INTREPID2_HGRAD_HEX_C2_FEM_DEF_HPP__
17#define __INTREPID2_HGRAD_HEX_C2_FEM_DEF_HPP__
18
19namespace Intrepid2 {
20
21 // -------------------------------------------------------------------------------------
22 namespace Impl {
23
24 template<bool serendipity>
25 template<EOperator opType>
26 template<typename OutputViewType,
27 typename inputViewType>
28 KOKKOS_INLINE_FUNCTION
29 void
30 Basis_HGRAD_HEX_DEG2_FEM<serendipity>::Serial<opType>::
31 getValues( OutputViewType output,
32 const inputViewType input ) {
33 switch (opType) {
34 case OPERATOR_VALUE : {
35 const auto x = input(0);
36 const auto y = input(1);
37 const auto z = input(2);
38
39 if constexpr (!serendipity) {
40 output.access( 0) = 0.125*(-1. + x)*x*(-1. + y)*y*(-1. + z)*z;
41 output.access( 1) = 0.125*x*(1.+ x)*(-1. + y)*y*(-1. + z)*z;
42 output.access( 2) = 0.125*x*(1.+ x)*y*(1.+ y)*(-1. + z)*z;
43 output.access( 3) = 0.125*(-1. + x)*x*y*(1.+ y)*(-1. + z)*z;
44 output.access( 4) = 0.125*(-1. + x)*x*(-1. + y)*y*z*(1.+ z);
45 output.access( 5) = 0.125*x*(1.+ x)*(-1. + y)*y*z*(1.+ z);
46 output.access( 6) = 0.125*x*(1.+ x)*y*(1.+ y)*z*(1.+ z);
47 output.access( 7) = 0.125*(-1. + x)*x*y*(1.+ y)*z*(1.+ z);
48 output.access( 8) = 0.25*(1. - x)*(1. + x)*(-1. + y)*y*(-1. + z)*z;
49 output.access( 9) = 0.25*x*(1.+ x)*(1. - y)*(1. + y)*(-1. + z)*z;
50 output.access(10) = 0.25*(1. - x)*(1. + x)*y*(1.+ y)*(-1. + z)*z;
51 output.access(11) = 0.25*(-1. + x)*x*(1. - y)*(1. + y)*(-1. + z)*z;
52 output.access(12) = 0.25*(-1. + x)*x*(-1. + y)*y*(1. - z)*(1. + z);
53 output.access(13) = 0.25*x*(1.+ x)*(-1. + y)*y*(1. - z)*(1. + z);
54 output.access(14) = 0.25*x*(1.+ x)*y*(1.+ y)*(1. - z)*(1. + z);
55 output.access(15) = 0.25*(-1. + x)*x*y*(1.+ y)*(1. - z)*(1. + z);
56 output.access(16) = 0.25*(1. - x)*(1. + x)*(-1. + y)*y*z*(1.+ z);
57 output.access(17) = 0.25*x*(1.+ x)*(1. - y)*(1. + y)*z*(1.+ z);
58 output.access(18) = 0.25*(1. - x)*(1. + x)*y*(1.+ y)*z*(1.+ z);
59 output.access(19) = 0.25*(-1. + x)*x*(1. - y)*(1. + y)*z*(1.+ z);
60
61 output.access(20) = (1. - x)*(1. + x)*(1. - y)*(1. + y)*(1. - z)*(1. + z);
62 output.access(21) = 0.5*(1. - x)*(1. + x)*(1. - y)*(1. + y)*(-1. + z)*z;
63 output.access(22) = 0.5*(1. - x)*(1. + x)*(1. - y)*(1. + y)*z*(1.+ z);
64 output.access(23) = 0.5*(-1. + x)*x*(1. - y)*(1. + y)*(1. - z)*(1. + z);
65 output.access(24) = 0.5*x*(1.+ x)*(1. - y)*(1. + y)*(1. - z)*(1. + z);
66 output.access(25) = 0.5*(1. - x)*(1. + x)*(-1. + y)*y*(1. - z)*(1. + z);
67 output.access(26) = 0.5*(1. - x)*(1. + x)*y*(1.+ y)*(1. - z)*(1. + z);
68
69 } else { //serendipity
70
71 output.access( 0) = 0.125*(1.0 - x)*(1.0 - y)*(1.0 - z)*(-x - y - z - 2.0);
72 output.access( 1) = 0.125*(1.0 + x)*(1.0 - y)*(1.0 - z)*( x - y - z - 2.0);
73 output.access( 2) = 0.125*(1.0 + x)*(1.0 + y)*(1.0 - z)*( x + y - z - 2.0);
74 output.access( 3) = 0.125*(1.0 - x)*(1.0 + y)*(1.0 - z)*(-x + y - z - 2.0);
75 output.access( 4) = 0.125*(1.0 - x)*(1.0 - y)*(1.0 + z)*(-x - y + z - 2.0);
76 output.access( 5) = 0.125*(1.0 + x)*(1.0 - y)*(1.0 + z)*( x - y + z - 2.0);
77 output.access( 6) = 0.125*(1.0 + x)*(1.0 + y)*(1.0 + z)*( x + y + z - 2.0);
78 output.access( 7) = 0.125*(1.0 - x)*(1.0 + y)*(1.0 + z)*(-x + y + z - 2.0);
79
80 output.access( 8) = 0.25*(1.0 - x*x)*(1.0 - y)*(1.0 - z);
81 output.access( 9) = 0.25*(1.0 + x)*(1.0 - y*y)*(1.0 - z);
82 output.access(10) = 0.25*(1.0 - x*x)*(1.0 + y)*(1.0 - z);
83 output.access(11) = 0.25*(1.0 - x)*(1.0 - y*y)*(1.0 - z);
84
85 output.access(12) = 0.25*(1.0 - x)*(1.0 - y)*(1.0 - z*z);
86 output.access(13) = 0.25*(1.0 + x)*(1.0 - y)*(1.0 - z*z);
87 output.access(14) = 0.25*(1.0 + x)*(1.0 + y)*(1.0 - z*z);
88 output.access(15) = 0.25*(1.0 - x)*(1.0 + y)*(1.0 - z*z);
89
90 output.access(16) = 0.25*(1.0 - x*x)*(1.0 - y)*(1.0 + z);
91 output.access(17) = 0.25*(1.0 + x)*(1.0 - y*y)*(1.0 + z);
92 output.access(18) = 0.25*(1.0 - x*x)*(1.0 + y)*(1.0 + z);
93 output.access(19) = 0.25*(1.0 - x)*(1.0 - y*y)*(1.0 + z);
94 }
95
96 break;
97 }
98 case OPERATOR_GRAD :
99 case OPERATOR_D1 : {
100 const auto x = input(0);
101 const auto y = input(1);
102 const auto z = input(2);
103
104 // output.access is a rank-2 array with dimensions (basisCardinality_, spaceDim)
105 if constexpr (!serendipity) {
106 output.access(0, 0) = (-0.125 + 0.25*x)*(-1. + y)*y*(-1. + z)*z;
107 output.access(0, 1) = (-1. + x)*x*(-0.125 + 0.25*y)*(-1. + z)*z;
108 output.access(0, 2) = (-1. + x)*x*(-1. + y)*y*(-0.125 + 0.25*z);
109
110 output.access(1, 0) = (0.125 + 0.25*x)*(-1. + y)*y*(-1. + z)*z;
111 output.access(1, 1) = x*(1. + x)*(-0.125 + 0.25*y)*(-1. + z)*z;
112 output.access(1, 2) = x*(1. + x)*(-1. + y)*y*(-0.125 + 0.25*z);
113
114 output.access(2, 0) = (0.125 + 0.25*x)*y*(1. + y)*(-1. + z)*z;
115 output.access(2, 1) = x*(1. + x)*(0.125 + 0.25*y)*(-1. + z)*z;
116 output.access(2, 2) = x*(1. + x)*y*(1. + y)*(-0.125 + 0.25*z);
117
118 output.access(3, 0) = (-0.125 + 0.25*x)*y*(1. + y)*(-1. + z)*z;
119 output.access(3, 1) = (-1. + x)*x*(0.125 + 0.25*y)*(-1. + z)*z;
120 output.access(3, 2) = (-1. + x)*x*y*(1. + y)*(-0.125 + 0.25*z);
121
122 output.access(4, 0) = (-0.125 + 0.25*x)*(-1. + y)*y*z*(1. + z);
123 output.access(4, 1) = (-1. + x)*x*(-0.125 + 0.25*y)*z*(1. + z);
124 output.access(4, 2) = (-1. + x)*x*(-1. + y)*y*(0.125 + 0.25*z);
125
126 output.access(5, 0) = (0.125 + 0.25*x)*(-1. + y)*y*z*(1. + z);
127 output.access(5, 1) = x*(1. + x)*(-0.125 + 0.25*y)*z*(1. + z);
128 output.access(5, 2) = x*(1. + x)*(-1. + y)*y*(0.125 + 0.25*z);
129
130 output.access(6, 0) = (0.125 + 0.25*x)*y*(1. + y)*z*(1. + z);
131 output.access(6, 1) = x*(1. + x)*(0.125 + 0.25*y)*z*(1. + z);
132 output.access(6, 2) = x*(1. + x)*y*(1. + y)*(0.125 + 0.25*z);
133
134 output.access(7, 0) = (-0.125 + 0.25*x)*y*(1. + y)*z*(1. + z);
135 output.access(7, 1) = (-1. + x)*x*(0.125 + 0.25*y)*z*(1. + z);
136 output.access(7, 2) = (-1. + x)*x*y*(1. + y)*(0.125 + 0.25*z);
137
138 output.access(8, 0) = -0.5*x*(-1. + y)*y*(-1. + z)*z;
139 output.access(8, 1) = (1. - x)*(1. + x)*(-0.25 + 0.5*y)*(-1. + z)*z;
140 output.access(8, 2) = (1. - x)*(1. + x)*(-1. + y)*y*(-0.25 + 0.5*z);
141
142 output.access(9, 0) = (0.25 + 0.5*x)*(1. - y)*(1. + y)*(-1. + z)*z;
143 output.access(9, 1) = x*(1. + x)*(-0.5*y)*(-1. + z)*z;
144 output.access(9, 2) = x*(1. + x)*(1. - y)*(1. + y)*(-0.25 + 0.5*z);
145
146 output.access(10, 0) = -0.5*x*y*(1. + y)*(-1. + z)*z;
147 output.access(10, 1) = (1. - x)*(1. + x)*(0.25 + 0.5*y)*(-1. + z)*z;
148 output.access(10, 2) = (1. - x)*(1. + x)*y*(1. + y)*(-0.25 + 0.5*z);
149
150 output.access(11, 0) = (-0.25 + 0.5*x)*(1. - y)*(1. + y)*(-1. + z)*z;
151 output.access(11, 1) = (-1. + x)*x*(-0.5*y)*(-1. + z)*z;
152 output.access(11, 2) = (-1. + x)*x*(1. - y)*(1. + y)*(-0.25 + 0.5*z);
153
154 output.access(12, 0) = (-0.25 + 0.5*x)*(-1. + y)*y*(1. - z)*(1. + z);
155 output.access(12, 1) = (-1. + x)*x*(-0.25 + 0.5*y)*(1. - z)*(1. + z);
156 output.access(12, 2) = (-1. + x)*x*(-1. + y)*y*(-0.5*z);
157
158 output.access(13, 0) = (0.25 + 0.5*x)*(-1. + y)*y*(1. - z)*(1. + z);
159 output.access(13, 1) = x*(1. + x)*(-0.25 + 0.5*y)*(1. - z)*(1. + z);
160 output.access(13, 2) = x*(1. + x)*(-1. + y)*y*(-0.5*z);
161
162 output.access(14, 0) = (0.25 + 0.5*x)*y*(1. + y)*(1. - z)*(1. + z);
163 output.access(14, 1) = x*(1. + x)*(0.25 + 0.5*y)*(1. - z)*(1. + z);
164 output.access(14, 2) = x*(1. + x)*y*(1. + y)*(-0.5*z);
165
166 output.access(15, 0) = (-0.25 + 0.5*x)*y*(1. + y)*(1. - z)*(1. + z);
167 output.access(15, 1) = (-1. + x)*x*(0.25 + 0.5*y)*(1. - z)*(1. + z);
168 output.access(15, 2) = (-1. + x)*x*y*(1. + y)*(-0.5*z);
169
170 output.access(16, 0) = -0.5*x*(-1. + y)*y*z*(1. + z);
171 output.access(16, 1) = (1. - x)*(1. + x)*(-0.25 + 0.5*y)*z*(1. + z);
172 output.access(16, 2) = (1. - x)*(1. + x)*(-1. + y)*y*(0.25 + 0.5*z);
173
174 output.access(17, 0) = (0.25 + 0.5*x)*(1. - y)*(1. + y)*z*(1. + z);
175 output.access(17, 1) = x*(1. + x)*(-0.5*y)*z*(1. + z);
176 output.access(17, 2) = x*(1. + x)*(1. - y)*(1. + y)*(0.25 + 0.5*z);
177
178 output.access(18, 0) = -0.5*x*y*(1. + y)*z*(1. + z);
179 output.access(18, 1) = (1. - x)*(1. + x)*(0.25 + 0.5*y)*z*(1. + z);
180 output.access(18, 2) = (1. - x)*(1. + x)*y*(1. + y)*(0.25 + 0.5*z);
181
182 output.access(19, 0) = (-0.25 + 0.5*x)*(1. - y)*(1. + y)*z*(1. + z);
183 output.access(19, 1) = (-1. + x)*x*(-0.5*y)*z*(1. + z);
184 output.access(19, 2) = (-1. + x)*x*(1. - y)*(1. + y)*(0.25 + 0.5*z);
185
186 output.access(20, 0) = -2.*x*(1. - y)*(1. + y)*(1. - z)*(1. + z);
187 output.access(20, 1) = (1. - x)*(1. + x)*(-2.*y)*(1. - z)*(1. + z);
188 output.access(20, 2) = (1. - x)*(1. + x)*(1. - y)*(1. + y)*(-2.*z);
189
190 output.access(21, 0) = -x*(1. - y)*(1. + y)*(-1. + z)*z;
191 output.access(21, 1) = (1. - x)*(1. + x)*(-y)*(-1. + z)*z;
192 output.access(21, 2) = (1. - x)*(1. + x)*(1. - y)*(1. + y)*(-0.5 + z);
193
194 output.access(22, 0) = -x*(1. - y)*(1. + y)*z*(1. + z);
195 output.access(22, 1) = (1. - x)*(1. + x)*(-y)*z*(1. + z);
196 output.access(22, 2) = (1. - x)*(1. + x)*(1. - y)*(1. + y)*(0.5 + z);
197
198 output.access(23, 0) = (-0.5 + x)*(1. - y)*(1. + y)*(1. - z)*(1. + z);
199 output.access(23, 1) = (-1. + x)*x*(-y)*(1. - z)*(1. + z);
200 output.access(23, 2) = (-1. + x)*x*(1. - y)*(1. + y)*(-z);
201
202 output.access(24, 0) = (0.5 + x)*(1. - y)*(1. + y)*(1. - z)*(1. + z);
203 output.access(24, 1) = x*(1. + x)*(-y)*(1. - z)*(1. + z);
204 output.access(24, 2) = x*(1. + x)*(1. - y)*(1. + y)*(-z);
205
206 output.access(25, 0) = -x*(-1. + y)*y*(1. - z)*(1. + z);
207 output.access(25, 1) = (1. - x)*(1. + x)*(-0.5 + y)*(1. - z)*(1. + z);
208 output.access(25, 2) = (1. - x)*(1. + x)*(-1. + y)*y*(-z);
209
210 output.access(26, 0) = -x*y*(1. + y)*(1. - z)*(1. + z);
211 output.access(26, 1) = (1. - x)*(1. + x)*(0.5 + y)*(1. - z)*(1. + z);
212 output.access(26, 2) = (1. - x)*(1. + x)*y*(1. + y)*(-z);
213
214 } else { //serendipity
215
216 output.access(0, 0) = -0.125*(1.0-y)*(1.0-z)*(-x-y-z-2.0) - 0.125*(1.0-x)*(1.0-y)*(1.0-z);
217 output.access(0, 1) = -0.125*(1.0-x)*(1.0-z)*(-x-y-z-2.0) - 0.125*(1.0-x)*(1.0-y)*(1.0-z);
218 output.access(0, 2) = -0.125*(1.0-x)*(1.0-y)*(-x-y-z-2.0) - 0.125*(1.0-x)*(1.0-y)*(1.0-z);
219
220 output.access(1, 0) = 0.125*(1.0-y)*(1.0-z)*( x-y-z-2.0) + 0.125*(1.0+x)*(1.0-y)*(1.0-z);
221 output.access(1, 1) = -0.125*(1.0+x)*(1.0-z)*( x-y-z-2.0) - 0.125*(1.0+x)*(1.0-y)*(1.0-z);
222 output.access(1, 2) = -0.125*(1.0+x)*(1.0-y)*( x-y-z-2.0) - 0.125*(1.0+x)*(1.0-y)*(1.0-z);
223
224 output.access(2, 0) = 0.125*(1.0+y)*(1.0-z)*( x+y-z-2.0) + 0.125*(1.0+x)*(1.0+y)*(1.0-z);
225 output.access(2, 1) = 0.125*(1.0+x)*(1.0-z)*( x+y-z-2.0) + 0.125*(1.0+x)*(1.0+y)*(1.0-z);
226 output.access(2, 2) = -0.125*(1.0+x)*(1.0+y)*( x+y-z-2.0) - 0.125*(1.0+x)*(1.0+y)*(1.0-z);
227
228 output.access(3, 0) = -0.125*(1.0+y)*(1.0-z)*(-x+y-z-2.0) - 0.125*(1.0-x)*(1.0+y)*(1.0-z);
229 output.access(3, 1) = 0.125*(1.0-x)*(1.0-z)*(-x+y-z-2.0) + 0.125*(1.0-x)*(1.0+y)*(1.0-z);
230 output.access(3, 2) = -0.125*(1.0-x)*(1.0+y)*(-x+y-z-2.0) - 0.125*(1.0-x)*(1.0+y)*(1.0-z);
231
232 output.access(4, 0) = -0.125*(1.0-y)*(1.0+z)*(-x-y+z-2.0) - 0.125*(1.0-x)*(1.0-y)*(1.0+z);
233 output.access(4, 1) = -0.125*(1.0-x)*(1.0+z)*(-x-y+z-2.0) - 0.125*(1.0-x)*(1.0-y)*(1.0+z);
234 output.access(4, 2) = 0.125*(1.0-x)*(1.0-y)*(-x-y+z-2.0) + 0.125*(1.0-x)*(1.0-y)*(1.0+z);
235
236 output.access(5, 0) = 0.125*(1.0-y)*(1.0+z)*( x-y+z-2.0) + 0.125*(1.0+x)*(1.0-y)*(1.0+z);
237 output.access(5, 1) = -0.125*(1.0+x)*(1.0+z)*( x-y+z-2.0) - 0.125*(1.0+x)*(1.0-y)*(1.0+z);
238 output.access(5, 2) = 0.125*(1.0+x)*(1.0-y)*( x-y+z-2.0) + 0.125*(1.0+x)*(1.0-y)*(1.0+z);
239
240 output.access(6, 0) = 0.125*(1.0+y)*(1.0+z)*( x+y+z-2.0) + 0.125*(1.0+x)*(1.0+y)*(1.0+z);
241 output.access(6, 1) = 0.125*(1.0+x)*(1.0+z)*( x+y+z-2.0) + 0.125*(1.0+x)*(1.0+y)*(1.0+z);
242 output.access(6, 2) = 0.125*(1.0+x)*(1.0+y)*( x+y+z-2.0) + 0.125*(1.0+x)*(1.0+y)*(1.0+z);
243
244 output.access(7, 0) = -0.125*(1.0+y)*(1.0+z)*(-x+y+z-2.0) - 0.125*(1.0-x)*(1.0+y)*(1.0+z);
245 output.access(7, 1) = 0.125*(1.0-x)*(1.0+z)*(-x+y+z-2.0) + 0.125*(1.0-x)*(1.0+y)*(1.0+z);
246 output.access(7, 2) = 0.125*(1.0-x)*(1.0+y)*(-x+y+z-2.0) + 0.125*(1.0-x)*(1.0+y)*(1.0+z);
247
248 output.access(8, 0) = -0.5*x*(1.0-y)*(1.0-z);
249 output.access(8, 1) = -0.25*(1.0-x*x)*(1.0-z);
250 output.access(8, 2) = -0.25*(1.0-x*x)*(1.0-y);
251
252 output.access(9, 0) = 0.25*(1.0-y*y)*(1.0-z);
253 output.access(9, 1) = -0.5*y*(1.0+x)*(1.0-z);
254 output.access(9, 2) = -0.25*(1.0+x)*(1.0-y*y);
255
256 output.access(10, 0) = -0.5*x*(1.0+y)*(1.0-z);
257 output.access(10, 1) = 0.25*(1.0-x*x)*(1.0-z);
258 output.access(10, 2) = -0.25*(1.0-x*x)*(1.0+y);
259
260 output.access(11, 0) = -0.25*(1.0-y*y)*(1.0-z);
261 output.access(11, 1) = -0.5*y*(1.0-x)*(1.0-z);
262 output.access(11, 2) = -0.25*(1.0-x)*(1.0-y*y);
263
264 output.access(12, 0) = -0.25*(1.0-y)*(1.0-z*z);
265 output.access(12, 1) = -0.25*(1.0-x)*(1.0-z*z);
266 output.access(12, 2) = -0.5*z*(1.0-x)*(1.0-y);
267
268 output.access(13, 0) = 0.25*(1.0-y)*(1.0-z*z);
269 output.access(13, 1) = -0.25*(1.0+x)*(1.0-z*z);
270 output.access(13, 2) = -0.5*z*(1.0+x)*(1.0-y);
271
272 output.access(14, 0) = 0.25*(1.0+y)*(1.0-z*z);
273 output.access(14, 1) = 0.25*(1.0+x)*(1.0-z*z);
274 output.access(14, 2) = -0.5*z*(1.0+x)*(1.0+y);
275
276 output.access(15, 0) = -0.25*(1.0+y)*(1.0-z*z);
277 output.access(15, 1) = 0.25*(1.0-x)*(1.0-z*z);
278 output.access(15, 2) = -0.5*z*(1.0-x)*(1.0+y);
279
280 output.access(16, 0) = -0.5*x*(1.0-y)*(1.0+z);
281 output.access(16, 1) = -0.25*(1.0-x*x)*(1.0+z);
282 output.access(16, 2) = 0.25*(1.0-x*x)*(1.0-y);
283
284 output.access(17, 0) = 0.25*(1.0-y*y)*(1.0+z);
285 output.access(17, 1) = -0.5*y*(1.0+x)*(1.0+z);
286 output.access(17, 2) = 0.25*(1.0+x)*(1.0-y*y);
287
288 output.access(18, 0) = -0.5*x*(1.0+y)*(1.0+z);
289 output.access(18, 1) = 0.25*(1.0-x*x)*(1.0+z);
290 output.access(18, 2) = 0.25*(1.0-x*x)*(1.0+y);
291
292 output.access(19, 0) = -0.25*(1.0-y*y)*(1.0+z);
293 output.access(19, 1) = -0.5*y*(1.0-x)*(1.0+z);
294 output.access(19, 2) = 0.25*(1.0-x)*(1.0-y*y);
295
296 }
297 break;
298 }
299 case OPERATOR_D2 : {
300 const auto x = input(0);
301 const auto y = input(1);
302 const auto z = input(2);
303
304 // output.access is a rank-2 array with dimensions (basisCardinality_, D2Cardinality = 6)
305 if constexpr (!serendipity) {
306 output.access(0, 0) = 0.25*(-1. + y)*y*(-1. + z)*z;
307 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;
308 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);
309 output.access(0, 3) = 0.25*(-1. + x)*x*(-1. + z)*z;
310 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);
311 output.access(0, 5) = 0.25*(-1. + x)*x*(-1. + y)*y;
312
313 output.access(1, 0) = 0.25*(-1. + y)*y*(-1. + z)*z;
314 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;
315 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);
316 output.access(1, 3) = 0.25*x*(1 + x)*(-1. + z)*z;
317 output.access(1, 4) = x*(1. + x)*(0.125 + y*(-0.25 + 0.5*z) - 0.25*z);
318 output.access(1, 5) = 0.25*x*(1 + x)*(-1. + y)*y;
319
320 output.access(2, 0) = 0.25*y*(1 + y)*(-1. + z)*z;
321 output.access(2, 1) = (0.125 + x*(0.25 + 0.5*y) + 0.25*y)*(-1. + z)*z;
322 output.access(2, 2) = y*(1. + y)*(-0.125 + x*(-0.25 + 0.5*z) + 0.25*z);
323 output.access(2, 3) = 0.25*x*(1 + x)*(-1. + z)*z;
324 output.access(2, 4) = x*(1. + x)*(-0.125 + y*(-0.25 + 0.5*z) + 0.25*z);
325 output.access(2, 5) = 0.25*x*(1 + x)*y*(1 + y);
326
327 output.access(3, 0) = 0.25*y*(1 + y)*(-1. + z)*z;
328 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;
329 output.access(3, 2) = y*(1. + y)*(0.125 + x*(-0.25 + 0.5*z) - 0.25*z);
330 output.access(3, 3) = 0.25*(-1. + x)*x*(-1. + z)*z;
331 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);
332 output.access(3, 5) = 0.25*(-1. + x)*x*y*(1 + y);
333
334 output.access(4, 0) = 0.25*(-1. + y)*y*z*(1 + z);
335 output.access(4, 1) = (0.125 + x*(-0.25 + 0.5*y) - 0.25*y)*z*(1. + z);
336 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);
337 output.access(4, 3) = 0.25*(-1. + x)*x*z*(1 + z);
338 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);
339 output.access(4, 5) = 0.25*(-1. + x)*x*(-1. + y)*y;
340
341 output.access(5, 0) = 0.25*(-1. + y)*y*z*(1 + z);
342 output.access(5, 1) = (-0.125 + x*(-0.25 + 0.5*y) + 0.25*y)*z*(1. + z);
343 output.access(5, 2) = (-1. + y)*y*(0.125 + x*(0.25 + 0.5*z) + 0.25*z);
344 output.access(5, 3) = 0.25*x*(1 + x)*z*(1 + z);
345 output.access(5, 4) = x*(1. + x)*(-0.125 + y*(0.25 + 0.5*z) - 0.25*z);
346 output.access(5, 5) = 0.25*x*(1 + x)*(-1. + y)*y;
347
348 output.access(6, 0) = 0.25*y*(1 + y)*z*(1 + z);
349 output.access(6, 1) = (0.125 + x*(0.25 + 0.5*y) + 0.25*y)*z*(1. + z);
350 output.access(6, 2) = y*(1. + y)*(0.125 + x*(0.25 + 0.5*z) + 0.25*z);
351 output.access(6, 3) = 0.25*x*(1 + x)*z*(1 + z);
352 output.access(6, 4) = x*(1. + x)*(0.125 + y*(0.25 + 0.5*z) + 0.25*z);
353 output.access(6, 5) = 0.25*x*(1 + x)*y*(1 + y);
354
355 output.access(7, 0) = 0.25*y*(1 + y)*z*(1 + z);
356 output.access(7, 1) = (-0.125 + x*(0.25 + 0.5*y) - 0.25*y)*z*(1. + z);
357 output.access(7, 2) = y*(1. + y)*(-0.125 + x*(0.25 + 0.5*z) - 0.25*z);
358 output.access(7, 3) = 0.25*(-1. + x)*x*z*(1 + z);
359 output.access(7, 4) = (-1. + x)*x*(0.125 + y*(0.25 + 0.5*z) + 0.25*z);
360 output.access(7, 5) = 0.25*(-1. + x)*x*y*(1 + y);
361
362 output.access(8, 0) = -0.5*(-1. + y)*y*(-1. + z)*z;
363 output.access(8, 1) = (0. + x*(-0.5 + y))*z + (x*(0.5 - y) )*(z*z);
364 output.access(8, 2) = (y*y)*(x*(0.5 - z) ) + y*(x*(-0.5 + z));
365 output.access(8, 3) = 0.5*(1. - x)*(1. + x)*(-1. + z)*z;
366 output.access(8, 4) = 0.25 + (x*x)*(-0.25 + y*(0.5 - z) + 0.5*z) - 0.5*z + y*(-0.5 + z);
367 output.access(8, 5) = 0.5*(1. - x)*(1. + x)*(-1. + y)*y;
368
369 output.access(9, 0) = 0.5*(1. - y)*(1. + y)*(-1. + z)*z;
370 output.access(9, 1) = (0.5*y + x*(y))*z + (x*(-y) - 0.5*y)*(z*z);
371 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);
372 output.access(9, 3) = -0.5*x*(1 + x)*(-1. + z)*z;
373 output.access(9, 4) = x*(y*(0.5 - z) ) + (x*x)*(y*(0.5 - z) );
374 output.access(9, 5) = 0.5*x*(1 + x)*(1. - y)*(1. + y);
375
376 output.access(10, 0) = -0.5*y*(1 + y)*(-1. + z)*z;
377 output.access(10, 1) = (0. + x*(0.5 + y))*z + (x*(-0.5 - y) )*(z*z);
378 output.access(10, 2) = y*(x*(0.5 - z) ) + (y*y)*(x*(0.5 - z) );
379 output.access(10, 3) = 0.5*(1. - x)*(1. + x)*(-1. + z)*z;
380 output.access(10, 4) = -0.25 + (x*x)*(0.25 + y*(0.5 - z) - 0.5*z) + 0.5*z + y*(-0.5 + z);
381 output.access(10, 5) = 0.5*(1. - x)*(1. + x)*y*(1 + y);
382
383 output.access(11, 0) = 0.5*(1. - y)*(1. + y)*(-1. + z)*z;
384 output.access(11, 1) = (-0.5*y + x*(y))*z + (x*(-y) + 0.5*y)*(z*z);
385 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);
386 output.access(11, 3) = -0.5*(-1. + x)*x*(-1. + z)*z;
387 output.access(11, 4) = (x*x)*(y*(0.5 - z) ) + x*(y*(-0.5 + z));
388 output.access(11, 5) = 0.5*(-1. + x)*x*(1. - y)*(1. + y);
389
390 output.access(12, 0) = 0.5*(-1. + y)*y*(1. - z)*(1. + z);
391 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)));
392 output.access(12, 2) = (y*y)*(x*(-z) + 0.5*z) + y*(-0.5*z + x*(z));
393 output.access(12, 3) = 0.5*(-1. + x)*x*(1. - z)*(1. + z);
394 output.access(12, 4) = (x*x)*(y*(-z) + 0.5*z) + x*(-0.5*z + y*(z));
395 output.access(12, 5) = -0.5*(-1. + x)*x*(-1. + y)*y;
396
397 output.access(13, 0) = 0.5*(-1. + y)*y*(1. - z)*(1. + z);
398 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)));
399 output.access(13, 2) = (y*y)*(x*(-z) - 0.5*z) + y*(0.5*z + x*(z));
400 output.access(13, 3) = 0.5*x*(1 + x)*(1. - z)*(1. + z);
401 output.access(13, 4) = x*(y*(-z) + 0.5*z) + (x*x)*(y*(-z) + 0.5*z);
402 output.access(13, 5) = -0.5*x*(1 + x)*(-1. + y)*y;
403
404 output.access(14, 0) = 0.5*y*(1 + y)*(1. - z)*(1. + z);
405 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)));
406 output.access(14, 2) = y*(x*(-z) - 0.5*z) + (y*y)*(x*(-z) - 0.5*z);
407 output.access(14, 3) = 0.5*x*(1 + x)*(1. - z)*(1. + z);
408 output.access(14, 4) = x*(y*(-z) - 0.5*z) + (x*x)*(y*(-z) - 0.5*z);
409 output.access(14, 5) = -0.5*x*(1 + x)*y*(1 + y);
410
411 output.access(15, 0) = 0.5*y*(1 + y)*(1. - z)*(1. + z);
412 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)));
413 output.access(15, 2) = y*(x*(-z) + 0.5*z) + (y*y)*(x*(-z) + 0.5*z);
414 output.access(15, 3) = 0.5*(-1. + x)*x*(1. - z)*(1. + z);
415 output.access(15, 4) = (x*x)*(y*(-z) - 0.5*z) + x*(0.5*z + y*(z));
416 output.access(15, 5) = -0.5*(-1. + x)*x*y*(1 + y);
417
418 output.access(16, 0) = -0.5*(-1. + y)*y*z*(1 + z);
419 output.access(16, 1) = (x*(0.5 - y) )*z + (x*(0.5 - y) )*(z*z);
420 output.access(16, 2) = (y*y)*(x*(-0.5 - z) ) + y*(x*(0.5 + z));
421 output.access(16, 3) = 0.5*(1. - x)*(1. + x)*z*(1 + z);
422 output.access(16, 4) = -0.25 + (x*x)*(0.25 + y*(-0.5 - z) + 0.5*z) - 0.5*z + y*(0.5 + z);
423 output.access(16, 5) = 0.5*(1. - x)*(1. + x)*(-1. + y)*y;
424
425 output.access(17, 0) = 0.5*(1. - y)*(1. + y)*z*(1 + z);
426 output.access(17, 1) = (x*(-y) - 0.5*y)*z + (x*(-y) - 0.5*y)*(z*z);
427 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);
428 output.access(17, 3) = -0.5*x*(1 + x)*z*(1 + z);
429 output.access(17, 4) = x*(y*(-0.5 - z) ) + (x*x)*(y*(-0.5 - z) );
430 output.access(17, 5) = 0.5*x*(1 + x)*(1. - y)*(1. + y);
431
432 output.access(18, 0) = -0.5*y*(1 + y)*z*(1 + z);
433 output.access(18, 1) = (x*(-0.5 - y) )*z + (x*(-0.5 - y) )*(z*z);
434 output.access(18, 2) = y*(x*(-0.5 - z) ) + (y*y)*(x*(-0.5 - z) );
435 output.access(18, 3) = 0.5*(1. - x)*(1. + x)*z*(1 + z);
436 output.access(18, 4) = 0.25 + (x*x)*(-0.25 + y*(-0.5 - z) - 0.5*z) + 0.5*z + y*(0.5 + z);
437 output.access(18, 5) = 0.5*(1. - x)*(1. + x)*y*(1 + y);
438
439 output.access(19, 0) = 0.5*(1. - y)*(1. + y)*z*(1 + z);
440 output.access(19, 1) = (x*(-y) + 0.5*y)*z + (x*(-y) + 0.5*y)*(z*z);
441 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);
442 output.access(19, 3) = -0.5*(-1. + x)*x*z*(1 + z);
443 output.access(19, 4) = (x*x)*(y*(-0.5 - z) ) + x*(y*(0.5 + z));
444 output.access(19, 5) = 0.5*(-1. + x)*x*(1. - y)*(1. + y);
445
446 output.access(20, 0) = -2.*(1. - y)*(1. + y)*(1. - z)*(1. + z);
447 output.access(20, 1) = -4.*x*y*(-1. + z*z);
448 output.access(20, 2) = x*((y*y)*(-4.*z) + 4.*z);
449 output.access(20, 3) = -2.*(1. - x)*(1. + x)*(1. - z)*(1. + z);
450 output.access(20, 4) = (x*x)*(y*(-4.*z) ) + y*(4.*z);
451 output.access(20, 5) = -2.*(1. - x)*(1. + x)*(1. - y)*(1. + y);
452
453 output.access(21, 0) = -(1. - y)*(1. + y)*(-1. + z)*z;
454 output.access(21, 1) = (x*(-2.*y) )*z + (0. + x*(2.*y))*(z*z);
455 output.access(21, 2) = x*(1. - 2.*z + (y*y)*(-1. + 2.*z));
456 output.access(21, 3) = -(1. - x)*(1. + x)*(-1. + z)*z;
457 output.access(21, 4) = y*(1. - 2.*z) + (x*x)*(y*(-1. + 2.*z));
458 output.access(21, 5) = (1. - x)*(1. + x)*(1. - y)*(1. + y);
459
460 output.access(22, 0) = -(1. - y)*(1. + y)*z*(1 + z);
461 output.access(22, 1) = (0. + x*(2.*y))*z + (0. + x*(2.*y))*(z*z);
462 output.access(22, 2) = x*(-1. - 2.*z + (y*y)*(1. + 2.*z));
463 output.access(22, 3) = -(1. - x)*(1. + x)*z*(1 + z);
464 output.access(22, 4) = y*(-1. - 2.*z) + (x*x)*(y*(1. + 2.*z));
465 output.access(22, 5) = (1. - x)*(1. + x)*(1. - y)*(1. + y);
466
467 output.access(23, 0) = (1. - y)*(1. + y)*(1. - z)*(1. + z);
468 output.access(23, 1) = (-1. + 2.*x)*y*(-1. + z*z);
469 output.access(23, 2) = (-1. + 2.*x)*(-1. + y*y)*z;
470 output.access(23, 3) =-(-1. + x)*x*(1. - z)*(1. + z);
471 output.access(23, 4) = 2.*(-1. + x)*x*y*z;
472 output.access(23, 5) =-(-1. + x)*x*(1. - y)*(1. + y);
473
474 output.access(24, 0) = (1. - y)*(1. + y)*(1. - z)*(1. + z);
475 output.access(24, 1) = (1. + 2.*x)*y*(-1. + z*z);
476 output.access(24, 2) = (1. + 2.*x)*(-1. + y*y)*z;
477 output.access(24, 3) = x*(1. + x)*(-1. + z)*(1. + z);
478 output.access(24, 4) = 2.*x*(1. + x)*y*z;
479 output.access(24, 5) = x*(1. + x)*(-1. + y)*(1. + y);
480
481 output.access(25, 0) = -(-1. + y)*y*(1. - z)*(1. + z);
482 output.access(25, 1) = x*(-1. + 2.*y)*(-1. + z*z);
483 output.access(25, 2) = 2.*x*(-1. + y)*y*z;
484 output.access(25, 3) = (1. - x)*(1. + x)*(1. - z)*(1. + z);
485 output.access(25, 4) = (-1. + x*x)*(-1. + 2.*y)*z;
486 output.access(25, 5) =-(1. - x)*(1. + x)*(-1. + y)*y;
487
488 output.access(26, 0) = y*(1. + y)*(-1. + z)*(1. + z);
489 output.access(26, 1) = x*(1. + 2.*y)*(-1. + z*z);
490 output.access(26, 2) = 2.*x*y*(1. + y)*z;
491 output.access(26, 3) = (-1. + x)*(1. + x)*(-1. + z)*(1. + z);
492 output.access(26, 4) = (-1. + x*x)*(1. + 2.*y)*z;
493 output.access(26, 5) = (-1. + x)*(1. + x)*y*(1. + y);
494
495 } else { //serendipity
496 output.access(0, 0) = 0.25*(1.0 - y)*(1.0 - z);
497 output.access(0, 1) = 0.125*(1.0 - z)*(-2.0*x - 2.0*y - z);
498 output.access(0, 2) = 0.125*(1.0 - y)*(-2.0*x - y - 2.0*z);
499 output.access(0, 3) = 0.25*(1.0 - x)*(1.0 - z);
500 output.access(0, 4) = 0.125*(1.0 - x)*(-x - 2.0*y - 2.0*z);
501 output.access(0, 5) = 0.25*(1.0 - x)*(1.0 - y);
502
503 output.access(1, 0) = 0.25*(1.0 - y)*(1.0 - z);
504 output.access(1, 1) = -0.125*(1.0 - z)*(2.0*x - 2.0*y - z);
505 output.access(1, 2) = -0.125*(1.0 - y)*(2.0*x - y - 2.0*z);
506 output.access(1, 3) = 0.25*(1.0 + x)*(1.0 - z);
507 output.access(1, 4) = 0.125*(1.0 + x)*(x - 2.0*y - 2.0*z);
508 output.access(1, 5) = 0.25*(1.0 + x)*(1.0 - y);
509
510 output.access(2, 0) = 0.25*(1.0 + y)*(1.0 - z);
511 output.access(2, 1) = 0.125*(1.0 - z)*(2.0*x + 2.0*y - z);
512 output.access(2, 2) = -0.125*(1.0 + y)*(2.0*x + y - 2.0*z);
513 output.access(2, 3) = 0.25*(1.0 + x)*(1.0 - z);
514 output.access(2, 4) = -0.125*(1.0 + x)*(x + 2.0*y - 2.0*z);
515 output.access(2, 5) = 0.25*(1.0 + x)*(1.0 + y);
516
517 output.access(3, 0) = 0.25*(1.0 + y)*(1.0 - z);
518 output.access(3, 1) = -0.125*(1.0 - z)*(-2.0*x + 2.0*y - z);
519 output.access(3, 2) = 0.125*(1.0 + y)*(-2.0*x + y - 2.0*z);
520 output.access(3, 3) = 0.25*(1.0 - x)*(1.0 - z);
521 output.access(3, 4) = -0.125*(1.0 - x)*(-x + 2.0*y - 2.0*z);
522 output.access(3, 5) = 0.25*(1.0 - x)*(1.0 + y);
523
524 output.access(4, 0) = 0.25*(1.0 - y)*(1.0 + z);
525 output.access(4, 1) = 0.125*(1.0 + z)*(-2.0*x - 2.0*y + z);
526 output.access(4, 2) = -0.125*(1.0 - y)*(-2.0*x - y + 2.0*z);
527 output.access(4, 3) = 0.25*(1.0 - x)*(1.0 + z);
528 output.access(4, 4) = -0.125*(1.0 - x)*(-x - 2.0*y + 2.0*z);
529 output.access(4, 5) = 0.25*(1.0 - x)*(1.0 - y);
530
531 output.access(5, 0) = 0.25*(1.0 - y)*(1.0 + z);
532 output.access(5, 1) = -0.125*(1.0 + z)*(2.0*x - 2.0*y + z);
533 output.access(5, 2) = 0.125*(1.0 - y)*(2.0*x - y + 2.0*z);
534 output.access(5, 3) = 0.25*(1.0 + x)*(1.0 + z);
535 output.access(5, 4) = -0.125*(1.0 + x)*(x - 2.0*y + 2.0*z);
536 output.access(5, 5) = 0.25*(1.0 + x)*(1.0 - y);
537
538 output.access(6, 0) = 0.25*(1.0 + y)*(1.0 + z);
539 output.access(6, 1) = 0.125*(1.0 + z)*(2.0*x + 2.0*y + z);
540 output.access(6, 2) = 0.125*(1.0 + y)*(2.0*x + y + 2.0*z);
541 output.access(6, 3) = 0.25*(1.0 + x)*(1.0 + z);
542 output.access(6, 4) = 0.125*(1.0 + x)*(x + 2.0*y + 2.0*z);
543 output.access(6, 5) = 0.25*(1.0 + x)*(1.0 + y);
544
545 output.access(7, 0) = 0.25*(1.0 + y)*(1.0 + z);
546 output.access(7, 1) = -0.125*(1.0 + z)*(-2.0*x + 2.0*y + z);
547 output.access(7, 2) = -0.125*(1.0 + y)*(-2.0*x + y + 2.0*z);
548 output.access(7, 3) = 0.25*(1.0 - x)*(1.0 + z);
549 output.access(7, 4) = 0.125*(1.0 - x)*(-x + 2.0*y + 2.0*z);
550 output.access(7, 5) = 0.25*(1.0 - x)*(1.0 + y);
551
552 output.access(8, 0) = -0.5*(1.0 - y)*(1.0 - z);
553 output.access(8, 1) = 0.5*x*(1.0 - z);
554 output.access(8, 2) = 0.5*x*(1.0 - y);
555 output.access(8, 3) = 0.0;
556 output.access(8, 4) = 0.25*(1.0 - x*x);
557 output.access(8, 5) = 0.0;
558
559 output.access(9, 0) = 0.0;
560 output.access(9, 1) = -0.5*y*(1.0 - z);
561 output.access(9, 2) = -0.25*(1.0 - y*y);
562 output.access(9, 3) = -0.5*(1.0 + x)*(1.0 - z);
563 output.access(9, 4) = 0.5*y*(1.0 + x);
564 output.access(9, 5) = 0.0;
565
566 output.access(10, 0) = -0.5*(1.0 + y)*(1.0 - z);
567 output.access(10, 1) = -0.5*x*(1.0 - z);
568 output.access(10, 2) = 0.5*x*(1.0 + y);
569 output.access(10, 3) = 0.0;
570 output.access(10, 4) = -0.25*(1.0 - x*x);
571 output.access(10, 5) = 0.0;
572
573 output.access(11, 0) = 0.0;
574 output.access(11, 1) = 0.5*y*(1.0 - z);
575 output.access(11, 2) = 0.25*(1.0 - y*y);
576 output.access(11, 3) = -0.5*(1.0 - x)*(1.0 - z);
577 output.access(11, 4) = 0.5*y*(1.0 - x);
578 output.access(11, 5) = 0.0;
579
580 output.access(12, 0) = 0.0;
581 output.access(12, 1) = 0.25*(1.0 - z*z);
582 output.access(12, 2) = 0.5*z*(1.0 - y);
583 output.access(12, 3) = 0.0;
584 output.access(12, 4) = 0.5*z*(1.0 - x);
585 output.access(12, 5) = -0.5*(1.0 - x)*(1.0 - y);
586
587 output.access(13, 0) = 0.0;
588 output.access(13, 1) = -0.25*(1.0 - z*z);
589 output.access(13, 2) = -0.5*z*(1.0 - y);
590 output.access(13, 3) = 0.0;
591 output.access(13, 4) = 0.5*z*(1.0 + x);
592 output.access(13, 5) = -0.5*(1.0 + x)*(1.0 - y);
593
594 output.access(14, 0) = 0.0;
595 output.access(14, 1) = 0.25*(1.0 - z*z);
596 output.access(14, 2) = -0.5*z*(1.0 + y);
597 output.access(14, 3) = 0.0;
598 output.access(14, 4) = -0.5*z*(1.0 + x);
599 output.access(14, 5) = -0.5*(1.0 + x)*(1.0 + y);
600
601 output.access(15, 0) = 0.0;
602 output.access(15, 1) = -0.25*(1.0 - z*z);
603 output.access(15, 2) = 0.5*z*(1.0 + y);
604 output.access(15, 3) = 0.0;
605 output.access(15, 4) = -0.5*z*(1.0 - x);
606 output.access(14, 5) = -0.5*(1.0 - x)*(1.0 + y);
607
608 output.access(16, 0) = -0.5*(1.0 - y)*(1.0 + z);
609 output.access(16, 1) = 0.5*x*(1.0 + z);
610 output.access(16, 2) = -0.5*x*(1.0 - y);
611 output.access(16, 3) = 0.0;
612 output.access(16, 4) = -0.25*(1.0 - x*x);
613 output.access(16, 5) = 0.0;
614
615 output.access(17, 0) = 0.0;
616 output.access(17, 1) = -0.5*y*(1.0 + z);
617 output.access(17, 2) = 0.25*(1.0 - y*y);
618 output.access(17, 3) = -0.5*(1.0 + x)*(1.0 + z);
619 output.access(17, 4) = -0.5*y*(1.0 + x);
620 output.access(17, 5) = 0.0;
621
622 output.access(18, 0) = -0.5*(1.0 + y)*(1.0 + z);
623 output.access(18, 1) = -0.5*x*(1.0 + z);
624 output.access(18, 2) = -0.5*x*(1.0 + y);
625 output.access(18, 3) = 0.0;
626 output.access(18, 4) = 0.25*(1.0 - x*x);
627 output.access(18, 5) = 0.0;
628
629 output.access(19, 0) = 0.0;
630 output.access(19, 1) = 0.5*y*(1.0 + z);
631 output.access(19, 2) = -0.25*(1.0 - y*y);
632 output.access(19, 3) = -0.5*(1.0 - x)*(1.0 + z);
633 output.access(19, 4) = -0.5*y*(1.0 - x);
634 output.access(19, 5) = 0.0;
635
636 }
637 break;
638 }
639 case OPERATOR_D3 : {
640 const auto x = input(0);
641 const auto y = input(1);
642 const auto z = input(2);
643
644 if constexpr(!serendipity) {
645
646 output.access(0, 0) = 0.;
647 output.access(0, 1) = ((-1.+ 2.*y)*(-1.+ z)*z)/4.;
648 output.access(0, 2) = ((-1.+ y)*y*(-1.+ 2.*z))/4.;
649 output.access(0, 3) = ((-1.+ 2.*x)*(-1.+ z)*z)/4.;
650 output.access(0, 4) = ((-1.+ 2.*x)*(-1.+ 2.*y)*(-1.+ 2.*z))/8.;
651 output.access(0, 5) = ((-1.+ 2.*x)*(-1.+ y)*y)/4.;
652 output.access(0, 6) = 0.;
653 output.access(0, 7) = ((-1.+ x)*x*(-1.+ 2.*z))/4.;
654 output.access(0, 8) = ((-1.+ x)*x*(-1.+ 2.*y))/4.;
655 output.access(0, 9) = 0.;
656
657 output.access(1, 0) = 0.;
658 output.access(1, 1) = ((-1.+ 2.*y)*(-1.+ z)*z)/4.;
659 output.access(1, 2) = ((-1.+ y)*y*(-1.+ 2.*z))/4.;
660 output.access(1, 3) = ((1.+ 2.*x)*(-1.+ z)*z)/4.;
661 output.access(1, 4) = ((1.+ 2.*x)*(-1.+ 2.*y)*(-1.+ 2.*z))/8.;
662 output.access(1, 5) = ((1.+ 2.*x)*(-1.+ y)*y)/4.;
663 output.access(1, 6) = 0.;
664 output.access(1, 7) = (x*(1.+ x)*(-1.+ 2.*z))/4.;
665 output.access(1, 8) = (x*(1.+ x)*(-1.+ 2.*y))/4.;
666 output.access(1, 9) = 0.;
667
668 output.access(2, 0) = 0.;
669 output.access(2, 1) = ((1.+ 2.*y)*(-1.+ z)*z)/4.;
670 output.access(2, 2) = (y*(1.+ y)*(-1.+ 2.*z))/4.;
671 output.access(2, 3) = ((1.+ 2.*x)*(-1.+ z)*z)/4.;
672 output.access(2, 4) = ((1.+ 2.*x)*(1.+ 2.*y)*(-1.+ 2.*z))/8.;
673 output.access(2, 5) = ((1.+ 2.*x)*y*(1.+ y))/4.;
674 output.access(2, 6) = 0.;
675 output.access(2, 7) = (x*(1.+ x)*(-1.+ 2.*z))/4.;
676 output.access(2, 8) = (x*(1.+ x)*(1.+ 2.*y))/4.;
677 output.access(2, 9) = 0.;
678
679 output.access(3, 0) = 0.;
680 output.access(3, 1) = ((1.+ 2.*y)*(-1.+ z)*z)/4.;
681 output.access(3, 2) = (y*(1.+ y)*(-1.+ 2.*z))/4.;
682 output.access(3, 3) = ((-1.+ 2.*x)*(-1.+ z)*z)/4.;
683 output.access(3, 4) = ((-1.+ 2.*x)*(1.+ 2.*y)*(-1.+ 2.*z))/8.;
684 output.access(3, 5) = ((-1.+ 2.*x)*y*(1.+ y))/4.;
685 output.access(3, 6) = 0.;
686 output.access(3, 7) = ((-1.+ x)*x*(-1.+ 2.*z))/4.;
687 output.access(3, 8) = ((-1.+ x)*x*(1.+ 2.*y))/4.;
688 output.access(3, 9) = 0.;
689
690 output.access(4, 0) = 0.;
691 output.access(4, 1) = ((-1.+ 2.*y)*z*(1.+ z))/4.;
692 output.access(4, 2) = ((-1.+ y)*y*(1.+ 2.*z))/4.;
693 output.access(4, 3) = ((-1.+ 2.*x)*z*(1.+ z))/4.;
694 output.access(4, 4) = ((-1.+ 2.*x)*(-1.+ 2.*y)*(1.+ 2.*z))/8.;
695 output.access(4, 5) = ((-1.+ 2.*x)*(-1.+ y)*y)/4.;
696 output.access(4, 6) = 0.;
697 output.access(4, 7) = ((-1.+ x)*x*(1.+ 2.*z))/4.;
698 output.access(4, 8) = ((-1.+ x)*x*(-1.+ 2.*y))/4.;
699 output.access(4, 9) = 0.;
700
701 output.access(5, 0) = 0.;
702 output.access(5, 1) = ((-1.+ 2.*y)*z*(1.+ z))/4.;
703 output.access(5, 2) = ((-1.+ y)*y*(1.+ 2.*z))/4.;
704 output.access(5, 3) = ((1.+ 2.*x)*z*(1.+ z))/4.;
705 output.access(5, 4) = ((1.+ 2.*x)*(-1.+ 2.*y)*(1.+ 2.*z))/8.;
706 output.access(5, 5) = ((1.+ 2.*x)*(-1.+ y)*y)/4.;
707 output.access(5, 6) = 0.;
708 output.access(5, 7) = (x*(1.+ x)*(1.+ 2.*z))/4.;
709 output.access(5, 8) = (x*(1.+ x)*(-1.+ 2.*y))/4.;
710 output.access(5, 9) = 0.;
711
712 output.access(6, 0) = 0.;
713 output.access(6, 1) = ((1.+ 2.*y)*z*(1.+ z))/4.;
714 output.access(6, 2) = (y*(1.+ y)*(1.+ 2.*z))/4.;
715 output.access(6, 3) = ((1.+ 2.*x)*z*(1.+ z))/4.;
716 output.access(6, 4) = ((1.+ 2.*x)*(1.+ 2.*y)*(1.+ 2.*z))/8.;
717 output.access(6, 5) = ((1.+ 2.*x)*y*(1.+ y))/4.;
718 output.access(6, 6) = 0.;
719 output.access(6, 7) = (x*(1.+ x)*(1.+ 2.*z))/4.;
720 output.access(6, 8) = (x*(1.+ x)*(1.+ 2.*y))/4.;
721 output.access(6, 9) = 0.;
722
723 output.access(7, 0) = 0.;
724 output.access(7, 1) = ((1.+ 2.*y)*z*(1.+ z))/4.;
725 output.access(7, 2) = (y*(1.+ y)*(1.+ 2.*z))/4.;
726 output.access(7, 3) = ((-1.+ 2.*x)*z*(1.+ z))/4.;
727 output.access(7, 4) = ((-1.+ 2.*x)*(1.+ 2.*y)*(1.+ 2.*z))/8.;
728 output.access(7, 5) = ((-1.+ 2.*x)*y*(1.+ y))/4.;
729 output.access(7, 6) = 0.;
730 output.access(7, 7) = ((-1.+ x)*x*(1.+ 2.*z))/4.;
731 output.access(7, 8) = ((-1.+ x)*x*(1.+ 2.*y))/4.;
732 output.access(7, 9) = 0.;
733
734 output.access(8, 0) = 0.;
735 output.access(8, 1) = -((-1.+ 2.*y)*(-1.+ z)*z)/2.;
736 output.access(8, 2) = -((-1.+ y)*y*(-1.+ 2.*z))/2.;
737 output.access(8, 3) = -(x*(-1.+ z)*z);
738 output.access(8, 4) = -(x*(-1.+ 2.*y)*(-1.+ 2.*z))/2.;
739 output.access(8, 5) = -(x*(-1.+ y)*y);
740 output.access(8, 6) = 0.;
741 output.access(8, 7) = -((-1.+ (x*x))*(-1.+ 2.*z))/2.;
742 output.access(8, 8) = -((-1.+ (x*x))*(-1.+ 2.*y))/2.;
743 output.access(8, 9) = 0.;
744
745 output.access(9, 0) = 0.;
746 output.access(9, 1) = -(y*(-1.+ z)*z);
747 output.access(9, 2) = -((-1.+ (y*y))*(-1.+ 2.*z))/2.;
748 output.access(9, 3) = -((1.+ 2.*x)*(-1.+ z)*z)/2.;
749 output.access(9, 4) = -((1.+ 2.*x)*y*(-1.+ 2.*z))/2.;
750 output.access(9, 5) = -((1.+ 2.*x)*(-1.+ (y*y)))/2.;
751 output.access(9, 6) = 0.;
752 output.access(9, 7) = -(x*(1.+ x)*(-1.+ 2.*z))/2.;
753 output.access(9, 8) = -(x*(1.+ x)*y);
754 output.access(9, 9) = 0.;
755
756 output.access(10, 0) = 0.;
757 output.access(10, 1) = -((1.+ 2.*y)*(-1.+ z)*z)/2.;
758 output.access(10, 2) = -(y*(1.+ y)*(-1.+ 2.*z))/2.;
759 output.access(10, 3) = -(x*(-1.+ z)*z);
760 output.access(10, 4) = -(x*(1.+ 2.*y)*(-1.+ 2.*z))/2.;
761 output.access(10, 5) = -(x*y*(1.+ y));
762 output.access(10, 6) = 0.;
763 output.access(10, 7) = -((-1.+ (x*x))*(-1.+ 2.*z))/2.;
764 output.access(10, 8) = -((-1.+ (x*x))*(1.+ 2.*y))/2.;
765 output.access(10, 9) = 0.;
766
767 output.access(11, 0) = 0.;
768 output.access(11, 1) = -(y*(-1.+ z)*z);
769 output.access(11, 2) = -((-1.+ (y*y))*(-1.+ 2.*z))/2.;
770 output.access(11, 3) = -((-1.+ 2.*x)*(-1.+ z)*z)/2.;
771 output.access(11, 4) = -((-1.+ 2.*x)*y*(-1.+ 2.*z))/2.;
772 output.access(11, 5) = -((-1.+ 2.*x)*(-1.+ (y*y)))/2.;
773 output.access(11, 6) = 0.;
774 output.access(11, 7) = -((-1.+ x)*x*(-1.+ 2.*z))/2.;
775 output.access(11, 8) = -((-1.+ x)*x*y);
776 output.access(11, 9) = 0.;
777
778 output.access(12, 0) = 0.;
779 output.access(12, 1) = -((-1.+ 2.*y)*(-1.+ (z*z)))/2.;
780 output.access(12, 2) = -((-1.+ y)*y*z);
781 output.access(12, 3) = -((-1.+ 2.*x)*(-1.+ (z*z)))/2.;
782 output.access(12, 4) = -((-1.+ 2.*x)*(-1.+ 2.*y)*z)/2.;
783 output.access(12, 5) = -((-1.+ 2.*x)*(-1.+ y)*y)/2.;
784 output.access(12, 6) = 0.;
785 output.access(12, 7) = -((-1.+ x)*x*z);
786 output.access(12, 8) = -((-1.+ x)*x*(-1.+ 2.*y))/2.;
787 output.access(12, 9) = 0.;
788
789 output.access(13, 0) = 0.;
790 output.access(13, 1) = -((-1.+ 2.*y)*(-1.+ (z*z)))/2.;
791 output.access(13, 2) = -((-1.+ y)*y*z);
792 output.access(13, 3) = -((1.+ 2.*x)*(-1.+ (z*z)))/2.;
793 output.access(13, 4) = -((1.+ 2.*x)*(-1.+ 2.*y)*z)/2.;
794 output.access(13, 5) = -((1.+ 2.*x)*(-1.+ y)*y)/2.;
795 output.access(13, 6) = 0.;
796 output.access(13, 7) = -(x*(1.+ x)*z);
797 output.access(13, 8) = -(x*(1.+ x)*(-1.+ 2.*y))/2.;
798 output.access(13, 9) = 0.;
799
800 output.access(14, 0) = 0.;
801 output.access(14, 1) = -((1.+ 2.*y)*(-1.+ (z*z)))/2.;
802 output.access(14, 2) = -(y*(1.+ y)*z);
803 output.access(14, 3) = -((1.+ 2.*x)*(-1.+ (z*z)))/2.;
804 output.access(14, 4) = -((1.+ 2.*x)*(1.+ 2.*y)*z)/2.;
805 output.access(14, 5) = -((1.+ 2.*x)*y*(1.+ y))/2.;
806 output.access(14, 6) = 0.;
807 output.access(14, 7) = -(x*(1.+ x)*z);
808 output.access(14, 8) = -(x*(1.+ x)*(1.+ 2.*y))/2.;
809 output.access(14, 9) = 0.;
810
811 output.access(15, 0) = 0.;
812 output.access(15, 1) = -((1.+ 2.*y)*(-1.+ (z*z)))/2.;
813 output.access(15, 2) = -(y*(1.+ y)*z);
814 output.access(15, 3) = -((-1.+ 2.*x)*(-1.+ (z*z)))/2.;
815 output.access(15, 4) = -((-1.+ 2.*x)*(1.+ 2.*y)*z)/2.;
816 output.access(15, 5) = -((-1.+ 2.*x)*y*(1.+ y))/2.;
817 output.access(15, 6) = 0.;
818 output.access(15, 7) = -((-1.+ x)*x*z);
819 output.access(15, 8) = -((-1.+ x)*x*(1.+ 2.*y))/2.;
820 output.access(15, 9) = 0.;
821
822 output.access(16, 0) = 0.;
823 output.access(16, 1) = -((-1.+ 2.*y)*z*(1.+ z))/2.;
824 output.access(16, 2) = -((-1.+ y)*y*(1.+ 2.*z))/2.;
825 output.access(16, 3) = -(x*z*(1.+ z));
826 output.access(16, 4) = -(x*(-1.+ 2.*y)*(1.+ 2.*z))/2.;
827 output.access(16, 5) = -(x*(-1.+ y)*y);
828 output.access(16, 6) = 0.;
829 output.access(16, 7) = -((-1.+ (x*x))*(1.+ 2.*z))/2.;
830 output.access(16, 8) = -((-1.+ (x*x))*(-1.+ 2.*y))/2.;
831 output.access(16, 9) = 0.;
832
833 output.access(17, 0) = 0.;
834 output.access(17, 1) = -(y*z*(1.+ z));
835 output.access(17, 2) = -((-1.+ (y*y))*(1.+ 2.*z))/2.;
836 output.access(17, 3) = -((1.+ 2.*x)*z*(1.+ z))/2.;
837 output.access(17, 4) = -((1.+ 2.*x)*y*(1.+ 2.*z))/2.;
838 output.access(17, 5) = -((1.+ 2.*x)*(-1.+ (y*y)))/2.;
839 output.access(17, 6) = 0.;
840 output.access(17, 7) = -(x*(1.+ x)*(1.+ 2.*z))/2.;
841 output.access(17, 8) = -(x*(1.+ x)*y);
842 output.access(17, 9) = 0.;
843
844 output.access(18, 0) = 0.;
845 output.access(18, 1) = -((1.+ 2.*y)*z*(1.+ z))/2.;
846 output.access(18, 2) = -(y*(1.+ y)*(1.+ 2.*z))/2.;
847 output.access(18, 3) = -(x*z*(1.+ z));
848 output.access(18, 4) = -(x*(1.+ 2.*y)*(1.+ 2.*z))/2.;
849 output.access(18, 5) = -(x*y*(1.+ y));
850 output.access(18, 6) = 0.;
851 output.access(18, 7) = -((-1.+ (x*x))*(1.+ 2.*z))/2.;
852 output.access(18, 8) = -((-1.+ (x*x))*(1.+ 2.*y))/2.;
853 output.access(18, 9) = 0.;
854
855 output.access(19, 0) = 0.;
856 output.access(19, 1) = -(y*z*(1.+ z));
857 output.access(19, 2) = -((-1.+ (y*y))*(1.+ 2.*z))/2.;
858 output.access(19, 3) = -((-1.+ 2.*x)*z*(1.+ z))/2.;
859 output.access(19, 4) = -((-1.+ 2.*x)*y*(1.+ 2.*z))/2.;
860 output.access(19, 5) = -((-1.+ 2.*x)*(-1.+ (y*y)))/2.;
861 output.access(19, 6) = 0.;
862 output.access(19, 7) = -((-1.+ x)*x*(1.+ 2.*z))/2.;
863 output.access(19, 8) = -((-1.+ x)*x*y);
864 output.access(19, 9) = 0.;
865
866 output.access(20, 0) = 0.;
867 output.access(20, 1) = -4*y*(-1.+ (z*z));
868 output.access(20, 2) = -4*(-1.+ (y*y))*z;
869 output.access(20, 3) = -4*x*(-1.+ (z*z));
870 output.access(20, 4) = -8*x*y*z;
871 output.access(20, 5) = -4*x*(-1.+ (y*y));
872 output.access(20, 6) = 0.;
873 output.access(20, 7) = -4*(-1.+ (x*x))*z;
874 output.access(20, 8) = -4*(-1.+ (x*x))*y;
875 output.access(20, 9) = 0.;
876
877 output.access(21, 0) = 0.;
878 output.access(21, 1) = 2.*y*(-1.+ z)*z;
879 output.access(21, 2) = (-1.+ (y*y))*(-1.+ 2.*z);
880 output.access(21, 3) = 2.*x*(-1.+ z)*z;
881 output.access(21, 4) = 2.*x*y*(-1.+ 2.*z);
882 output.access(21, 5) = 2.*x*(-1.+ (y*y));
883 output.access(21, 6) = 0.;
884 output.access(21, 7) = (-1.+ (x*x))*(-1.+ 2.*z);
885 output.access(21, 8) = 2.*(-1.+ (x*x))*y;
886 output.access(21, 9) = 0.;
887
888 output.access(22, 0) = 0.;
889 output.access(22, 1) = 2.*y*z*(1.+ z);
890 output.access(22, 2) = (-1.+ (y*y))*(1.+ 2.*z);
891 output.access(22, 3) = 2.*x*z*(1.+ z);
892 output.access(22, 4) = 2.*x*y*(1.+ 2.*z);
893 output.access(22, 5) = 2.*x*(-1.+ (y*y));
894 output.access(22, 6) = 0.;
895 output.access(22, 7) = (-1.+ (x*x))*(1.+ 2.*z);
896 output.access(22, 8) = 2.*(-1.+ (x*x))*y;
897 output.access(22, 9) = 0.;
898
899 output.access(23, 0) = 0.;
900 output.access(23, 1) = 2.*y*(-1.+ (z*z));
901 output.access(23, 2) = 2.*(-1.+ (y*y))*z;
902 output.access(23, 3) = (-1.+ 2.*x)*(-1.+ (z*z));
903 output.access(23, 4) = 2.*(-1.+ 2.*x)*y*z;
904 output.access(23, 5) = (-1.+ 2.*x)*(-1.+ (y*y));
905 output.access(23, 6) = 0.;
906 output.access(23, 7) = 2.*(-1.+ x)*x*z;
907 output.access(23, 8) = 2.*(-1.+ x)*x*y;
908 output.access(23, 9) = 0.;
909
910 output.access(24, 0) = 0.;
911 output.access(24, 1) = 2.*y*(-1.+ (z*z));
912 output.access(24, 2) = 2.*(-1.+ (y*y))*z;
913 output.access(24, 3) = (1.+ 2.*x)*(-1.+ (z*z));
914 output.access(24, 4) = 2.*(1.+ 2.*x)*y*z;
915 output.access(24, 5) = (1.+ 2.*x)*(-1.+ (y*y));
916 output.access(24, 6) = 0.;
917 output.access(24, 7) = 2.*x*(1.+ x)*z;
918 output.access(24, 8) = 2.*x*(1.+ x)*y;
919 output.access(24, 9) = 0.;
920
921 output.access(25, 0) = 0.;
922 output.access(25, 1) = (-1.+ 2.*y)*(-1.+ (z*z));
923 output.access(25, 2) = 2.*(-1.+ y)*y*z;
924 output.access(25, 3) = 2.*x*(-1.+ (z*z));
925 output.access(25, 4) = 2.*x*(-1.+ 2.*y)*z;
926 output.access(25, 5) = 2.*x*(-1.+ y)*y;
927 output.access(25, 6) = 0.;
928 output.access(25, 7) = 2.*(-1.+ (x*x))*z;
929 output.access(25, 8) = (-1.+ (x*x))*(-1.+ 2.*y);
930 output.access(25, 9) = 0.;
931
932 output.access(26, 0) = 0.;
933 output.access(26, 1) = (1.+ 2.*y)*(-1.+ (z*z));
934 output.access(26, 2) = 2.*y*(1.+ y)*z;
935 output.access(26, 3) = 2.*x*(-1.+ (z*z));
936 output.access(26, 4) = 2.*x*(1.+ 2.*y)*z;
937 output.access(26, 5) = 2.*x*y*(1.+ y);
938 output.access(26, 6) = 0.;
939 output.access(26, 7) = 2.*(-1.+ (x*x))*z;
940 output.access(26, 8) = (-1.+ (x*x))*(1.+ 2.*y);
941 output.access(26, 9) = 0.;
942
943 } else { //serendipity
944
945 output.access(0, 0) = 0.0;
946 output.access(0, 1) = -0.25*(1.0 - z);
947 output.access(0, 2) = -0.25*(1.0 - y);
948 output.access(0, 3) = -0.25*(1.0 - z);
949 output.access(0, 4) = -0.125*(-2.0*x - 2.0*y - 2.0*z + 1.0);
950 output.access(0, 5) = -0.25*(1.0 - y);
951 output.access(0, 6) = 0.0;
952 output.access(0, 7) = -0.25*(1.0 - x);
953 output.access(0, 8) = -0.25*(1.0 - x);
954 output.access(0, 9) = 0.0;
955
956 output.access(1, 0) = 0.0;
957 output.access(1, 1) = -0.25*(1.0 - z);
958 output.access(1, 2) = -0.25*(1.0 - y);
959 output.access(1, 3) = 0.25*(1.0 - z);
960 output.access(1, 4) = 0.125*(2.0*x - 2.0*y - 2.0*z + 1.0);
961 output.access(1, 5) = 0.25*(1.0 - y);
962 output.access(1, 6) = 0.0;
963 output.access(1, 7) = -0.25*(1.0 + x);
964 output.access(1, 8) = -0.25*(1.0 + x);
965 output.access(1, 9) = 0.0;
966
967 output.access(2, 0) = 0.0;
968 output.access(2, 1) = 0.25*(1.0 - z);
969 output.access(2, 2) = -0.25*(1.0 + y);
970 output.access(2, 3) = 0.25*(1.0 - z);
971 output.access(2, 4) = -0.125*(2.0*x + 2.0*y - 2.0*z + 1.0);
972 output.access(2, 5) = 0.25*(1.0 + y);
973 output.access(2, 6) = 0.0;
974 output.access(2, 7) = -0.25*(1.0 + x);
975 output.access(2, 8) = 0.25*(1.0 + x);
976 output.access(2, 9) = 0.0;
977
978 output.access(3, 0) = 0.0;
979 output.access(3, 1) = 0.25*(1.0 - z);
980 output.access(3, 2) = -0.25*(1.0 + y);
981 output.access(3, 3) = -0.25*(1.0 - z);
982 output.access(3, 4) = 0.125*(-2.0*x + 2.0*y - 2.0*z + 1.0);
983 output.access(3, 5) = -0.25*(1.0 + y);
984 output.access(3, 6) = 0.0;
985 output.access(3, 7) = -0.25*(1.0 - x);
986 output.access(3, 8) = 0.25*(1.0 - x);
987 output.access(3, 9) = 0.0;
988
989 output.access(4, 0) = 0.0;
990 output.access(4, 1) = -0.25*(1.0 + z);
991 output.access(4, 2) = 0.25*(1.0 - y);
992 output.access(4, 3) = -0.25*(1.0 + z);
993 output.access(4, 4) = 0.125*(-2.0*x - 2.0*y + 2.0*z + 1.0);
994 output.access(4, 5) = -0.25*(1.0 - y);
995 output.access(4, 6) = 0.0;
996 output.access(4, 7) = 0.25*(1.0 - x);
997 output.access(4, 8) = -0.25*(1.0 - x);
998 output.access(4, 9) = 0.0;
999
1000 output.access(5, 0) = 0.0;
1001 output.access(5, 1) = -0.25*(1.0 + z);
1002 output.access(5, 2) = 0.25*(1.0 - y);
1003 output.access(5, 3) = 0.25*(1.0 + z);
1004 output.access(5, 4) = -0.125*(2.0*x - 2.0*y + 2.0*z + 1.0);
1005 output.access(5, 5) = 0.25*(1.0 - y);
1006 output.access(5, 6) = 0.0;
1007 output.access(5, 7) = 0.25*(1.0 + x);
1008 output.access(5, 8) = -0.25*(1.0 + x);
1009 output.access(5, 9) = 0.0;
1010
1011 output.access(6, 0) = 0.0;
1012 output.access(6, 1) = 0.25*(1.0 + z);
1013 output.access(6, 2) = 0.25*(1.0 + y);
1014 output.access(6, 3) = 0.25*(1.0 + z);
1015 output.access(6, 4) = 0.125*(2.0*x + 2.0*y + 2.0*z + 1.0);
1016 output.access(6, 5) = 0.25*(1.0 + y);
1017 output.access(6, 6) = 0.0;
1018 output.access(6, 7) = 0.25*(1.0 + x);
1019 output.access(6, 8) = 0.25*(1.0 + x);
1020 output.access(6, 9) = 0.0;
1021
1022 output.access(7, 0) = 0.0;
1023 output.access(7, 1) = 0.25*(1.0 + z);
1024 output.access(7, 2) = 0.25*(1.0 + y);
1025 output.access(7, 3) = -0.25*(1.0 + z);
1026 output.access(7, 4) = -0.125*(-2.0*x + 2.0*y + 2.0*z + 1.0);
1027 output.access(7, 5) = -0.25*(1.0 + y);
1028 output.access(7, 6) = 0.0;
1029 output.access(7, 7) = 0.25*(1.0 - x);
1030 output.access(7, 8) = 0.25*(1.0 - x);
1031 output.access(7, 9) = 0.0;
1032
1033 output.access(8, 0) = 0.0;
1034 output.access(8, 1) = 0.5*(1.0 - z);
1035 output.access(8, 2) = 0.5*(1.0 - y);
1036 output.access(8, 3) = 0.0;
1037 output.access(8, 4) = -0.5*x;
1038 output.access(8, 5) = 0.0;
1039 output.access(8, 6) = 0.0;
1040 output.access(8, 7) = 0.0;
1041 output.access(8, 8) = 0.0;
1042 output.access(8, 9) = 0.0;
1043
1044 output.access(9, 0) = 0.0;
1045 output.access(9, 1) = 0.0;
1046 output.access(9, 2) = 0.0;
1047 output.access(9, 3) = -0.5*(1.0 - z);
1048 output.access(9, 4) = 0.5*y;
1049 output.access(9, 5) = 0.0;
1050 output.access(9, 6) = 0.0;
1051 output.access(9, 7) = 0.5*(1.0 + x);
1052 output.access(9, 8) = 0.0;
1053 output.access(9, 9) = 0.0;
1054
1055 output.access(10, 0) = 0.0;
1056 output.access(10, 1) = -0.5*(1.0 - z);
1057 output.access(10, 2) = 0.5*(1.0 + y);
1058 output.access(10, 3) = 0.0;
1059 output.access(10, 4) = 0.5*x;
1060 output.access(10, 5) = 0.0;
1061 output.access(10, 6) = 0.0;
1062 output.access(10, 7) = 0.0;
1063 output.access(10, 8) = 0.0;
1064 output.access(10, 9) = 0.0;
1065
1066 output.access(11, 0) = 0.0;
1067 output.access(11, 1) = 0.0;
1068 output.access(11, 2) = 0.0;
1069 output.access(11, 3) = 0.5*(1.0 - z);
1070 output.access(11, 4) = -0.5*y;
1071 output.access(11, 5) = 0.0;
1072 output.access(11, 6) = 0.0;
1073 output.access(11, 7) = 0.5*(1.0 - x);
1074 output.access(11, 8) = 0.0;
1075 output.access(11, 9) = 0.0;
1076
1077 output.access(12, 0) = 0.0;
1078 output.access(12, 1) = 0.0;
1079 output.access(12, 2) = 0.0;
1080 output.access(12, 3) = 0.0;
1081 output.access(12, 4) = -0.5*z;
1082 output.access(12, 5) = 0.5*(1.0 - y);
1083 output.access(12, 6) = 0.0;
1084 output.access(12, 7) = 0.0;
1085 output.access(12, 8) = 0.5*(1.0 - x);
1086 output.access(12, 9) = 0.0;
1087
1088 output.access(13, 0) = 0.0;
1089 output.access(13, 1) = 0.0;
1090 output.access(13, 2) = 0.0;
1091 output.access(13, 3) = 0.0;
1092 output.access(13, 4) = 0.5*z;
1093 output.access(13, 5) = -0.5*(1.0 - y);
1094 output.access(13, 6) = 0.0;
1095 output.access(13, 7) = 0.0;
1096 output.access(13, 8) = 0.5*(1.0 + x);
1097 output.access(13, 9) = 0.0;
1098
1099 output.access(14, 0) = 0.0;
1100 output.access(14, 1) = 0.0;
1101 output.access(14, 2) = 0.0;
1102 output.access(14, 3) = 0.0;
1103 output.access(14, 4) = -0.5*z;
1104 output.access(14, 5) = -0.5*(1.0 + y);
1105 output.access(14, 6) = 0.0;
1106 output.access(14, 7) = 0.0;
1107 output.access(14, 8) = -0.5*(1.0 + x);
1108 output.access(14, 9) = 0.0;
1109
1110 output.access(15, 0) = 0.0;
1111 output.access(15, 1) = 0.0;
1112 output.access(15, 2) = 0.0;
1113 output.access(15, 3) = 0.0;
1114 output.access(15, 4) = 0.5*z;
1115 output.access(15, 5) = 0.5*(1.0 + y);
1116 output.access(15, 6) = 0.0;
1117 output.access(15, 7) = 0.0;
1118 output.access(15, 8) = -0.5*(1.0 - x);
1119 output.access(15, 9) = 0.0;
1120
1121 output.access(16, 0) = 0.0;
1122 output.access(16, 1) = 0.5*(1.0 + z);
1123 output.access(16, 2) = -0.5*(1.0 - y);
1124 output.access(16, 3) = 0.0;
1125 output.access(16, 4) = 0.5*x;
1126 output.access(16, 5) = 0.0;
1127 output.access(16, 6) = 0.0;
1128 output.access(16, 7) = 0.0;
1129 output.access(16, 8) = 0.0;
1130 output.access(16, 9) = 0.0;
1131
1132 output.access(17, 0) = 0.0;
1133 output.access(17, 1) = 0.0;
1134 output.access(17, 2) = 0.0;
1135 output.access(17, 3) = -0.5*(1.0 + z);
1136 output.access(17, 4) = -0.5*y;
1137 output.access(17, 5) = 0.0;
1138 output.access(17, 6) = 0.0;
1139 output.access(17, 7) = -0.5*(1.0 + x);
1140 output.access(17, 8) = 0.0;
1141 output.access(17, 9) = 0.0;
1142
1143 output.access(18, 0) = 0.0;
1144 output.access(18, 1) = -0.5*(1.0 + z);
1145 output.access(18, 2) = -0.5*(1.0 + y);
1146 output.access(18, 3) = 0.0;
1147 output.access(18, 4) = -0.5*x;
1148 output.access(18, 5) = 0.0;
1149 output.access(18, 6) = 0.0;
1150 output.access(18, 7) = 0.0;
1151 output.access(18, 8) = 0.0;
1152 output.access(18, 9) = 0.0;
1153
1154 output.access(19, 0) = 0.0;
1155 output.access(19, 1) = 0.0;
1156 output.access(19, 2) = 0.0;
1157 output.access(19, 3) = 0.5*(1.0 + z);
1158 output.access(19, 4) = 0.5*y;
1159 output.access(19, 5) = 0.0;
1160 output.access(19, 6) = 0.0;
1161 output.access(19, 7) = -0.5*(1.0 - x);
1162 output.access(19, 8) = 0.0;
1163 output.access(19, 9) = 0.0;
1164 }
1165 break;
1166 }
1167 case OPERATOR_D4 : {
1168 // Non-zero entries have Dk (derivative cardinality) indices {3,4,5,7,8,12}, all other entries are 0.
1169 // Intitialize array by zero and then fill only non-zero entries.
1170 const ordinal_type jend = output.extent(1);
1171 const ordinal_type iend = output.extent(0);
1172
1173 for (ordinal_type j=0;j<jend;++j)
1174 for (ordinal_type i=0;i<iend;++i)
1175 output.access(i, j) = 0.0;
1176
1177 if constexpr (!serendipity) {
1178 const auto x = input(0);
1179 const auto y = input(1);
1180 const auto z = input(2);
1181
1182 output.access(0, 3) = ((-1.+ z)*z)/2.;
1183 output.access(0, 4) = ((-1.+ 2.*y)*(-1.+ 2.*z))/4.;
1184 output.access(0, 5) = ((-1.+ y)*y)/2.;
1185 output.access(0, 7) = ((-1.+ 2.*x)*(-1.+ 2.*z))/4.;
1186 output.access(0, 8) = ((-1.+ 2.*x)*(-1.+ 2.*y))/4.;
1187 output.access(0, 12)= ((-1.+ x)*x)/2.;
1188
1189 output.access(1, 3) = ((-1.+ z)*z)/2.;
1190 output.access(1, 4) = ((-1.+ 2.*y)*(-1.+ 2.*z))/4.;
1191 output.access(1, 5) = ((-1.+ y)*y)/2.;
1192 output.access(1, 7) = ((1. + 2.*x)*(-1.+ 2.*z))/4.;
1193 output.access(1, 8) = ((1. + 2.*x)*(-1.+ 2.*y))/4.;
1194 output.access(1, 12)= (x*(1. + x))/2.;
1195
1196 output.access(2, 3) = ((-1.+ z)*z)/2.;
1197 output.access(2, 4) = ((1. + 2.*y)*(-1.+ 2.*z))/4.;
1198 output.access(2, 5) = (y*(1. + y))/2.;
1199 output.access(2, 7) = ((1. + 2.*x)*(-1.+ 2.*z))/4.;
1200 output.access(2, 8) = ((1. + 2.*x)*(1. + 2.*y))/4.;
1201 output.access(2, 12)= (x*(1. + x))/2.;
1202
1203 output.access(3, 3) = ((-1.+ z)*z)/2.;
1204 output.access(3, 4) = ((1. + 2.*y)*(-1.+ 2.*z))/4.;
1205 output.access(3, 5) = (y*(1. + y))/2.;
1206 output.access(3, 7) = ((-1.+ 2.*x)*(-1.+ 2.*z))/4.;
1207 output.access(3, 8) = ((-1.+ 2.*x)*(1. + 2.*y))/4.;
1208 output.access(3, 12)= ((-1.+ x)*x)/2.;
1209
1210 output.access(4, 3) = (z*(1. + z))/2.;
1211 output.access(4, 4) = ((-1.+ 2.*y)*(1. + 2.*z))/4.;
1212 output.access(4, 5) = ((-1.+ y)*y)/2.;
1213 output.access(4, 7) = ((-1.+ 2.*x)*(1. + 2.*z))/4.;
1214 output.access(4, 8) = ((-1.+ 2.*x)*(-1.+ 2.*y))/4.;
1215 output.access(4, 12)= ((-1.+ x)*x)/2.;
1216
1217 output.access(5, 3) = (z*(1. + z))/2.;
1218 output.access(5, 4) = ((-1.+ 2.*y)*(1. + 2.*z))/4.;
1219 output.access(5, 5) = ((-1.+ y)*y)/2.;
1220 output.access(5, 7) = ((1. + 2.*x)*(1. + 2.*z))/4.;
1221 output.access(5, 8) = ((1. + 2.*x)*(-1.+ 2.*y))/4.;
1222 output.access(5, 12)= (x*(1. + x))/2.;
1223
1224 output.access(6, 3) = (z*(1. + z))/2.;
1225 output.access(6, 4) = ((1. + 2.*y)*(1. + 2.*z))/4.;
1226 output.access(6, 5) = (y*(1. + y))/2.;
1227 output.access(6, 7) = ((1. + 2.*x)*(1. + 2.*z))/4.;
1228 output.access(6, 8) = ((1. + 2.*x)*(1. + 2.*y))/4.;
1229 output.access(6, 12)= (x*(1. + x))/2.;
1230
1231 output.access(7, 3) = (z*(1. + z))/2.;
1232 output.access(7, 4) = ((1. + 2.*y)*(1. + 2.*z))/4.;
1233 output.access(7, 5) = (y*(1. + y))/2.;
1234 output.access(7, 7) = ((-1.+ 2.*x)*(1. + 2.*z))/4.;
1235 output.access(7, 8) = ((-1.+ 2.*x)*(1. + 2.*y))/4.;
1236 output.access(7, 12)= ((-1.+ x)*x)/2.;
1237
1238 output.access(8, 3) = -((-1.+ z)*z);
1239 output.access(8, 4) = -0.5 + y + z - 2.*y*z;
1240 output.access(8, 5) = -((-1.+ y)*y);
1241 output.access(8, 7) = x - 2.*x*z;
1242 output.access(8, 8) = x - 2.*x*y;
1243 output.access(8, 12)= 1. - x*x;
1244
1245 output.access(9, 3) = -((-1.+ z)*z);
1246 output.access(9, 4) = y - 2.*y*z;
1247 output.access(9, 5) = 1 - y*y;
1248 output.access(9, 7) = 0.5 + x - z - 2.*x*z;
1249 output.access(9, 8) = -((1. + 2.*x)*y);
1250 output.access(9, 12)= -(x*(1. + x));
1251
1252 output.access(10, 3) = -((-1.+ z)*z);
1253 output.access(10, 4) = 0.5 + y - z - 2.*y*z;
1254 output.access(10, 5) = -(y*(1. + y));
1255 output.access(10, 7) = x - 2.*x*z;
1256 output.access(10, 8) = -(x*(1. + 2.*y));
1257 output.access(10, 12)= 1. - x*x;
1258
1259 output.access(11, 3) = -((-1.+ z)*z);
1260 output.access(11, 4) = y - 2.*y*z;
1261 output.access(11, 5) = 1. - y*y;
1262 output.access(11, 7) = -0.5 + x + z - 2.*x*z;
1263 output.access(11, 8) = y - 2.*x*y;
1264 output.access(11, 12)= -((-1.+ x)*x);
1265
1266 output.access(12, 3) = 1. - z*z;
1267 output.access(12, 4) = z - 2.*y*z;
1268 output.access(12, 5) = -((-1.+ y)*y);
1269 output.access(12, 7) = z - 2.*x*z;
1270 output.access(12, 8) = -0.5 + x + y - 2.*x*y;
1271 output.access(12, 12)= -((-1.+ x)*x);
1272
1273 output.access(13, 3) = 1. - z*z;
1274 output.access(13, 4) = z - 2.*y*z;
1275 output.access(13, 5) = -((-1.+ y)*y);
1276 output.access(13, 7) = -((1. + 2.*x)*z);
1277 output.access(13, 8) = 0.5 + x - y - 2.*x*y;
1278 output.access(13, 12)= -(x*(1. + x));
1279
1280 output.access(14, 3) = 1. - z*z;
1281 output.access(14, 4) = -((1. + 2.*y)*z);
1282 output.access(14, 5) = -(y*(1. + y));
1283 output.access(14, 7) = -((1. + 2.*x)*z);
1284 output.access(14, 8) = -((1. + 2.*x)*(1. + 2.*y))/2.;
1285 output.access(14, 12)= -(x*(1. + x));
1286
1287 output.access(15, 3) = 1. - z*z;
1288 output.access(15, 4) = -((1. + 2.*y)*z);
1289 output.access(15, 5) = -(y*(1. + y));
1290 output.access(15, 7) = z - 2.*x*z;
1291 output.access(15, 8) = 0.5 + y - x*(1. + 2.*y);
1292 output.access(15, 12)= -((-1.+ x)*x);
1293
1294 output.access(16, 3) = -(z*(1. + z));
1295 output.access(16, 4) = 0.5 + z - y*(1. + 2.*z);
1296 output.access(16, 5) = -((-1.+ y)*y);
1297 output.access(16, 7) = -(x*(1. + 2.*z));
1298 output.access(16, 8) = x - 2.*x*y;
1299 output.access(16, 12)= 1. - x*x;
1300
1301 output.access(17, 3) = -(z*(1. + z));
1302 output.access(17, 4) = -(y*(1. + 2.*z));
1303 output.access(17, 5) = 1. - y*y;
1304 output.access(17, 7) = -((1. + 2.*x)*(1. + 2.*z))/2.;
1305 output.access(17, 8) = -((1. + 2.*x)*y);
1306 output.access(17, 12)= -(x*(1. + x));
1307
1308 output.access(18, 3) = -(z*(1. + z));
1309 output.access(18, 4) = -((1. + 2.*y)*(1. + 2.*z))/2.;
1310 output.access(18, 5) = -(y*(1. + y));
1311 output.access(18, 7) = -(x*(1. + 2.*z));
1312 output.access(18, 8) = -(x*(1. + 2.*y));
1313 output.access(18, 12)= 1. - x*x;
1314
1315 output.access(19, 3) = -(z*(1. + z));
1316 output.access(19, 4) = -(y*(1. + 2.*z));
1317 output.access(19, 5) = 1. - y*y;
1318 output.access(19, 7) = 0.5 + z - x*(1. + 2.*z);
1319 output.access(19, 8) = y - 2.*x*y;
1320 output.access(19, 12)= -((-1.+ x)*x);
1321
1322 output.access(20, 3) = 4. - 4.*z*z;
1323 output.access(20, 4) = -8.*y*z;
1324 output.access(20, 5) = 4. - 4.*y*y;
1325 output.access(20, 7) = -8.*x*z;
1326 output.access(20, 8) = -8.*x*y;
1327 output.access(20, 12)= 4. - 4.*x*x;
1328
1329 output.access(21, 3) = 2.*(-1.+ z)*z;
1330 output.access(21, 4) = 2.*y*(-1.+ 2.*z);
1331 output.access(21, 5) = 2.*(-1.+ y*y);
1332 output.access(21, 7) = 2.*x*(-1.+ 2.*z);
1333 output.access(21, 8) = 4.*x*y;
1334 output.access(21, 12)= 2.*(-1.+ x*x);
1335
1336 output.access(22, 3) = 2.*z*(1. + z);
1337 output.access(22, 4) = 2.*(y + 2.*y*z);
1338 output.access(22, 5) = 2.*(-1.+ y*y);
1339 output.access(22, 7) = 2.*(x + 2.*x*z);
1340 output.access(22, 8) = 4.*x*y;
1341 output.access(22, 12)= 2.*(-1.+ x*x);
1342
1343 output.access(23, 3) = 2.*(-1.+ z*z);
1344 output.access(23, 4) = 4.*y*z;
1345 output.access(23, 5) = 2.*(-1.+ y*y);
1346 output.access(23, 7) = 2.*(-1.+ 2.*x)*z;
1347 output.access(23, 8) = 2.*(-1.+ 2.*x)*y;
1348 output.access(23, 12)= 2.*(-1.+ x)*x;
1349
1350 output.access(24, 3) = 2.*(-1.+ z*z);
1351 output.access(24, 4) = 4.*y*z;
1352 output.access(24, 5) = 2.*(-1.+ y*y);
1353 output.access(24, 7) = 2.*(z + 2.*x*z);
1354 output.access(24, 8) = 2.*(y + 2.*x*y);
1355 output.access(24, 12)= 2.*x*(1. + x);
1356
1357 output.access(25, 3) = 2.*(-1.+ z*z);
1358 output.access(25, 4) = 2.*(-1.+ 2.*y)*z;
1359 output.access(25, 5) = 2.*(-1.+ y)*y;
1360 output.access(25, 7) = 4.*x*z;
1361 output.access(25, 8) = 2.*x*(-1.+ 2.*y);
1362 output.access(25, 12)= 2.*(-1.+ x*x);
1363
1364 output.access(26, 3) = 2.*(-1.+ z*z);
1365 output.access(26, 4) = 2.*(z + 2.*y*z);
1366 output.access(26, 5) = 2.*y*(1. + y);
1367 output.access(26, 7) = 4.*x*z;
1368 output.access(26, 8) = 2.*(x + 2.*x*y);
1369 output.access(26, 12)= 2.*(-1.+ x*x);
1370
1371 } else { //serendipity
1372
1373 output.access( 0, 4) = 0.25;
1374 output.access( 0, 7) = 0.25;
1375 output.access( 0, 8) = 0.25;
1376
1377 output.access( 1, 4) = 0.25;
1378 output.access( 1, 7) = -0.25;
1379 output.access( 1, 8) = -0.25;
1380
1381 output.access( 2, 4) = -0.25;
1382 output.access( 2, 7) = -0.25;
1383 output.access( 2, 8) = 0.25;
1384
1385 output.access( 3, 4) = -0.25;
1386 output.access( 3, 7) = 0.25;
1387 output.access( 3, 8) = -0.25;
1388
1389 output.access( 4, 4) = -0.25;
1390 output.access( 4, 7) = -0.25;
1391 output.access( 4, 8) = 0.25;
1392
1393 output.access( 5, 4) = -0.25;
1394 output.access( 5, 7) = 0.25;
1395 output.access( 5, 8) = -0.25;
1396
1397 output.access( 6, 4) = 0.25;
1398 output.access( 6, 7) = 0.25;
1399 output.access( 6, 8) = 0.25;
1400
1401 output.access( 7, 4) = 0.25;
1402 output.access( 7, 7) = -0.25;
1403 output.access( 7, 8) = -0.25;
1404
1405 output.access( 8, 4) = -0.5;
1406 output.access( 9, 7) = 0.5;
1407 output.access(10, 4) = 0.5;
1408 output.access(11, 7) = -0.5;
1409 output.access(12, 8) = -0.5;
1410 output.access(13, 8) = 0.5;
1411 output.access(14, 8) = -0.5;
1412 output.access(15, 8) = 0.5;
1413 output.access(16, 4) = 0.5;
1414 output.access(17, 7) = -0.5;
1415 output.access(18, 4) = -0.5;
1416 output.access(19, 7) = 0.5;
1417 }
1418
1419 break;
1420 }
1421 case OPERATOR_MAX : {
1422 const ordinal_type jend = output.extent(1);
1423 const ordinal_type iend = output.extent(0);
1424
1425 for (ordinal_type j=0;j<jend;++j)
1426 for (ordinal_type i=0;i<iend;++i)
1427 output.access(i, j) = 0.0;
1428 break;
1429 }
1430 default: {
1431 INTREPID2_TEST_FOR_ABORT( opType != OPERATOR_VALUE &&
1432 opType != OPERATOR_GRAD &&
1433 opType != OPERATOR_CURL &&
1434 opType != OPERATOR_D1 &&
1435 opType != OPERATOR_D2 &&
1436 opType != OPERATOR_D3 &&
1437 opType != OPERATOR_D4 &&
1438 opType != OPERATOR_MAX,
1439 ">>> ERROR: (Intrepid2::Basis_HGRAD_HEX_DEG2_FEM::Serial::getValues) operator is not supported");
1440
1441 }
1442 }
1443 }
1444
1445 template<bool serendipity>
1446 template<typename DT,
1447 typename outputValueValueType, class ...outputValueProperties,
1448 typename inputPointValueType, class ...inputPointProperties>
1449 void
1450 Basis_HGRAD_HEX_DEG2_FEM<serendipity>::
1451 getValues( const typename DT::execution_space& space,
1452 Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
1453 const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
1454 const EOperator operatorType ) {
1455 typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
1456 typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
1457 typedef typename ExecSpace<typename inputPointViewType::execution_space,typename DT::execution_space>::ExecSpaceType ExecSpaceType;
1458
1459 // Number of evaluation points = dim 0 of inputPoints
1460 const auto loopSize = inputPoints.extent(0);
1461 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(space, 0, loopSize);
1462
1463 switch (operatorType) {
1464
1465 case OPERATOR_VALUE: {
1466 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_VALUE> FunctorType;
1467 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
1468 break;
1469 }
1470 case OPERATOR_GRAD:
1471 case OPERATOR_D1: {
1472 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_GRAD> FunctorType;
1473 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
1474 break;
1475 }
1476 case OPERATOR_CURL:
1477 INTREPID2_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_CURL), std::invalid_argument,
1478 ">>> ERROR (Basis_HGRAD_HEX_DEG2_FEM): CURL is invalid operator for rank-0 (scalar) functions in 3D");
1479 break;
1480
1481 case OPERATOR_DIV:
1482 INTREPID2_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
1483 ">>> ERROR (Basis_HGRAD_HEX_DEG2_FEM): DIV is invalid operator for rank-0 (scalar) functions in 3D");
1484 break;
1485
1486 case OPERATOR_D2: {
1487 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_D2> FunctorType;
1488 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
1489 break;
1490 }
1491 case OPERATOR_D3: {
1492 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_D3> FunctorType;
1493 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
1494 break;
1495 }
1496 case OPERATOR_D4: {
1497 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_D4> FunctorType;
1498 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
1499 break;
1500 }
1501 case OPERATOR_D5:
1502 case OPERATOR_D6: {
1503 INTREPID2_TEST_FOR_EXCEPTION( true, std::invalid_argument,
1504 ">>> ERROR (Basis_HGRAD_HEX_DEG2_FEM): operator not supported");
1505 // break; commented out becuase this always throws
1506 }
1507 case OPERATOR_D7:
1508 case OPERATOR_D8:
1509 case OPERATOR_D9:
1510 case OPERATOR_D10: {
1511 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_MAX> FunctorType;
1512 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
1513 break;
1514 }
1515 default: {
1516 INTREPID2_TEST_FOR_EXCEPTION( !( Intrepid2::isValidOperator(operatorType) ), std::invalid_argument,
1517 ">>> ERROR (Basis_HGRAD_HEX_DEG2_FEM): Invalid operator type");
1518 }
1519 }
1520 }
1521
1522 }
1523
1524 // -------------------------------------------------------------------------------------
1525
1526 template<bool serendipity, typename DT, typename OT, typename PT>
1527 Basis_HGRAD_HEX_DEG2_FEM<serendipity, DT,OT,PT>::
1528 Basis_HGRAD_HEX_DEG2_FEM() {
1529 const ordinal_type spaceDim = 3;
1530 this -> basisCardinality_ = serendipity ? 20 : 27;
1531 this -> basisDegree_ = 2;
1532 this -> basisCellTopologyKey_ = shards::Hexahedron<8>::key;
1533 this -> basisType_ = BASIS_FEM_DEFAULT;
1534 this -> basisCoordinates_ = COORDINATES_CARTESIAN;
1535 this->functionSpace_ = FUNCTION_SPACE_HGRAD;
1536
1537 // initialize tags
1538 {
1539 // Basis-dependent intializations
1540 const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
1541 const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
1542 const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
1543 const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
1544
1545 // An array with local DoF tags assigned to basis functions, in the order of their local enumeration
1546 ordinal_type tags[108] = { 0, 0, 0, 1, // Nodes 0 to 7 follow vertex order of the topology
1547 0, 1, 0, 1,
1548 0, 2, 0, 1,
1549 0, 3, 0, 1,
1550 0, 4, 0, 1,
1551 0, 5, 0, 1,
1552 0, 6, 0, 1,
1553 0, 7, 0, 1,
1554 1, 0, 0, 1, // Node 8 -> edge 0
1555 1, 1, 0, 1, // Node 9 -> edge 1
1556 1, 2, 0, 1, // Node 10 -> edge 2
1557 1, 3, 0, 1, // Node 11 -> edge 3
1558 1, 8, 0, 1, // Node 12 -> edge 8
1559 1, 9, 0, 1, // Node 13 -> edge 9
1560 1,10, 0, 1, // Node 14 -> edge 10
1561 1,11, 0, 1, // Node 15 -> edge 11
1562 1, 4, 0, 1, // Node 16 -> edge 4
1563 1, 5, 0, 1, // Node 17 -> edge 5
1564 1, 6, 0, 1, // Node 18 -> edge 6
1565 1, 7, 0, 1, // Node 19 -> edge 7
1566
1567 // following entries not used for serendipity elements
1568 3, 0, 0, 1, // Node 20 -> Hexahedron
1569 2, 4, 0, 1, // Node 21 -> face 4
1570 2, 5, 0, 1, // Node 22 -> face 5
1571 2, 3, 0, 1, // Node 23 -> face 3
1572 2, 1, 0, 1, // Node 24 -> face 1
1573 2, 0, 0, 1, // Node 25 -> face 0
1574 2, 2, 0, 1, // Node 26 -> face 2
1575 };
1576
1577 // host tags
1578 OrdinalTypeArray1DHost tagView(&tags[0], serendipity ? 80 : 108);
1579
1580 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
1581 this->setOrdinalTagData(this -> tagToOrdinal_,
1582 this -> ordinalToTag_,
1583 tagView,
1584 this -> basisCardinality_,
1585 tagSize,
1586 posScDim,
1587 posScOrd,
1588 posDfOrd);
1589 }
1590 // dofCoords on host and create its mirror view to device
1591 Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
1592 dofCoords("dofCoordsHost", this->basisCardinality_, spaceDim);
1593
1594 dofCoords(0,0) = -1.0; dofCoords(0,1) = -1.0; dofCoords(0,2) = -1.0;
1595 dofCoords(1,0) = 1.0; dofCoords(1,1) = -1.0; dofCoords(1,2) = -1.0;
1596 dofCoords(2,0) = 1.0; dofCoords(2,1) = 1.0; dofCoords(2,2) = -1.0;
1597 dofCoords(3,0) = -1.0; dofCoords(3,1) = 1.0; dofCoords(3,2) = -1.0;
1598 dofCoords(4,0) = -1.0; dofCoords(4,1) = -1.0; dofCoords(4,2) = 1.0;
1599 dofCoords(5,0) = 1.0; dofCoords(5,1) = -1.0; dofCoords(5,2) = 1.0;
1600 dofCoords(6,0) = 1.0; dofCoords(6,1) = 1.0; dofCoords(6,2) = 1.0;
1601 dofCoords(7,0) = -1.0; dofCoords(7,1) = 1.0; dofCoords(7,2) = 1.0;
1602
1603 dofCoords(8,0) = 0.0; dofCoords(8,1) = -1.0; dofCoords(8,2) = -1.0;
1604 dofCoords(9,0) = 1.0; dofCoords(9,1) = 0.0; dofCoords(9,2) = -1.0;
1605 dofCoords(10,0) = 0.0; dofCoords(10,1) = 1.0; dofCoords(10,2) = -1.0;
1606 dofCoords(11,0) = -1.0; dofCoords(11,1) = 0.0; dofCoords(11,2) = -1.0;
1607 dofCoords(12,0) = -1.0; dofCoords(12,1) = -1.0; dofCoords(12,2) = 0.0;
1608 dofCoords(13,0) = 1.0; dofCoords(13,1) = -1.0; dofCoords(13,2) = 0.0;
1609 dofCoords(14,0) = 1.0; dofCoords(14,1) = 1.0; dofCoords(14,2) = 0.0;
1610 dofCoords(15,0) = -1.0; dofCoords(15,1) = 1.0; dofCoords(15,2) = 0.0;
1611 dofCoords(16,0) = 0.0; dofCoords(16,1) = -1.0; dofCoords(16,2) = 1.0;
1612 dofCoords(17,0) = 1.0; dofCoords(17,1) = 0.0; dofCoords(17,2) = 1.0;
1613 dofCoords(18,0) = 0.0; dofCoords(18,1) = 1.0; dofCoords(18,2) = 1.0;
1614 dofCoords(19,0) = -1.0; dofCoords(19,1) = 0.0; dofCoords(19,2) = 1.0;
1615
1616 if constexpr (!serendipity) {
1617 dofCoords(20,0) = 0.0; dofCoords(20,1) = 0.0; dofCoords(20,2) = 0.0;
1618
1619 dofCoords(21,0) = 0.0; dofCoords(21,1) = 0.0; dofCoords(21,2) = -1.0;
1620 dofCoords(22,0) = 0.0; dofCoords(22,1) = 0.0; dofCoords(22,2) = 1.0;
1621 dofCoords(23,0) = -1.0; dofCoords(23,1) = 0.0; dofCoords(23,2) = 0.0;
1622 dofCoords(24,0) = 1.0; dofCoords(24,1) = 0.0; dofCoords(24,2) = 0.0;
1623 dofCoords(25,0) = 0.0; dofCoords(25,1) = -1.0; dofCoords(25,2) = 0.0;
1624 dofCoords(26,0) = 0.0; dofCoords(26,1) = 1.0; dofCoords(26,2) = 0.0;
1625 }
1626
1627 this->dofCoords_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoords);
1628 Kokkos::deep_copy(this->dofCoords_, dofCoords);
1629 }
1630
1631 template<bool serendipity, typename DT, typename OT, typename PT>
1632 void
1633 Basis_HGRAD_HEX_DEG2_FEM<serendipity,DT,OT,PT>::getScratchSpaceSize(
1634 ordinal_type& perTeamSpaceSize,
1635 ordinal_type& perThreadSpaceSize,
1636 const PointViewType inputPoints,
1637 const EOperator operatorType) const {
1638 perTeamSpaceSize = 0;
1639 perThreadSpaceSize = 0;
1640 }
1641
1642 template<bool serendipity, typename DT, typename OT, typename PT>
1643 KOKKOS_INLINE_FUNCTION
1644 void
1645 Basis_HGRAD_HEX_DEG2_FEM<serendipity,DT,OT,PT>::getValues(
1646 OutputViewType outputValues,
1647 const PointViewType inputPoints,
1648 const EOperator operatorType,
1649 const typename Kokkos::TeamPolicy<typename DT::execution_space>::member_type& team_member,
1650 const typename DT::execution_space::scratch_memory_space & scratchStorage,
1651 const ordinal_type subcellDim,
1652 const ordinal_type subcellOrdinal) const {
1653
1654 INTREPID2_TEST_FOR_ABORT( !((subcellDim <= 0) && (subcellOrdinal == -1)),
1655 ">>> ERROR: (Intrepid2::Basis_HGRAD_HEX_DEG2_FEM::getValues), The capability of selecting subsets of basis functions has not been implemented yet.");
1656
1657 (void) scratchStorage; //avoid unused variable warning
1658
1659 const int numPoints = inputPoints.extent(0);
1660
1661 switch(operatorType) {
1662 case OPERATOR_VALUE:
1663 Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=] (ordinal_type& pt) {
1664 auto output = Kokkos::subview( outputValues, Kokkos::ALL(), pt, Kokkos::ALL() );
1665 const auto input = Kokkos::subview( inputPoints, pt, Kokkos::ALL() );
1666 using SerialValue = typename Impl::Basis_HGRAD_HEX_DEG2_FEM<serendipity>::template Serial<OPERATOR_VALUE>;
1667 SerialValue::getValues( output, input);
1668 });
1669 break;
1670 case OPERATOR_GRAD:
1671 Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=] (ordinal_type& pt) {
1672 auto output = Kokkos::subview( outputValues, Kokkos::ALL(), pt, Kokkos::ALL() );
1673 const auto input = Kokkos::subview( inputPoints, pt, Kokkos::ALL() );
1674 using SerialGrad = typename Impl::Basis_HGRAD_HEX_DEG2_FEM<serendipity>::template Serial<OPERATOR_GRAD>;
1675 SerialGrad::getValues( output, input);
1676 });
1677 break;
1678 default: {}
1679 }
1680 }
1681}// namespace Intrepid2
1682#endif
Intrepid2::Basis_HGRAD_HEX_DEG2_FEM::getScratchSpaceSize
virtual void getScratchSpaceSize(ordinal_type &perTeamSpaceSize, ordinal_type &perThreadSpaceSize, const PointViewType inputPointsconst, const EOperator operatorType=OPERATOR_VALUE) const override
Return the size of the scratch space, in bytes, needed for using the team-level implementation of get...
Definition Intrepid2_HGRAD_HEX_C2_FEMDef.hpp:1633
Intrepid2::Basis_HGRAD_HEX_DEG2_FEM::Basis_HGRAD_HEX_DEG2_FEM
Basis_HGRAD_HEX_DEG2_FEM()
Constructor.
Definition Intrepid2_HGRAD_HEX_C2_FEMDef.hpp:1528
Intrepid2::Basis_HGRAD_HEX_DEG2_FEM::getValues
virtual void getValues(const ExecutionSpace &space, OutputViewType outputValues, const PointViewType inputPoints, const EOperator operatorType=OPERATOR_VALUE) const override
Evaluation of a FEM basis on a reference cell.
Definition Intrepid2_HGRAD_HEX_C2_FEM.hpp:222
Intrepid2::Basis< DeviceType, outputValueType, pointValueType >::OrdinalTypeArray1DHost
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
Definition Intrepid2_Basis.hpp:125
Intrepid2::Impl::Basis_HGRAD_HEX_DEG2_FEM
See Intrepid2::Basis_HGRAD_HEX_DEG2_FEM.
Definition Intrepid2_HGRAD_HEX_C2_FEM.hpp:120
Generated by doxygen 1.9.8