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