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 
19 namespace 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::OrdinalTypeArray1DHost
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
Definition: Intrepid2_Basis.hpp:125
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::Impl::Basis_HGRAD_HEX_DEG2_FEM
See Intrepid2::Basis_HGRAD_HEX_DEG2_FEM.
Definition: Intrepid2_HGRAD_HEX_C2_FEM.hpp:120
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::PointViewType
Kokkos::DynRankView< PointValueType, Kokkos::LayoutStride, DeviceType > PointViewType
View type for input points.
Definition: Intrepid2_Basis.hpp:338
Generated by   doxygen 1.8.5