Intrepid
Intrepid_HGRAD_HEX_C2_FEMDef.hpp
Go to the documentation of this file.
1 #ifndef INTREPID_HGRAD_HEX_C2_FEMDEF_HPP
2 #define INTREPID_HGRAD_HEX_C2_FEMDEF_HPP
3 // @HEADER
4 // ************************************************************************
5 //
6 // Intrepid Package
7 // Copyright (2007) Sandia Corporation
8 //
9 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10 // license for use of this work by or on behalf of the U.S. Government.
11 //
12 // Redistribution and use in source and binary forms, with or without
13 // modification, are permitted provided that the following conditions are
14 // met:
15 //
16 // 1. Redistributions of source code must retain the above copyright
17 // notice, this list of conditions and the following disclaimer.
18 //
19 // 2. Redistributions in binary form must reproduce the above copyright
20 // notice, this list of conditions and the following disclaimer in the
21 // documentation and/or other materials provided with the distribution.
22 //
23 // 3. Neither the name of the Corporation nor the names of the
24 // contributors may be used to endorse or promote products derived from
25 // this software without specific prior written permission.
26 //
27 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38 //
39 // Questions? Contact Pavel Bochev (pbboche@sandia.gov)
40 // Denis Ridzal (dridzal@sandia.gov), or
41 // Kara Peterson (kjpeter@sandia.gov)
42 //
43 // ************************************************************************
44 // @HEADER
45 
51 namespace Intrepid {
52 
53 
54 template<class Scalar, class ArrayScalar>
55 Basis_HGRAD_HEX_C2_FEM<Scalar,ArrayScalar>::Basis_HGRAD_HEX_C2_FEM()
56  {
57  this -> basisCardinality_ = 27;
58  this -> basisDegree_ = 2;
59  this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >() );
60  this -> basisType_ = BASIS_FEM_DEFAULT;
61  this -> basisCoordinates_ = COORDINATES_CARTESIAN;
62  this -> basisTagsAreSet_ = false;
63  }
64 
65 
66 
67 template<class Scalar, class ArrayScalar>
68 void Basis_HGRAD_HEX_C2_FEM<Scalar, ArrayScalar>::initializeTags() {
69 
70  // Basis-dependent intializations
71  int tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
72  int posScDim = 0; // position in the tag, counting from 0, of the subcell dim
73  int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
74  int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
75 
76  // An array with local DoF tags assigned to basis functions, in the order of their local enumeration
77  int tags[] = { 0, 0, 0, 1, // Nodes 0 to 7 follow vertex order of the topology
78  0, 1, 0, 1,
79  0, 2, 0, 1,
80  0, 3, 0, 1,
81  0, 4, 0, 1,
82  0, 5, 0, 1,
83  0, 6, 0, 1,
84  0, 7, 0, 1,
85  1, 0, 0, 1, // Node 8 -> edge 0
86  1, 1, 0, 1, // Node 9 -> edge 1
87  1, 2, 0, 1, // Node 10 -> edge 2
88  1, 3, 0, 1, // Node 11 -> edge 3
89  1, 8, 0, 1, // Node 12 -> edge 8
90  1, 9, 0, 1, // Node 13 -> edge 9
91  1,10, 0, 1, // Node 14 -> edge 10
92  1,11, 0, 1, // Node 15 -> edge 11
93  1, 4, 0, 1, // Node 16 -> edge 4
94  1, 5, 0, 1, // Node 17 -> edge 5
95  1, 6, 0, 1, // Node 18 -> edge 6
96  1, 7, 0, 1, // Node 19 -> edge 7
97  3, 0, 0, 1, // Node 20 -> Hexahedron
98  2, 4, 0, 1, // Node 21 -> face 4
99  2, 5, 0, 1, // Node 22 -> face 5
100  2, 3, 0, 1, // Node 23 -> face 3
101  2, 1, 0, 1, // Node 24 -> face 1
102  2, 0, 0, 1, // Node 25 -> face 0
103  2, 2, 0, 1, // Node 26 -> face 2
104  };
105 
106  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
107  Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
108  this -> ordinalToTag_,
109  tags,
110  this -> basisCardinality_,
111  tagSize,
112  posScDim,
113  posScOrd,
114  posDfOrd);
115 }
116 
117 
118 
119 template<class Scalar, class ArrayScalar>
120 void Basis_HGRAD_HEX_C2_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar & outputValues,
121  const ArrayScalar & inputPoints,
122  const EOperator operatorType) const {
123 
124  // Verify arguments
125 #ifdef HAVE_INTREPID_DEBUG
126  Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
127  inputPoints,
128  operatorType,
129  this -> getBaseCellTopology(),
130  this -> getCardinality() );
131 #endif
132 
133  // Number of evaluation points = dim 0 of inputPoints
134  int dim0 = inputPoints.dimension(0);
135 
136  // Temporaries: (x,y,z) coordinates of the evaluation point
137  Scalar x = 0.0;
138  Scalar y = 0.0;
139  Scalar z = 0.0;
140 
141  switch (operatorType) {
142 
143  case OPERATOR_VALUE:
144  for (int i0 = 0; i0 < dim0; i0++) {
145  x = inputPoints(i0, 0);
146  y = inputPoints(i0, 1);
147  z = inputPoints(i0, 2);
148 
149  // outputValues is a rank-2 array with dimensions (basisCardinality_, dim0)
150  outputValues( 0, i0) = 0.125*(-1. + x)*x*(-1. + y)*y*(-1. + z)*z;
151  outputValues( 1, i0) = 0.125*x*(1.+ x)*(-1. + y)*y*(-1. + z)*z;
152  outputValues( 2, i0) = 0.125*x*(1.+ x)*y*(1.+ y)*(-1. + z)*z;
153  outputValues( 3, i0) = 0.125*(-1. + x)*x*y*(1.+ y)*(-1. + z)*z;
154  outputValues( 4, i0) = 0.125*(-1. + x)*x*(-1. + y)*y*z*(1.+ z);
155  outputValues( 5, i0) = 0.125*x*(1.+ x)*(-1. + y)*y*z*(1.+ z);
156  outputValues( 6, i0) = 0.125*x*(1.+ x)*y*(1.+ y)*z*(1.+ z);
157  outputValues( 7, i0) = 0.125*(-1. + x)*x*y*(1.+ y)*z*(1.+ z);
158  outputValues( 8, i0) = 0.25*(1. - x)*(1. + x)*(-1. + y)*y*(-1. + z)*z;
159  outputValues( 9, i0) = 0.25*x*(1.+ x)*(1. - y)*(1. + y)*(-1. + z)*z;
160  outputValues(10, i0) = 0.25*(1. - x)*(1. + x)*y*(1.+ y)*(-1. + z)*z;
161  outputValues(11, i0) = 0.25*(-1. + x)*x*(1. - y)*(1. + y)*(-1. + z)*z;
162  outputValues(12, i0) = 0.25*(-1. + x)*x*(-1. + y)*y*(1. - z)*(1. + z);
163  outputValues(13, i0) = 0.25*x*(1.+ x)*(-1. + y)*y*(1. - z)*(1. + z);
164  outputValues(14, i0) = 0.25*x*(1.+ x)*y*(1.+ y)*(1. - z)*(1. + z);
165  outputValues(15, i0) = 0.25*(-1. + x)*x*y*(1.+ y)*(1. - z)*(1. + z);
166  outputValues(16, i0) = 0.25*(1. - x)*(1. + x)*(-1. + y)*y*z*(1.+ z);
167  outputValues(17, i0) = 0.25*x*(1.+ x)*(1. - y)*(1. + y)*z*(1.+ z);
168  outputValues(18, i0) = 0.25*(1. - x)*(1. + x)*y*(1.+ y)*z*(1.+ z);
169  outputValues(19, i0) = 0.25*(-1. + x)*x*(1. - y)*(1. + y)*z*(1.+ z);
170  outputValues(20, i0) = (1. - x)*(1. + x)*(1. - y)*(1. + y)*(1. - z)*(1. + z);
171  outputValues(21, i0) = 0.5*(1. - x)*(1. + x)*(1. - y)*(1. + y)*(-1. + z)*z;
172  outputValues(22, i0) = 0.5*(1. - x)*(1. + x)*(1. - y)*(1. + y)*z*(1.+ z);
173  outputValues(23, i0) = 0.5*(-1. + x)*x*(1. - y)*(1. + y)*(1. - z)*(1. + z);
174  outputValues(24, i0) = 0.5*x*(1.+ x)*(1. - y)*(1. + y)*(1. - z)*(1. + z);
175  outputValues(25, i0) = 0.5*(1. - x)*(1. + x)*(-1. + y)*y*(1. - z)*(1. + z);
176  outputValues(26, i0) = 0.5*(1. - x)*(1. + x)*y*(1.+ y)*(1. - z)*(1. + z);
177  }
178  break;
179 
180  case OPERATOR_GRAD:
181  case OPERATOR_D1:
182  for (int i0 = 0; i0 < dim0; i0++) {
183  x = inputPoints(i0,0);
184  y = inputPoints(i0,1);
185  z = inputPoints(i0,2);
186 
187  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
188  outputValues(0, i0, 0) = (-0.125 + 0.25*x)*(-1. + y)*y*(-1. + z)*z;
189  outputValues(0, i0, 1) = (-1. + x)*x*(-0.125 + 0.25*y)*(-1. + z)*z;
190  outputValues(0, i0, 2) = (-1. + x)*x*(-1. + y)*y*(-0.125 + 0.25*z);
191 
192  outputValues(1, i0, 0) = (0.125 + 0.25*x)*(-1. + y)*y*(-1. + z)*z;
193  outputValues(1, i0, 1) = x*(1. + x)*(-0.125 + 0.25*y)*(-1. + z)*z;
194  outputValues(1, i0, 2) = x*(1. + x)*(-1. + y)*y*(-0.125 + 0.25*z);
195 
196  outputValues(2, i0, 0) = (0.125 + 0.25*x)*y*(1. + y)*(-1. + z)*z;
197  outputValues(2, i0, 1) = x*(1. + x)*(0.125 + 0.25*y)*(-1. + z)*z;
198  outputValues(2, i0, 2) = x*(1. + x)*y*(1. + y)*(-0.125 + 0.25*z);
199 
200  outputValues(3, i0, 0) = (-0.125 + 0.25*x)*y*(1. + y)*(-1. + z)*z;
201  outputValues(3, i0, 1) = (-1. + x)*x*(0.125 + 0.25*y)*(-1. + z)*z;
202  outputValues(3, i0, 2) = (-1. + x)*x*y*(1. + y)*(-0.125 + 0.25*z);
203 
204  outputValues(4, i0, 0) = (-0.125 + 0.25*x)*(-1. + y)*y*z*(1. + z);
205  outputValues(4, i0, 1) = (-1. + x)*x*(-0.125 + 0.25*y)*z*(1. + z);
206  outputValues(4, i0, 2) = (-1. + x)*x*(-1. + y)*y*(0.125 + 0.25*z);
207 
208  outputValues(5, i0, 0) = (0.125 + 0.25*x)*(-1. + y)*y*z*(1. + z);
209  outputValues(5, i0, 1) = x*(1. + x)*(-0.125 + 0.25*y)*z*(1. + z);
210  outputValues(5, i0, 2) = x*(1. + x)*(-1. + y)*y*(0.125 + 0.25*z);
211 
212  outputValues(6, i0, 0) = (0.125 + 0.25*x)*y*(1. + y)*z*(1. + z);
213  outputValues(6, i0, 1) = x*(1. + x)*(0.125 + 0.25*y)*z*(1. + z);
214  outputValues(6, i0, 2) = x*(1. + x)*y*(1. + y)*(0.125 + 0.25*z);
215 
216  outputValues(7, i0, 0) = (-0.125 + 0.25*x)*y*(1. + y)*z*(1. + z);
217  outputValues(7, i0, 1) = (-1. + x)*x*(0.125 + 0.25*y)*z*(1. + z);
218  outputValues(7, i0, 2) = (-1. + x)*x*y*(1. + y)*(0.125 + 0.25*z);
219 
220  outputValues(8, i0, 0) = -0.5*x*(-1. + y)*y*(-1. + z)*z;
221  outputValues(8, i0, 1) = (1. - x)*(1. + x)*(-0.25 + 0.5*y)*(-1. + z)*z;
222  outputValues(8, i0, 2) = (1. - x)*(1. + x)*(-1. + y)*y*(-0.25 + 0.5*z);
223 
224  outputValues(9, i0, 0) = (0.25 + 0.5*x)*(1. - y)*(1. + y)*(-1. + z)*z;
225  outputValues(9, i0, 1) = x*(1. + x)*(-0.5*y)*(-1. + z)*z;
226  outputValues(9, i0, 2) = x*(1. + x)*(1. - y)*(1. + y)*(-0.25 + 0.5*z);
227 
228  outputValues(10,i0, 0) = -0.5*x*y*(1. + y)*(-1. + z)*z;
229  outputValues(10,i0, 1) = (1. - x)*(1. + x)*(0.25 + 0.5*y)*(-1. + z)*z;
230  outputValues(10,i0, 2) = (1. - x)*(1. + x)*y*(1. + y)*(-0.25 + 0.5*z);
231 
232  outputValues(11,i0, 0) = (-0.25 + 0.5*x)*(1. - y)*(1. + y)*(-1. + z)*z;
233  outputValues(11,i0, 1) = (-1. + x)*x*(-0.5*y)*(-1. + z)*z;
234  outputValues(11,i0, 2) = (-1. + x)*x*(1. - y)*(1. + y)*(-0.25 + 0.5*z);
235 
236  outputValues(12,i0, 0) = (-0.25 + 0.5*x)*(-1. + y)*y*(1. - z)*(1. + z);
237  outputValues(12,i0, 1) = (-1. + x)*x*(-0.25 + 0.5*y)*(1. - z)*(1. + z);
238  outputValues(12,i0, 2) = (-1. + x)*x*(-1. + y)*y*(-0.5*z);
239 
240  outputValues(13,i0, 0) = (0.25 + 0.5*x)*(-1. + y)*y*(1. - z)*(1. + z);
241  outputValues(13,i0, 1) = x*(1. + x)*(-0.25 + 0.5*y)*(1. - z)*(1. + z);
242  outputValues(13,i0, 2) = x*(1. + x)*(-1. + y)*y*(-0.5*z);
243 
244  outputValues(14,i0, 0) = (0.25 + 0.5*x)*y*(1. + y)*(1. - z)*(1. + z);
245  outputValues(14,i0, 1) = x*(1. + x)*(0.25 + 0.5*y)*(1. - z)*(1. + z);
246  outputValues(14,i0, 2) = x*(1. + x)*y*(1. + y)*(-0.5*z);
247 
248  outputValues(15,i0, 0) = (-0.25 + 0.5*x)*y*(1. + y)*(1. - z)*(1. + z);
249  outputValues(15,i0, 1) = (-1. + x)*x*(0.25 + 0.5*y)*(1. - z)*(1. + z);
250  outputValues(15,i0, 2) = (-1. + x)*x*y*(1. + y)*(-0.5*z);
251 
252  outputValues(16,i0, 0) = -0.5*x*(-1. + y)*y*z*(1. + z);
253  outputValues(16,i0, 1) = (1. - x)*(1. + x)*(-0.25 + 0.5*y)*z*(1. + z);
254  outputValues(16,i0, 2) = (1. - x)*(1. + x)*(-1. + y)*y*(0.25 + 0.5*z);
255 
256  outputValues(17,i0, 0) = (0.25 + 0.5*x)*(1. - y)*(1. + y)*z*(1. + z);
257  outputValues(17,i0, 1) = x*(1. + x)*(-0.5*y)*z*(1. + z);
258  outputValues(17,i0, 2) = x*(1. + x)*(1. - y)*(1. + y)*(0.25 + 0.5*z);
259 
260  outputValues(18,i0, 0) = -0.5*x*y*(1. + y)*z*(1. + z);
261  outputValues(18,i0, 1) = (1. - x)*(1. + x)*(0.25 + 0.5*y)*z*(1. + z);
262  outputValues(18,i0, 2) = (1. - x)*(1. + x)*y*(1. + y)*(0.25 + 0.5*z);
263 
264  outputValues(19,i0, 0) = (-0.25 + 0.5*x)*(1. - y)*(1. + y)*z*(1. + z);
265  outputValues(19,i0, 1) = (-1. + x)*x*(-0.5*y)*z*(1. + z);
266  outputValues(19,i0, 2) = (-1. + x)*x*(1. - y)*(1. + y)*(0.25 + 0.5*z);
267 
268  outputValues(20,i0, 0) = -2.*x*(1. - y)*(1. + y)*(1. - z)*(1. + z);
269  outputValues(20,i0, 1) = (1. - x)*(1. + x)*(-2.*y)*(1. - z)*(1. + z);
270  outputValues(20,i0, 2) = (1. - x)*(1. + x)*(1. - y)*(1. + y)*(-2.*z);
271 
272  outputValues(21,i0, 0) = -x*(1. - y)*(1. + y)*(-1. + z)*z;
273  outputValues(21,i0, 1) = (1. - x)*(1. + x)*(-y)*(-1. + z)*z;
274  outputValues(21,i0, 2) = (1. - x)*(1. + x)*(1. - y)*(1. + y)*(-0.5 + z);
275 
276  outputValues(22,i0, 0) = -x*(1. - y)*(1. + y)*z*(1. + z);
277  outputValues(22,i0, 1) = (1. - x)*(1. + x)*(-y)*z*(1. + z);
278  outputValues(22,i0, 2) = (1. - x)*(1. + x)*(1. - y)*(1. + y)*(0.5 + z);
279 
280  outputValues(23,i0, 0) = (-0.5 + x)*(1. - y)*(1. + y)*(1. - z)*(1. + z);
281  outputValues(23,i0, 1) = (-1. + x)*x*(-y)*(1. - z)*(1. + z);
282  outputValues(23,i0, 2) = (-1. + x)*x*(1. - y)*(1. + y)*(-z);
283 
284  outputValues(24,i0, 0) = (0.5 + x)*(1. - y)*(1. + y)*(1. - z)*(1. + z);
285  outputValues(24,i0, 1) = x*(1. + x)*(-y)*(1. - z)*(1. + z);
286  outputValues(24,i0, 2) = x*(1. + x)*(1. - y)*(1. + y)*(-z);
287 
288  outputValues(25,i0, 0) = -x*(-1. + y)*y*(1. - z)*(1. + z);
289  outputValues(25,i0, 1) = (1. - x)*(1. + x)*(-0.5 + y)*(1. - z)*(1. + z);
290  outputValues(25,i0, 2) = (1. - x)*(1. + x)*(-1. + y)*y*(-z);
291 
292  outputValues(26,i0, 0) = -x*y*(1. + y)*(1. - z)*(1. + z);
293  outputValues(26,i0, 1) = (1. - x)*(1. + x)*(0.5 + y)*(1. - z)*(1. + z);
294  outputValues(26,i0, 2) = (1. - x)*(1. + x)*y*(1. + y)*(-z);
295  }
296  break;
297 
298  case OPERATOR_CURL:
299  TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_CURL), std::invalid_argument,
300  ">>> ERROR (Basis_HGRAD_HEX_C2_FEM): CURL is invalid operator for rank-0 (scalar) functions in 3D");
301  break;
302 
303  case OPERATOR_DIV:
304  TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
305  ">>> ERROR (Basis_HGRAD_HEX_C2_FEM): DIV is invalid operator for rank-0 (scalar) functions in 3D");
306  break;
307 
308  case OPERATOR_D2:
309  for (int i0 = 0; i0 < dim0; i0++) {
310  x = inputPoints(i0,0);
311  y = inputPoints(i0,1);
312  z = inputPoints(i0,2);
313 
314  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, D2Cardinality = 6)
315  outputValues(0, i0, 0) = 0.25*(-1. + y)*y*(-1. + z)*z;
316  outputValues(0, i0, 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;
317  outputValues(0, i0, 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);
318  outputValues(0, i0, 3) = 0.25*(-1. + x)*x*(-1. + z)*z;
319  outputValues(0, i0, 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);
320  outputValues(0, i0, 5) = 0.25*(-1. + x)*x*(-1. + y)*y;
321 
322  outputValues(1, i0, 0) = 0.25*(-1. + y)*y*(-1. + z)*z;
323  outputValues(1, i0, 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;
324  outputValues(1, i0, 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);
325  outputValues(1, i0, 3) = 0.25*x*(1 + x)*(-1. + z)*z;
326  outputValues(1, i0, 4) = x*(1. + x)*(0.125 + y*(-0.25 + 0.5*z) - 0.25*z);
327  outputValues(1, i0, 5) = 0.25*x*(1 + x)*(-1. + y)*y;
328 
329  outputValues(2, i0, 0) = 0.25*y*(1 + y)*(-1. + z)*z;
330  outputValues(2, i0, 1) = (0.125 + x*(0.25 + 0.5*y) + 0.25*y)*(-1. + z)*z;
331  outputValues(2, i0, 2) = y*(1. + y)*(-0.125 + x*(-0.25 + 0.5*z) + 0.25*z);
332  outputValues(2, i0, 3) = 0.25*x*(1 + x)*(-1. + z)*z;
333  outputValues(2, i0, 4) = x*(1. + x)*(-0.125 + y*(-0.25 + 0.5*z) + 0.25*z);
334  outputValues(2, i0, 5) = 0.25*x*(1 + x)*y*(1 + y);
335 
336  outputValues(3, i0, 0) = 0.25*y*(1 + y)*(-1. + z)*z;
337  outputValues(3, i0, 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;
338  outputValues(3, i0, 2) = y*(1. + y)*(0.125 + x*(-0.25 + 0.5*z) - 0.25*z);
339  outputValues(3, i0, 3) = 0.25*(-1. + x)*x*(-1. + z)*z;
340  outputValues(3, i0, 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);
341  outputValues(3, i0, 5) = 0.25*(-1. + x)*x*y*(1 + y);
342 
343  outputValues(4, i0, 0) = 0.25*(-1. + y)*y*z*(1 + z);
344  outputValues(4, i0, 1) = (0.125 + x*(-0.25 + 0.5*y) - 0.25*y)*z*(1. + z);
345  outputValues(4, i0, 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);
346  outputValues(4, i0, 3) = 0.25*(-1. + x)*x*z*(1 + z);
347  outputValues(4, i0, 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);
348  outputValues(4, i0, 5) = 0.25*(-1. + x)*x*(-1. + y)*y;
349 
350  outputValues(5, i0, 0) = 0.25*(-1. + y)*y*z*(1 + z);
351  outputValues(5, i0, 1) = (-0.125 + x*(-0.25 + 0.5*y) + 0.25*y)*z*(1. + z);
352  outputValues(5, i0, 2) = (-1. + y)*y*(0.125 + x*(0.25 + 0.5*z) + 0.25*z);
353  outputValues(5, i0, 3) = 0.25*x*(1 + x)*z*(1 + z);
354  outputValues(5, i0, 4) = x*(1. + x)*(-0.125 + y*(0.25 + 0.5*z) - 0.25*z);
355  outputValues(5, i0, 5) = 0.25*x*(1 + x)*(-1. + y)*y;
356 
357  outputValues(6, i0, 0) = 0.25*y*(1 + y)*z*(1 + z);
358  outputValues(6, i0, 1) = (0.125 + x*(0.25 + 0.5*y) + 0.25*y)*z*(1. + z);
359  outputValues(6, i0, 2) = y*(1. + y)*(0.125 + x*(0.25 + 0.5*z) + 0.25*z);
360  outputValues(6, i0, 3) = 0.25*x*(1 + x)*z*(1 + z);
361  outputValues(6, i0, 4) = x*(1. + x)*(0.125 + y*(0.25 + 0.5*z) + 0.25*z);
362  outputValues(6, i0, 5) = 0.25*x*(1 + x)*y*(1 + y);
363 
364  outputValues(7, i0, 0) = 0.25*y*(1 + y)*z*(1 + z);
365  outputValues(7, i0, 1) = (-0.125 + x*(0.25 + 0.5*y) - 0.25*y)*z*(1. + z);
366  outputValues(7, i0, 2) = y*(1. + y)*(-0.125 + x*(0.25 + 0.5*z) - 0.25*z);
367  outputValues(7, i0, 3) = 0.25*(-1. + x)*x*z*(1 + z);
368  outputValues(7, i0, 4) = (-1. + x)*x*(0.125 + y*(0.25 + 0.5*z) + 0.25*z);
369  outputValues(7, i0, 5) = 0.25*(-1. + x)*x*y*(1 + y);
370 
371  outputValues(8, i0, 0) = -0.5*(-1. + y)*y*(-1. + z)*z;
372  outputValues(8, i0, 1) = (0. + x*(-0.5 + y))*z + (x*(0.5 - y) )*(z*z);
373  outputValues(8, i0, 2) = (y*y)*(x*(0.5 - z) ) + y*(x*(-0.5 + z));
374  outputValues(8, i0, 3) = 0.5*(1. - x)*(1. + x)*(-1. + z)*z;
375  outputValues(8, i0, 4) = 0.25 + (x*x)*(-0.25 + y*(0.5 - z) + 0.5*z) - 0.5*z + y*(-0.5 + z);
376  outputValues(8, i0, 5) = 0.5*(1. - x)*(1. + x)*(-1. + y)*y;
377 
378  outputValues(9, i0, 0) = 0.5*(1. - y)*(1. + y)*(-1. + z)*z;
379  outputValues(9, i0, 1) = (0.5*y + x*(y))*z + (x*(-y) - 0.5*y)*(z*z);
380  outputValues(9, i0, 2) = -0.25 + (y*y)*(0.25 - 0.5*z) + 0.5*z + x*(-0.5 + (y*y)*(0.5 - z) + z);
381  outputValues(9, i0, 3) = -0.5*x*(1 + x)*(-1. + z)*z;
382  outputValues(9, i0, 4) = x*(y*(0.5 - z) ) + (x*x)*(y*(0.5 - z) );
383  outputValues(9, i0, 5) = 0.5*x*(1 + x)*(1. - y)*(1. + y);
384 
385  outputValues(10,i0, 0) = -0.5*y*(1 + y)*(-1. + z)*z;
386  outputValues(10,i0, 1) = (0. + x*(0.5 + y))*z + (x*(-0.5 - y) )*(z*z);
387  outputValues(10,i0, 2) = y*(x*(0.5 - z) ) + (y*y)*(x*(0.5 - z) );
388  outputValues(10,i0, 3) = 0.5*(1. - x)*(1. + x)*(-1. + z)*z;
389  outputValues(10,i0, 4) = -0.25 + (x*x)*(0.25 + y*(0.5 - z) - 0.5*z) + 0.5*z + y*(-0.5 + z);
390  outputValues(10,i0, 5) = 0.5*(1. - x)*(1. + x)*y*(1 + y);
391 
392  outputValues(11,i0, 0) = 0.5*(1. - y)*(1. + y)*(-1. + z)*z;
393  outputValues(11,i0, 1) = (-0.5*y + x*(y))*z + (x*(-y) + 0.5*y)*(z*z);
394  outputValues(11,i0, 2) = 0.25 + (y*y)*(-0.25 + 0.5*z) - 0.5*z + x*(-0.5 + (y*y)*(0.5 - z) + z);
395  outputValues(11,i0, 3) = -0.5*(-1. + x)*x*(-1. + z)*z;
396  outputValues(11,i0, 4) = (x*x)*(y*(0.5 - z) ) + x*(y*(-0.5 + z));
397  outputValues(11,i0, 5) = 0.5*(-1. + x)*x*(1. - y)*(1. + y);
398 
399  outputValues(12,i0, 0) = 0.5*(-1. + y)*y*(1. - z)*(1. + z);
400  outputValues(12,i0, 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)));
401  outputValues(12,i0, 2) = (y*y)*(x*(-z) + 0.5*z) + y*(-0.5*z + x*(z));
402  outputValues(12,i0, 3) = 0.5*(-1. + x)*x*(1. - z)*(1. + z);
403  outputValues(12,i0, 4) = (x*x)*(y*(-z) + 0.5*z) + x*(-0.5*z + y*(z));
404  outputValues(12,i0, 5) = -0.5*(-1. + x)*x*(-1. + y)*y;
405 
406  outputValues(13,i0, 0) = 0.5*(-1. + y)*y*(1. - z)*(1. + z);
407  outputValues(13,i0, 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)));
408  outputValues(13,i0, 2) = (y*y)*(x*(-z) - 0.5*z) + y*(0.5*z + x*(z));
409  outputValues(13,i0, 3) = 0.5*x*(1 + x)*(1. - z)*(1. + z);
410  outputValues(13,i0, 4) = x*(y*(-z) + 0.5*z) + (x*x)*(y*(-z) + 0.5*z);
411  outputValues(13,i0, 5) = -0.5*x*(1 + x)*(-1. + y)*y;
412 
413  outputValues(14,i0, 0) = 0.5*y*(1 + y)*(1. - z)*(1. + z);
414  outputValues(14,i0, 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)));
415  outputValues(14,i0, 2) = y*(x*(-z) - 0.5*z) + (y*y)*(x*(-z) - 0.5*z);
416  outputValues(14,i0, 3) = 0.5*x*(1 + x)*(1. - z)*(1. + z);
417  outputValues(14,i0, 4) = x*(y*(-z) - 0.5*z) + (x*x)*(y*(-z) - 0.5*z);
418  outputValues(14,i0, 5) = -0.5*x*(1 + x)*y*(1 + y);
419 
420  outputValues(15,i0, 0) = 0.5*y*(1 + y)*(1. - z)*(1. + z);
421  outputValues(15,i0, 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)));
422  outputValues(15,i0, 2) = y*(x*(-z) + 0.5*z) + (y*y)*(x*(-z) + 0.5*z);
423  outputValues(15,i0, 3) = 0.5*(-1. + x)*x*(1. - z)*(1. + z);
424  outputValues(15,i0, 4) = (x*x)*(y*(-z) - 0.5*z) + x*(0.5*z + y*(z));
425  outputValues(15,i0, 5) = -0.5*(-1. + x)*x*y*(1 + y);
426 
427  outputValues(16,i0, 0) = -0.5*(-1. + y)*y*z*(1 + z);
428  outputValues(16,i0, 1) = (x*(0.5 - y) )*z + (x*(0.5 - y) )*(z*z);
429  outputValues(16,i0, 2) = (y*y)*(x*(-0.5 - z) ) + y*(x*(0.5 + z));
430  outputValues(16,i0, 3) = 0.5*(1. - x)*(1. + x)*z*(1 + z);
431  outputValues(16,i0, 4) = -0.25 + (x*x)*(0.25 + y*(-0.5 - z) + 0.5*z) - 0.5*z + y*(0.5 + z);
432  outputValues(16,i0, 5) = 0.5*(1. - x)*(1. + x)*(-1. + y)*y;
433 
434  outputValues(17,i0, 0) = 0.5*(1. - y)*(1. + y)*z*(1 + z);
435  outputValues(17,i0, 1) = (x*(-y) - 0.5*y)*z + (x*(-y) - 0.5*y)*(z*z);
436  outputValues(17,i0, 2) = 0.25 + (y*y)*(-0.25 - 0.5*z) + 0.5*z + x*(0.5 + (y*y)*(-0.5 - z) + z);
437  outputValues(17,i0, 3) = -0.5*x*(1 + x)*z*(1 + z);
438  outputValues(17,i0, 4) = x*(y*(-0.5 - z) ) + (x*x)*(y*(-0.5 - z) );
439  outputValues(17,i0, 5) = 0.5*x*(1 + x)*(1. - y)*(1. + y);
440 
441  outputValues(18,i0, 0) = -0.5*y*(1 + y)*z*(1 + z);
442  outputValues(18,i0, 1) = (x*(-0.5 - y) )*z + (x*(-0.5 - y) )*(z*z);
443  outputValues(18,i0, 2) = y*(x*(-0.5 - z) ) + (y*y)*(x*(-0.5 - z) );
444  outputValues(18,i0, 3) = 0.5*(1. - x)*(1. + x)*z*(1 + z);
445  outputValues(18,i0, 4) = 0.25 + (x*x)*(-0.25 + y*(-0.5 - z) - 0.5*z) + 0.5*z + y*(0.5 + z);
446  outputValues(18,i0, 5) = 0.5*(1. - x)*(1. + x)*y*(1 + y);
447 
448  outputValues(19,i0, 0) = 0.5*(1. - y)*(1. + y)*z*(1 + z);
449  outputValues(19,i0, 1) = (x*(-y) + 0.5*y)*z + (x*(-y) + 0.5*y)*(z*z);
450  outputValues(19,i0, 2) = -0.25 + (y*y)*(0.25 + 0.5*z) - 0.5*z + x*(0.5 + (y*y)*(-0.5 - z) + z);
451  outputValues(19,i0, 3) = -0.5*(-1. + x)*x*z*(1 + z);
452  outputValues(19,i0, 4) = (x*x)*(y*(-0.5 - z) ) + x*(y*(0.5 + z));
453  outputValues(19,i0, 5) = 0.5*(-1. + x)*x*(1. - y)*(1. + y);
454 
455  outputValues(20,i0, 0) = -2.*(1. - y)*(1. + y)*(1. - z)*(1. + z);
456  outputValues(20,i0, 1) = -4.*x*y*(-1. + z*z);
457  outputValues(20,i0, 2) = x*((y*y)*(-4.*z) + 4.*z);
458  outputValues(20,i0, 3) = -2.*(1. - x)*(1. + x)*(1. - z)*(1. + z);
459  outputValues(20,i0, 4) = (x*x)*(y*(-4.*z) ) + y*(4.*z);
460  outputValues(20,i0, 5) = -2.*(1. - x)*(1. + x)*(1. - y)*(1. + y);
461 
462  outputValues(21,i0, 0) = -(1. - y)*(1. + y)*(-1. + z)*z;
463  outputValues(21,i0, 1) = (x*(-2.*y) )*z + (0. + x*(2.*y))*(z*z);
464  outputValues(21,i0, 2) = x*(1. - 2.*z + (y*y)*(-1. + 2.*z));
465  outputValues(21,i0, 3) = -(1. - x)*(1. + x)*(-1. + z)*z;
466  outputValues(21,i0, 4) = y*(1. - 2.*z) + (x*x)*(y*(-1. + 2.*z));
467  outputValues(21,i0, 5) = (1. - x)*(1. + x)*(1. - y)*(1. + y);
468 
469  outputValues(22,i0, 0) = -(1. - y)*(1. + y)*z*(1 + z);
470  outputValues(22,i0, 1) = (0. + x*(2.*y))*z + (0. + x*(2.*y))*(z*z);
471  outputValues(22,i0, 2) = x*(-1. - 2.*z + (y*y)*(1. + 2.*z));
472  outputValues(22,i0, 3) = -(1. - x)*(1. + x)*z*(1 + z);
473  outputValues(22,i0, 4) = y*(-1. - 2.*z) + (x*x)*(y*(1. + 2.*z));
474  outputValues(22,i0, 5) = (1. - x)*(1. + x)*(1. - y)*(1. + y);
475 
476  outputValues(23,i0, 0) = (1. - y)*(1. + y)*(1. - z)*(1. + z);
477  outputValues(23,i0, 1) = (-1. + 2.*x)*y*(-1. + z*z);
478  outputValues(23,i0, 2) = (-1. + 2.*x)*(-1. + y*y)*z;
479  outputValues(23,i0, 3) =-(-1. + x)*x*(1. - z)*(1. + z);
480  outputValues(23,i0, 4) = 2.*(-1. + x)*x*y*z;
481  outputValues(23,i0, 5) =-(-1. + x)*x*(1. - y)*(1. + y);
482 
483  outputValues(24,i0, 0) = (1. - y)*(1. + y)*(1. - z)*(1. + z);
484  outputValues(24,i0, 1) = (1. + 2.*x)*y*(-1. + z*z);
485  outputValues(24,i0, 2) = (1. + 2.*x)*(-1. + y*y)*z;
486  outputValues(24,i0, 3) = x*(1. + x)*(-1. + z)*(1. + z);
487  outputValues(24,i0, 4) = 2.*x*(1. + x)*y*z;
488  outputValues(24,i0, 5) = x*(1. + x)*(-1. + y)*(1. + y);
489 
490  outputValues(25,i0, 0) = -(-1. + y)*y*(1. - z)*(1. + z);
491  outputValues(25,i0, 1) = x*(-1. + 2.*y)*(-1. + z*z);
492  outputValues(25,i0, 2) = 2.*x*(-1. + y)*y*z;
493  outputValues(25,i0, 3) = (1. - x)*(1. + x)*(1. - z)*(1. + z);
494  outputValues(25,i0, 4) = (-1. + x*x)*(-1. + 2.*y)*z;
495  outputValues(25,i0, 5) =-(1. - x)*(1. + x)*(-1. + y)*y;
496 
497  outputValues(26,i0, 0) = y*(1. + y)*(-1. + z)*(1. + z);
498  outputValues(26,i0, 1) = x*(1. + 2.*y)*(-1. + z*z);
499  outputValues(26,i0, 2) = 2.*x*y*(1. + y)*z;
500  outputValues(26,i0, 3) = (-1. + x)*(1. + x)*(-1. + z)*(1. + z);
501  outputValues(26,i0, 4) = (-1. + x*x)*(1. + 2.*y)*z;
502  outputValues(26,i0, 5) = (-1. + x)*(1. + x)*y*(1. + y);
503  }
504  break;
505 
506  case OPERATOR_D3:
507  for (int i0 = 0; i0 < dim0; i0++) {
508  x = inputPoints(i0,0);
509  y = inputPoints(i0,1);
510  z = inputPoints(i0,2);
511 
512  outputValues(0,i0, 0) = 0.;
513  outputValues(0,i0, 1) = ((-1.+ 2.*y)*(-1.+ z)*z)/4.;
514  outputValues(0,i0, 2) = ((-1.+ y)*y*(-1.+ 2.*z))/4.;
515  outputValues(0,i0, 3) = ((-1.+ 2.*x)*(-1.+ z)*z)/4.;
516  outputValues(0,i0, 4) = ((-1.+ 2.*x)*(-1.+ 2.*y)*(-1.+ 2.*z))/8.;
517  outputValues(0,i0, 5) = ((-1.+ 2.*x)*(-1.+ y)*y)/4.;
518  outputValues(0,i0, 6) = 0.;
519  outputValues(0,i0, 7) = ((-1.+ x)*x*(-1.+ 2.*z))/4.;
520  outputValues(0,i0, 8) = ((-1.+ x)*x*(-1.+ 2.*y))/4.;
521  outputValues(0,i0, 9) = 0.;
522 
523  outputValues(1, i0, 0) = 0.;
524  outputValues(1, i0, 1) = ((-1.+ 2.*y)*(-1.+ z)*z)/4.;
525  outputValues(1, i0, 2) = ((-1.+ y)*y*(-1.+ 2.*z))/4.;
526  outputValues(1, i0, 3) = ((1.+ 2.*x)*(-1.+ z)*z)/4.;
527  outputValues(1, i0, 4) = ((1.+ 2.*x)*(-1.+ 2.*y)*(-1.+ 2.*z))/8.;
528  outputValues(1, i0, 5) = ((1.+ 2.*x)*(-1.+ y)*y)/4.;
529  outputValues(1, i0, 6) = 0.;
530  outputValues(1, i0, 7) = (x*(1.+ x)*(-1.+ 2.*z))/4.;
531  outputValues(1, i0, 8) = (x*(1.+ x)*(-1.+ 2.*y))/4.;
532  outputValues(1, i0, 9) = 0.;
533 
534  outputValues(2, i0, 0) = 0.;
535  outputValues(2, i0, 1) = ((1.+ 2.*y)*(-1.+ z)*z)/4.;
536  outputValues(2, i0, 2) = (y*(1.+ y)*(-1.+ 2.*z))/4.;
537  outputValues(2, i0, 3) = ((1.+ 2.*x)*(-1.+ z)*z)/4.;
538  outputValues(2, i0, 4) = ((1.+ 2.*x)*(1.+ 2.*y)*(-1.+ 2.*z))/8.;
539  outputValues(2, i0, 5) = ((1.+ 2.*x)*y*(1.+ y))/4.;
540  outputValues(2, i0, 6) = 0.;
541  outputValues(2, i0, 7) = (x*(1.+ x)*(-1.+ 2.*z))/4.;
542  outputValues(2, i0, 8) = (x*(1.+ x)*(1.+ 2.*y))/4.;
543  outputValues(2, i0, 9) = 0.;
544 
545  outputValues(3, i0, 0) = 0.;
546  outputValues(3, i0, 1) = ((1.+ 2.*y)*(-1.+ z)*z)/4.;
547  outputValues(3, i0, 2) = (y*(1.+ y)*(-1.+ 2.*z))/4.;
548  outputValues(3, i0, 3) = ((-1.+ 2.*x)*(-1.+ z)*z)/4.;
549  outputValues(3, i0, 4) = ((-1.+ 2.*x)*(1.+ 2.*y)*(-1.+ 2.*z))/8.;
550  outputValues(3, i0, 5) = ((-1.+ 2.*x)*y*(1.+ y))/4.;
551  outputValues(3, i0, 6) = 0.;
552  outputValues(3, i0, 7) = ((-1.+ x)*x*(-1.+ 2.*z))/4.;
553  outputValues(3, i0, 8) = ((-1.+ x)*x*(1.+ 2.*y))/4.;
554  outputValues(3, i0, 9) = 0.;
555 
556  outputValues(4, i0, 0) = 0.;
557  outputValues(4, i0, 1) = ((-1.+ 2.*y)*z*(1.+ z))/4.;
558  outputValues(4, i0, 2) = ((-1.+ y)*y*(1.+ 2.*z))/4.;
559  outputValues(4, i0, 3) = ((-1.+ 2.*x)*z*(1.+ z))/4.;
560  outputValues(4, i0, 4) = ((-1.+ 2.*x)*(-1.+ 2.*y)*(1.+ 2.*z))/8.;
561  outputValues(4, i0, 5) = ((-1.+ 2.*x)*(-1.+ y)*y)/4.;
562  outputValues(4, i0, 6) = 0.;
563  outputValues(4, i0, 7) = ((-1.+ x)*x*(1.+ 2.*z))/4.;
564  outputValues(4, i0, 8) = ((-1.+ x)*x*(-1.+ 2.*y))/4.;
565  outputValues(4, i0, 9) = 0.;
566 
567  outputValues(5, i0, 0) = 0.;
568  outputValues(5, i0, 1) = ((-1.+ 2.*y)*z*(1.+ z))/4.;
569  outputValues(5, i0, 2) = ((-1.+ y)*y*(1.+ 2.*z))/4.;
570  outputValues(5, i0, 3) = ((1.+ 2.*x)*z*(1.+ z))/4.;
571  outputValues(5, i0, 4) = ((1.+ 2.*x)*(-1.+ 2.*y)*(1.+ 2.*z))/8.;
572  outputValues(5, i0, 5) = ((1.+ 2.*x)*(-1.+ y)*y)/4.;
573  outputValues(5, i0, 6) = 0.;
574  outputValues(5, i0, 7) = (x*(1.+ x)*(1.+ 2.*z))/4.;
575  outputValues(5, i0, 8) = (x*(1.+ x)*(-1.+ 2.*y))/4.;
576  outputValues(5, i0, 9) = 0.;
577 
578  outputValues(6, i0, 0) = 0.;
579  outputValues(6, i0, 1) = ((1.+ 2.*y)*z*(1.+ z))/4.;
580  outputValues(6, i0, 2) = (y*(1.+ y)*(1.+ 2.*z))/4.;
581  outputValues(6, i0, 3) = ((1.+ 2.*x)*z*(1.+ z))/4.;
582  outputValues(6, i0, 4) = ((1.+ 2.*x)*(1.+ 2.*y)*(1.+ 2.*z))/8.;
583  outputValues(6, i0, 5) = ((1.+ 2.*x)*y*(1.+ y))/4.;
584  outputValues(6, i0, 6) = 0.;
585  outputValues(6, i0, 7) = (x*(1.+ x)*(1.+ 2.*z))/4.;
586  outputValues(6, i0, 8) = (x*(1.+ x)*(1.+ 2.*y))/4.;
587  outputValues(6, i0, 9) = 0.;
588 
589  outputValues(7, i0, 0) = 0.;
590  outputValues(7, i0, 1) = ((1.+ 2.*y)*z*(1.+ z))/4.;
591  outputValues(7, i0, 2) = (y*(1.+ y)*(1.+ 2.*z))/4.;
592  outputValues(7, i0, 3) = ((-1.+ 2.*x)*z*(1.+ z))/4.;
593  outputValues(7, i0, 4) = ((-1.+ 2.*x)*(1.+ 2.*y)*(1.+ 2.*z))/8.;
594  outputValues(7, i0, 5) = ((-1.+ 2.*x)*y*(1.+ y))/4.;
595  outputValues(7, i0, 6) = 0.;
596  outputValues(7, i0, 7) = ((-1.+ x)*x*(1.+ 2.*z))/4.;
597  outputValues(7, i0, 8) = ((-1.+ x)*x*(1.+ 2.*y))/4.;
598  outputValues(7, i0, 9) = 0.;
599 
600  outputValues(8, i0, 0) = 0.;
601  outputValues(8, i0, 1) = -((-1.+ 2.*y)*(-1.+ z)*z)/2.;
602  outputValues(8, i0, 2) = -((-1.+ y)*y*(-1.+ 2.*z))/2.;
603  outputValues(8, i0, 3) = -(x*(-1.+ z)*z);
604  outputValues(8, i0, 4) = -(x*(-1.+ 2.*y)*(-1.+ 2.*z))/2.;
605  outputValues(8, i0, 5) = -(x*(-1.+ y)*y);
606  outputValues(8, i0, 6) = 0.;
607  outputValues(8, i0, 7) = -((-1.+ (x*x))*(-1.+ 2.*z))/2.;
608  outputValues(8, i0, 8) = -((-1.+ (x*x))*(-1.+ 2.*y))/2.;
609  outputValues(8, i0, 9) = 0.;
610 
611  outputValues(9, i0, 0) = 0.;
612  outputValues(9, i0, 1) = -(y*(-1.+ z)*z);
613  outputValues(9, i0, 2) = -((-1.+ (y*y))*(-1.+ 2.*z))/2.;
614  outputValues(9, i0, 3) = -((1.+ 2.*x)*(-1.+ z)*z)/2.;
615  outputValues(9, i0, 4) = -((1.+ 2.*x)*y*(-1.+ 2.*z))/2.;
616  outputValues(9, i0, 5) = -((1.+ 2.*x)*(-1.+ (y*y)))/2.;
617  outputValues(9, i0, 6) = 0.;
618  outputValues(9, i0, 7) = -(x*(1.+ x)*(-1.+ 2.*z))/2.;
619  outputValues(9, i0, 8) = -(x*(1.+ x)*y);
620  outputValues(9, i0, 9) = 0.;
621 
622  outputValues(10,i0, 0) = 0.;
623  outputValues(10,i0, 1) = -((1.+ 2.*y)*(-1.+ z)*z)/2.;
624  outputValues(10,i0, 2) = -(y*(1.+ y)*(-1.+ 2.*z))/2.;
625  outputValues(10,i0, 3) = -(x*(-1.+ z)*z);
626  outputValues(10,i0, 4) = -(x*(1.+ 2.*y)*(-1.+ 2.*z))/2.;
627  outputValues(10,i0, 5) = -(x*y*(1.+ y));
628  outputValues(10,i0, 6) = 0.;
629  outputValues(10,i0, 7) = -((-1.+ (x*x))*(-1.+ 2.*z))/2.;
630  outputValues(10,i0, 8) = -((-1.+ (x*x))*(1.+ 2.*y))/2.;
631  outputValues(10,i0, 9) = 0.;
632 
633  outputValues(11,i0, 0) = 0.;
634  outputValues(11,i0, 1) = -(y*(-1.+ z)*z);
635  outputValues(11,i0, 2) = -((-1.+ (y*y))*(-1.+ 2.*z))/2.;
636  outputValues(11,i0, 3) = -((-1.+ 2.*x)*(-1.+ z)*z)/2.;
637  outputValues(11,i0, 4) = -((-1.+ 2.*x)*y*(-1.+ 2.*z))/2.;
638  outputValues(11,i0, 5) = -((-1.+ 2.*x)*(-1.+ (y*y)))/2.;
639  outputValues(11,i0, 6) = 0.;
640  outputValues(11,i0, 7) = -((-1.+ x)*x*(-1.+ 2.*z))/2.;
641  outputValues(11,i0, 8) = -((-1.+ x)*x*y);
642  outputValues(11,i0, 9) = 0.;
643 
644  outputValues(12,i0, 0) = 0.;
645  outputValues(12,i0, 1) = -((-1.+ 2.*y)*(-1.+ (z*z)))/2.;
646  outputValues(12,i0, 2) = -((-1.+ y)*y*z);
647  outputValues(12,i0, 3) = -((-1.+ 2.*x)*(-1.+ (z*z)))/2.;
648  outputValues(12,i0, 4) = -((-1.+ 2.*x)*(-1.+ 2.*y)*z)/2.;
649  outputValues(12,i0, 5) = -((-1.+ 2.*x)*(-1.+ y)*y)/2.;
650  outputValues(12,i0, 6) = 0.;
651  outputValues(12,i0, 7) = -((-1.+ x)*x*z);
652  outputValues(12,i0, 8) = -((-1.+ x)*x*(-1.+ 2.*y))/2.;
653  outputValues(12,i0, 9) = 0.;
654 
655  outputValues(13,i0, 0) = 0.;
656  outputValues(13,i0, 1) = -((-1.+ 2.*y)*(-1.+ (z*z)))/2.;
657  outputValues(13,i0, 2) = -((-1.+ y)*y*z);
658  outputValues(13,i0, 3) = -((1.+ 2.*x)*(-1.+ (z*z)))/2.;
659  outputValues(13,i0, 4) = -((1.+ 2.*x)*(-1.+ 2.*y)*z)/2.;
660  outputValues(13,i0, 5) = -((1.+ 2.*x)*(-1.+ y)*y)/2.;
661  outputValues(13,i0, 6) = 0.;
662  outputValues(13,i0, 7) = -(x*(1.+ x)*z);
663  outputValues(13,i0, 8) = -(x*(1.+ x)*(-1.+ 2.*y))/2.;
664  outputValues(13,i0, 9) = 0.;
665 
666  outputValues(14,i0, 0) = 0.;
667  outputValues(14,i0, 1) = -((1.+ 2.*y)*(-1.+ (z*z)))/2.;
668  outputValues(14,i0, 2) = -(y*(1.+ y)*z);
669  outputValues(14,i0, 3) = -((1.+ 2.*x)*(-1.+ (z*z)))/2.;
670  outputValues(14,i0, 4) = -((1.+ 2.*x)*(1.+ 2.*y)*z)/2.;
671  outputValues(14,i0, 5) = -((1.+ 2.*x)*y*(1.+ y))/2.;
672  outputValues(14,i0, 6) = 0.;
673  outputValues(14,i0, 7) = -(x*(1.+ x)*z);
674  outputValues(14,i0, 8) = -(x*(1.+ x)*(1.+ 2.*y))/2.;
675  outputValues(14,i0, 9) = 0.;
676 
677  outputValues(15,i0, 0) = 0.;
678  outputValues(15,i0, 1) = -((1.+ 2.*y)*(-1.+ (z*z)))/2.;
679  outputValues(15,i0, 2) = -(y*(1.+ y)*z);
680  outputValues(15,i0, 3) = -((-1.+ 2.*x)*(-1.+ (z*z)))/2.;
681  outputValues(15,i0, 4) = -((-1.+ 2.*x)*(1.+ 2.*y)*z)/2.;
682  outputValues(15,i0, 5) = -((-1.+ 2.*x)*y*(1.+ y))/2.;
683  outputValues(15,i0, 6) = 0.;
684  outputValues(15,i0, 7) = -((-1.+ x)*x*z);
685  outputValues(15,i0, 8) = -((-1.+ x)*x*(1.+ 2.*y))/2.;
686  outputValues(15,i0, 9) = 0.;
687 
688  outputValues(16,i0, 0) = 0.;
689  outputValues(16,i0, 1) = -((-1.+ 2.*y)*z*(1.+ z))/2.;
690  outputValues(16,i0, 2) = -((-1.+ y)*y*(1.+ 2.*z))/2.;
691  outputValues(16,i0, 3) = -(x*z*(1.+ z));
692  outputValues(16,i0, 4) = -(x*(-1.+ 2.*y)*(1.+ 2.*z))/2.;
693  outputValues(16,i0, 5) = -(x*(-1.+ y)*y);
694  outputValues(16,i0, 6) = 0.;
695  outputValues(16,i0, 7) = -((-1.+ (x*x))*(1.+ 2.*z))/2.;
696  outputValues(16,i0, 8) = -((-1.+ (x*x))*(-1.+ 2.*y))/2.;
697  outputValues(16,i0, 9) = 0.;
698 
699  outputValues(17,i0, 0) = 0.;
700  outputValues(17,i0, 1) = -(y*z*(1.+ z));
701  outputValues(17,i0, 2) = -((-1.+ (y*y))*(1.+ 2.*z))/2.;
702  outputValues(17,i0, 3) = -((1.+ 2.*x)*z*(1.+ z))/2.;
703  outputValues(17,i0, 4) = -((1.+ 2.*x)*y*(1.+ 2.*z))/2.;
704  outputValues(17,i0, 5) = -((1.+ 2.*x)*(-1.+ (y*y)))/2.;
705  outputValues(17,i0, 6) = 0.;
706  outputValues(17,i0, 7) = -(x*(1.+ x)*(1.+ 2.*z))/2.;
707  outputValues(17,i0, 8) = -(x*(1.+ x)*y);
708  outputValues(17,i0, 9) = 0.;
709 
710  outputValues(18,i0, 0) = 0.;
711  outputValues(18,i0, 1) = -((1.+ 2.*y)*z*(1.+ z))/2.;
712  outputValues(18,i0, 2) = -(y*(1.+ y)*(1.+ 2.*z))/2.;
713  outputValues(18,i0, 3) = -(x*z*(1.+ z));
714  outputValues(18,i0, 4) = -(x*(1.+ 2.*y)*(1.+ 2.*z))/2.;
715  outputValues(18,i0, 5) = -(x*y*(1.+ y));
716  outputValues(18,i0, 6) = 0.;
717  outputValues(18,i0, 7) = -((-1.+ (x*x))*(1.+ 2.*z))/2.;
718  outputValues(18,i0, 8) = -((-1.+ (x*x))*(1.+ 2.*y))/2.;
719  outputValues(18,i0, 9) = 0.;
720 
721  outputValues(19,i0, 0) = 0.;
722  outputValues(19,i0, 1) = -(y*z*(1.+ z));
723  outputValues(19,i0, 2) = -((-1.+ (y*y))*(1.+ 2.*z))/2.;
724  outputValues(19,i0, 3) = -((-1.+ 2.*x)*z*(1.+ z))/2.;
725  outputValues(19,i0, 4) = -((-1.+ 2.*x)*y*(1.+ 2.*z))/2.;
726  outputValues(19,i0, 5) = -((-1.+ 2.*x)*(-1.+ (y*y)))/2.;
727  outputValues(19,i0, 6) = 0.;
728  outputValues(19,i0, 7) = -((-1.+ x)*x*(1.+ 2.*z))/2.;
729  outputValues(19,i0, 8) = -((-1.+ x)*x*y);
730  outputValues(19,i0, 9) = 0.;
731 
732  outputValues(20,i0, 0) = 0.;
733  outputValues(20,i0, 1) = -4*y*(-1.+ (z*z));
734  outputValues(20,i0, 2) = -4*(-1.+ (y*y))*z;
735  outputValues(20,i0, 3) = -4*x*(-1.+ (z*z));
736  outputValues(20,i0, 4) = -8*x*y*z;
737  outputValues(20,i0, 5) = -4*x*(-1.+ (y*y));
738  outputValues(20,i0, 6) = 0.;
739  outputValues(20,i0, 7) = -4*(-1.+ (x*x))*z;
740  outputValues(20,i0, 8) = -4*(-1.+ (x*x))*y;
741  outputValues(20,i0, 9) = 0.;
742 
743  outputValues(21,i0, 0) = 0.;
744  outputValues(21,i0, 1) = 2.*y*(-1.+ z)*z;
745  outputValues(21,i0, 2) = (-1.+ (y*y))*(-1.+ 2.*z);
746  outputValues(21,i0, 3) = 2.*x*(-1.+ z)*z;
747  outputValues(21,i0, 4) = 2.*x*y*(-1.+ 2.*z);
748  outputValues(21,i0, 5) = 2.*x*(-1.+ (y*y));
749  outputValues(21,i0, 6) = 0.;
750  outputValues(21,i0, 7) = (-1.+ (x*x))*(-1.+ 2.*z);
751  outputValues(21,i0, 8) = 2.*(-1.+ (x*x))*y;
752  outputValues(21,i0, 9) = 0.;
753 
754  outputValues(22,i0, 0) = 0.;
755  outputValues(22,i0, 1) = 2.*y*z*(1.+ z);
756  outputValues(22,i0, 2) = (-1.+ (y*y))*(1.+ 2.*z);
757  outputValues(22,i0, 3) = 2.*x*z*(1.+ z);
758  outputValues(22,i0, 4) = 2.*x*y*(1.+ 2.*z);
759  outputValues(22,i0, 5) = 2.*x*(-1.+ (y*y));
760  outputValues(22,i0, 6) = 0.;
761  outputValues(22,i0, 7) = (-1.+ (x*x))*(1.+ 2.*z);
762  outputValues(22,i0, 8) = 2.*(-1.+ (x*x))*y;
763  outputValues(22,i0, 9) = 0.;
764 
765  outputValues(23,i0, 0) = 0.;
766  outputValues(23,i0, 1) = 2.*y*(-1.+ (z*z));
767  outputValues(23,i0, 2) = 2.*(-1.+ (y*y))*z;
768  outputValues(23,i0, 3) = (-1.+ 2.*x)*(-1.+ (z*z));
769  outputValues(23,i0, 4) = 2.*(-1.+ 2.*x)*y*z;
770  outputValues(23,i0, 5) = (-1.+ 2.*x)*(-1.+ (y*y));
771  outputValues(23,i0, 6) = 0.;
772  outputValues(23,i0, 7) = 2.*(-1.+ x)*x*z;
773  outputValues(23,i0, 8) = 2.*(-1.+ x)*x*y;
774  outputValues(23,i0, 9) = 0.;
775 
776  outputValues(24,i0, 0) = 0.;
777  outputValues(24,i0, 1) = 2.*y*(-1.+ (z*z));
778  outputValues(24,i0, 2) = 2.*(-1.+ (y*y))*z;
779  outputValues(24,i0, 3) = (1.+ 2.*x)*(-1.+ (z*z));
780  outputValues(24,i0, 4) = 2.*(1.+ 2.*x)*y*z;
781  outputValues(24,i0, 5) = (1.+ 2.*x)*(-1.+ (y*y));
782  outputValues(24,i0, 6) = 0.;
783  outputValues(24,i0, 7) = 2.*x*(1.+ x)*z;
784  outputValues(24,i0, 8) = 2.*x*(1.+ x)*y;
785  outputValues(24,i0, 9) = 0.;
786 
787  outputValues(25,i0, 0) = 0.;
788  outputValues(25,i0, 1) = (-1.+ 2.*y)*(-1.+ (z*z));
789  outputValues(25,i0, 2) = 2.*(-1.+ y)*y*z;
790  outputValues(25,i0, 3) = 2.*x*(-1.+ (z*z));
791  outputValues(25,i0, 4) = 2.*x*(-1.+ 2.*y)*z;
792  outputValues(25,i0, 5) = 2.*x*(-1.+ y)*y;
793  outputValues(25,i0, 6) = 0.;
794  outputValues(25,i0, 7) = 2.*(-1.+ (x*x))*z;
795  outputValues(25,i0, 8) = (-1.+ (x*x))*(-1.+ 2.*y);
796  outputValues(25,i0, 9) = 0.;
797 
798  outputValues(26,i0, 0) = 0.;
799  outputValues(26,i0, 1) = (1.+ 2.*y)*(-1.+ (z*z));
800  outputValues(26,i0, 2) = 2.*y*(1.+ y)*z;
801  outputValues(26,i0, 3) = 2.*x*(-1.+ (z*z));
802  outputValues(26,i0, 4) = 2.*x*(1.+ 2.*y)*z;
803  outputValues(26,i0, 5) = 2.*x*y*(1.+ y);
804  outputValues(26,i0, 6) = 0.;
805  outputValues(26,i0, 7) = 2.*(-1.+ (x*x))*z;
806  outputValues(26,i0, 8) = (-1.+ (x*x))*(1.+ 2.*y);
807  outputValues(26,i0, 9) = 0.;
808  }
809  break;
810 
811  case OPERATOR_D4:
812  {
813  // Non-zero entries have Dk (derivative cardinality) indices {3,4,5,7,8,12}, all other entries are 0.
814  // Intitialize array by zero and then fill only non-zero entries.
815  int DkCardinality = Intrepid::getDkCardinality(operatorType, this -> basisCellTopology_.getDimension() );
816  for(int dofOrd = 0; dofOrd < this -> basisCardinality_; dofOrd++) {
817  for (int i0 = 0; i0 < dim0; i0++) {
818  for(int dkOrd = 0; dkOrd < DkCardinality; dkOrd++){
819  outputValues(dofOrd, i0, dkOrd) = 0.0;
820  }
821  }
822  }
823 
824  for (int i0 = 0; i0 < dim0; i0++) {
825  x = inputPoints(i0,0);
826  y = inputPoints(i0,1);
827  z = inputPoints(i0,2);
828 
829  outputValues(0, i0, 3) = ((-1.+ z)*z)/2.;
830  outputValues(0, i0, 4) = ((-1.+ 2.*y)*(-1.+ 2.*z))/4.;
831  outputValues(0, i0, 5) = ((-1.+ y)*y)/2.;
832  outputValues(0, i0, 7) = ((-1.+ 2.*x)*(-1.+ 2.*z))/4.;
833  outputValues(0, i0, 8) = ((-1.+ 2.*x)*(-1.+ 2.*y))/4.;
834  outputValues(0, i0, 12)= ((-1.+ x)*x)/2.;
835 
836  outputValues(1, i0, 3) = ((-1.+ z)*z)/2.;
837  outputValues(1, i0, 4) = ((-1.+ 2.*y)*(-1.+ 2.*z))/4.;
838  outputValues(1, i0, 5) = ((-1.+ y)*y)/2.;
839  outputValues(1, i0, 7) = ((1. + 2.*x)*(-1.+ 2.*z))/4.;
840  outputValues(1, i0, 8) = ((1. + 2.*x)*(-1.+ 2.*y))/4.;
841  outputValues(1, i0, 12)= (x*(1. + x))/2.;
842 
843  outputValues(2, i0, 3) = ((-1.+ z)*z)/2.;
844  outputValues(2, i0, 4) = ((1. + 2.*y)*(-1.+ 2.*z))/4.;
845  outputValues(2, i0, 5) = (y*(1. + y))/2.;
846  outputValues(2, i0, 7) = ((1. + 2.*x)*(-1.+ 2.*z))/4.;
847  outputValues(2, i0, 8) = ((1. + 2.*x)*(1. + 2.*y))/4.;
848  outputValues(2, i0, 12)= (x*(1. + x))/2.;
849 
850  outputValues(3, i0, 3) = ((-1.+ z)*z)/2.;
851  outputValues(3, i0, 4) = ((1. + 2.*y)*(-1.+ 2.*z))/4.;
852  outputValues(3, i0, 5) = (y*(1. + y))/2.;
853  outputValues(3, i0, 7) = ((-1.+ 2.*x)*(-1.+ 2.*z))/4.;
854  outputValues(3, i0, 8) = ((-1.+ 2.*x)*(1. + 2.*y))/4.;
855  outputValues(3, i0, 12)= ((-1.+ x)*x)/2.;
856 
857  outputValues(4, i0, 3) = (z*(1. + z))/2.;
858  outputValues(4, i0, 4) = ((-1.+ 2.*y)*(1. + 2.*z))/4.;
859  outputValues(4, i0, 5) = ((-1.+ y)*y)/2.;
860  outputValues(4, i0, 7) = ((-1.+ 2.*x)*(1. + 2.*z))/4.;
861  outputValues(4, i0, 8) = ((-1.+ 2.*x)*(-1.+ 2.*y))/4.;
862  outputValues(4, i0, 12)= ((-1.+ x)*x)/2.;
863 
864  outputValues(5, i0, 3) = (z*(1. + z))/2.;
865  outputValues(5, i0, 4) = ((-1.+ 2.*y)*(1. + 2.*z))/4.;
866  outputValues(5, i0, 5) = ((-1.+ y)*y)/2.;
867  outputValues(5, i0, 7) = ((1. + 2.*x)*(1. + 2.*z))/4.;
868  outputValues(5, i0, 8) = ((1. + 2.*x)*(-1.+ 2.*y))/4.;
869  outputValues(5, i0, 12)= (x*(1. + x))/2.;
870 
871  outputValues(6, i0, 3) = (z*(1. + z))/2.;
872  outputValues(6, i0, 4) = ((1. + 2.*y)*(1. + 2.*z))/4.;
873  outputValues(6, i0, 5) = (y*(1. + y))/2.;
874  outputValues(6, i0, 7) = ((1. + 2.*x)*(1. + 2.*z))/4.;
875  outputValues(6, i0, 8) = ((1. + 2.*x)*(1. + 2.*y))/4.;
876  outputValues(6, i0, 12)= (x*(1. + x))/2.;
877 
878  outputValues(7, i0, 3) = (z*(1. + z))/2.;
879  outputValues(7, i0, 4) = ((1. + 2.*y)*(1. + 2.*z))/4.;
880  outputValues(7, i0, 5) = (y*(1. + y))/2.;
881  outputValues(7, i0, 7) = ((-1.+ 2.*x)*(1. + 2.*z))/4.;
882  outputValues(7, i0, 8) = ((-1.+ 2.*x)*(1. + 2.*y))/4.;
883  outputValues(7, i0, 12)= ((-1.+ x)*x)/2.;
884 
885  outputValues(8, i0, 3) = -((-1.+ z)*z);
886  outputValues(8, i0, 4) = -0.5 + y + z - 2.*y*z;
887  outputValues(8, i0, 5) = -((-1.+ y)*y);
888  outputValues(8, i0, 7) = x - 2.*x*z;
889  outputValues(8, i0, 8) = x - 2.*x*y;
890  outputValues(8, i0, 12)= 1. - x*x;
891 
892  outputValues(9, i0, 3) = -((-1.+ z)*z);
893  outputValues(9, i0, 4) = y - 2.*y*z;
894  outputValues(9, i0, 5) = 1 - y*y;
895  outputValues(9, i0, 7) = 0.5 + x - z - 2.*x*z;
896  outputValues(9, i0, 8) = -((1. + 2.*x)*y);
897  outputValues(9, i0, 12)= -(x*(1. + x));
898 
899  outputValues(10,i0, 3) = -((-1.+ z)*z);
900  outputValues(10,i0, 4) = 0.5 + y - z - 2.*y*z;
901  outputValues(10,i0, 5) = -(y*(1. + y));
902  outputValues(10,i0, 7) = x - 2.*x*z;
903  outputValues(10,i0, 8) = -(x*(1. + 2.*y));
904  outputValues(10,i0, 12)= 1. - x*x;
905 
906  outputValues(11,i0, 3) = -((-1.+ z)*z);
907  outputValues(11,i0, 4) = y - 2.*y*z;
908  outputValues(11,i0, 5) = 1. - y*y;
909  outputValues(11,i0, 7) = -0.5 + x + z - 2.*x*z;
910  outputValues(11,i0, 8) = y - 2.*x*y;
911  outputValues(11,i0, 12)= -((-1.+ x)*x);
912 
913  outputValues(12,i0, 3) = 1. - z*z;
914  outputValues(12,i0, 4) = z - 2.*y*z;
915  outputValues(12,i0, 5) = -((-1.+ y)*y);
916  outputValues(12,i0, 7) = z - 2.*x*z;
917  outputValues(12,i0, 8) = -0.5 + x + y - 2.*x*y;
918  outputValues(12,i0, 12)= -((-1.+ x)*x);
919 
920  outputValues(13,i0, 3) = 1. - z*z;
921  outputValues(13,i0, 4) = z - 2.*y*z;
922  outputValues(13,i0, 5) = -((-1.+ y)*y);
923  outputValues(13,i0, 7) = -((1. + 2.*x)*z);
924  outputValues(13,i0, 8) = 0.5 + x - y - 2.*x*y;
925  outputValues(13,i0, 12)= -(x*(1. + x));
926 
927  outputValues(14,i0, 3) = 1. - z*z;
928  outputValues(14,i0, 4) = -((1. + 2.*y)*z);
929  outputValues(14,i0, 5) = -(y*(1. + y));
930  outputValues(14,i0, 7) = -((1. + 2.*x)*z);
931  outputValues(14,i0, 8) = -((1. + 2.*x)*(1. + 2.*y))/2.;
932  outputValues(14,i0, 12)= -(x*(1. + x));
933 
934  outputValues(15,i0, 3) = 1. - z*z;
935  outputValues(15,i0, 4) = -((1. + 2.*y)*z);
936  outputValues(15,i0, 5) = -(y*(1. + y));
937  outputValues(15,i0, 7) = z - 2.*x*z;
938  outputValues(15,i0, 8) = 0.5 + y - x*(1. + 2.*y);
939  outputValues(15,i0, 12)= -((-1.+ x)*x);
940 
941  outputValues(16,i0, 3) = -(z*(1. + z));
942  outputValues(16,i0, 4) = 0.5 + z - y*(1. + 2.*z);
943  outputValues(16,i0, 5) = -((-1.+ y)*y);
944  outputValues(16,i0, 7) = -(x*(1. + 2.*z));
945  outputValues(16,i0, 8) = x - 2.*x*y;
946  outputValues(16,i0, 12)= 1. - x*x;
947 
948  outputValues(17,i0, 3) = -(z*(1. + z));
949  outputValues(17,i0, 4) = -(y*(1. + 2.*z));
950  outputValues(17,i0, 5) = 1. - y*y;
951  outputValues(17,i0, 7) = -((1. + 2.*x)*(1. + 2.*z))/2.;
952  outputValues(17,i0, 8) = -((1. + 2.*x)*y);
953  outputValues(17,i0, 12)= -(x*(1. + x));
954 
955  outputValues(18,i0, 3) = -(z*(1. + z));
956  outputValues(18,i0, 4) = -((1. + 2.*y)*(1. + 2.*z))/2.;
957  outputValues(18,i0, 5) = -(y*(1. + y));
958  outputValues(18,i0, 7) = -(x*(1. + 2.*z));
959  outputValues(18,i0, 8) = -(x*(1. + 2.*y));
960  outputValues(18,i0, 12)= 1. - x*x;
961 
962  outputValues(19,i0, 3) = -(z*(1. + z));
963  outputValues(19,i0, 4) = -(y*(1. + 2.*z));
964  outputValues(19,i0, 5) = 1. - y*y;
965  outputValues(19,i0, 7) = 0.5 + z - x*(1. + 2.*z);
966  outputValues(19,i0, 8) = y - 2.*x*y;
967  outputValues(19,i0, 12)= -((-1.+ x)*x);
968 
969  outputValues(20,i0, 3) = 4. - 4.*z*z;
970  outputValues(20,i0, 4) = -8.*y*z;
971  outputValues(20,i0, 5) = 4. - 4.*y*y;
972  outputValues(20,i0, 7) = -8.*x*z;
973  outputValues(20,i0, 8) = -8.*x*y;
974  outputValues(20,i0, 12)= 4. - 4.*x*x;
975 
976  outputValues(21,i0, 3) = 2.*(-1.+ z)*z;
977  outputValues(21,i0, 4) = 2.*y*(-1.+ 2.*z);
978  outputValues(21,i0, 5) = 2.*(-1.+ y*y);
979  outputValues(21,i0, 7) = 2.*x*(-1.+ 2.*z);
980  outputValues(21,i0, 8) = 4.*x*y;
981  outputValues(21,i0, 12)= 2.*(-1.+ x*x);
982 
983  outputValues(22,i0, 3) = 2.*z*(1. + z);
984  outputValues(22,i0, 4) = 2.*(y + 2.*y*z);
985  outputValues(22,i0, 5) = 2.*(-1.+ y*y);
986  outputValues(22,i0, 7) = 2.*(x + 2.*x*z);
987  outputValues(22,i0, 8) = 4.*x*y;
988  outputValues(22,i0, 12)= 2.*(-1.+ x*x);
989 
990  outputValues(23,i0, 3) = 2.*(-1.+ z*z);
991  outputValues(23,i0, 4) = 4.*y*z;
992  outputValues(23,i0, 5) = 2.*(-1.+ y*y);
993  outputValues(23,i0, 7) = 2.*(-1.+ 2.*x)*z;
994  outputValues(23,i0, 8) = 2.*(-1.+ 2.*x)*y;
995  outputValues(23,i0, 12)= 2.*(-1.+ x)*x;
996 
997  outputValues(24,i0, 3) = 2.*(-1.+ z*z);
998  outputValues(24,i0, 4) = 4.*y*z;
999  outputValues(24,i0, 5) = 2.*(-1.+ y*y);
1000  outputValues(24,i0, 7) = 2.*(z + 2.*x*z);
1001  outputValues(24,i0, 8) = 2.*(y + 2.*x*y);
1002  outputValues(24,i0, 12)= 2.*x*(1. + x);
1003 
1004  outputValues(25,i0, 3) = 2.*(-1.+ z*z);
1005  outputValues(25,i0, 4) = 2.*(-1.+ 2.*y)*z;
1006  outputValues(25,i0, 5) = 2.*(-1.+ y)*y;
1007  outputValues(25,i0, 7) = 4.*x*z;
1008  outputValues(25,i0, 8) = 2.*x*(-1.+ 2.*y);
1009  outputValues(25,i0, 12)= 2.*(-1.+ x*x);
1010 
1011  outputValues(26,i0, 3) = 2.*(-1.+ z*z);
1012  outputValues(26,i0, 4) = 2.*(z + 2.*y*z);
1013  outputValues(26,i0, 5) = 2.*y*(1. + y);
1014  outputValues(26,i0, 7) = 4.*x*z;
1015  outputValues(26,i0, 8) = 2.*(x + 2.*x*y);
1016  outputValues(26,i0, 12)= 2.*(-1.+ x*x);
1017  }
1018  }
1019  break;
1020 
1021  case OPERATOR_D5:
1022  case OPERATOR_D6:
1023  TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument,
1024  ">>> ERROR (Basis_HGRAD_HEX_C2_FEM): operator not supported");
1025  break;
1026 
1027  case OPERATOR_D7:
1028  case OPERATOR_D8:
1029  case OPERATOR_D9:
1030  case OPERATOR_D10:
1031  {
1032  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, DkCardinality)
1033  int DkCardinality = Intrepid::getDkCardinality(operatorType,
1034  this -> basisCellTopology_.getDimension() );
1035  for(int dofOrd = 0; dofOrd < this -> basisCardinality_; dofOrd++) {
1036  for (int i0 = 0; i0 < dim0; i0++) {
1037  for(int dkOrd = 0; dkOrd < DkCardinality; dkOrd++){
1038  outputValues(dofOrd, i0, dkOrd) = 0.0;
1039  }
1040  }
1041  }
1042  }
1043  break;
1044 
1045  default:
1046  TEUCHOS_TEST_FOR_EXCEPTION( !( Intrepid::isValidOperator(operatorType) ), std::invalid_argument,
1047  ">>> ERROR (Basis_HGRAD_HEX_C2_FEM): Invalid operator type");
1048  }
1049 }
1050 
1051 
1052 
1053 template<class Scalar, class ArrayScalar>
1054 void Basis_HGRAD_HEX_C2_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar& outputValues,
1055  const ArrayScalar & inputPoints,
1056  const ArrayScalar & cellVertices,
1057  const EOperator operatorType) const {
1058  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
1059  ">>> ERROR (Basis_HGRAD_HEX_C2_FEM): FEM Basis calling an FVD member function");
1060  }
1061 
1062 template<class Scalar, class ArrayScalar>
1063 void Basis_HGRAD_HEX_C2_FEM<Scalar, ArrayScalar>::getDofCoords(ArrayScalar & DofCoords) const {
1064 #ifdef HAVE_INTREPID_DEBUG
1065  // Verify rank of output array.
1066  TEUCHOS_TEST_FOR_EXCEPTION( !(DofCoords.rank() == 2), std::invalid_argument,
1067  ">>> ERROR: (Intrepid::Basis_HGRAD_HEX_C2_FEM::getDofCoords) rank = 2 required for DofCoords array");
1068  // Verify 0th dimension of output array.
1069  TEUCHOS_TEST_FOR_EXCEPTION( !( DofCoords.dimension(0) == this -> basisCardinality_ ), std::invalid_argument,
1070  ">>> ERROR: (Intrepid::Basis_HGRAD_HEX_C2_FEM::getDofCoords) mismatch in number of DoF and 0th dimension of DofCoords array");
1071  // Verify 1st dimension of output array.
1072  TEUCHOS_TEST_FOR_EXCEPTION( !( DofCoords.dimension(1) == (int)(this -> basisCellTopology_.getDimension()) ), std::invalid_argument,
1073  ">>> ERROR: (Intrepid::Basis_HGRAD_HEX_C2_FEM::getDofCoords) incorrect reference cell (1st) dimension in DofCoords array");
1074 #endif
1075 
1076  DofCoords(0,0) = -1.0; DofCoords(0,1) = -1.0; DofCoords(0,2) = -1.0;
1077  DofCoords(1,0) = 1.0; DofCoords(1,1) = -1.0; DofCoords(1,2) = -1.0;
1078  DofCoords(2,0) = 1.0; DofCoords(2,1) = 1.0; DofCoords(2,2) = -1.0;
1079  DofCoords(3,0) = -1.0; DofCoords(3,1) = 1.0; DofCoords(3,2) = -1.0;
1080  DofCoords(4,0) = -1.0; DofCoords(4,1) = -1.0; DofCoords(4,2) = 1.0;
1081  DofCoords(5,0) = 1.0; DofCoords(5,1) = -1.0; DofCoords(5,2) = 1.0;
1082  DofCoords(6,0) = 1.0; DofCoords(6,1) = 1.0; DofCoords(6,2) = 1.0;
1083  DofCoords(7,0) = -1.0; DofCoords(7,1) = 1.0; DofCoords(7,2) = 1.0;
1084 
1085  DofCoords(8,0) = 0.0; DofCoords(8,1) = -1.0; DofCoords(8,2) = -1.0;
1086  DofCoords(9,0) = 1.0; DofCoords(9,1) = 0.0; DofCoords(9,2) = -1.0;
1087  DofCoords(10,0) = 0.0; DofCoords(10,1) = 1.0; DofCoords(10,2) = -1.0;
1088  DofCoords(11,0) = -1.0; DofCoords(11,1) = 0.0; DofCoords(11,2) = -1.0;
1089  DofCoords(12,0) = -1.0; DofCoords(12,1) = -1.0; DofCoords(12,2) = 0.0;
1090  DofCoords(13,0) = 1.0; DofCoords(13,1) = -1.0; DofCoords(13,2) = 0.0;
1091  DofCoords(14,0) = 1.0; DofCoords(14,1) = 1.0; DofCoords(14,2) = 0.0;
1092  DofCoords(15,0) = -1.0; DofCoords(15,1) = 1.0; DofCoords(15,2) = 0.0;
1093  DofCoords(16,0) = 0.0; DofCoords(16,1) = -1.0; DofCoords(16,2) = 1.0;
1094  DofCoords(17,0) = 1.0; DofCoords(17,1) = 0.0; DofCoords(17,2) = 1.0;
1095  DofCoords(18,0) = 0.0; DofCoords(18,1) = 1.0; DofCoords(18,2) = 1.0;
1096  DofCoords(19,0) = -1.0; DofCoords(19,1) = 0.0; DofCoords(19,2) = 1.0;
1097 
1098  DofCoords(20,0) = 0.0; DofCoords(20,1) = 0.0; DofCoords(20,2) = 0.0;
1099 
1100  DofCoords(21,0) = 0.0; DofCoords(21,1) = 0.0; DofCoords(21,2) = -1.0;
1101  DofCoords(22,0) = 0.0; DofCoords(22,1) = 0.0; DofCoords(22,2) = 1.0;
1102  DofCoords(23,0) = -1.0; DofCoords(23,1) = 0.0; DofCoords(23,2) = 0.0;
1103  DofCoords(24,0) = 1.0; DofCoords(24,1) = 0.0; DofCoords(24,2) = 0.0;
1104  DofCoords(25,0) = 0.0; DofCoords(25,1) = -1.0; DofCoords(25,2) = 0.0;
1105  DofCoords(26,0) = 0.0; DofCoords(26,1) = 1.0; DofCoords(26,2) = 0.0;
1106 }
1107 
1108 }// namespace Intrepid
1109 #endif
Intrepid::Basis_HGRAD_HEX_C2_FEM::getValues
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Hexahedron cell.
Definition: Intrepid_HGRAD_HEX_C2_FEMDef.hpp:120
Intrepid::Basis_HGRAD_HEX_C2_FEM::initializeTags
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
Definition: Intrepid_HGRAD_HEX_C2_FEMDef.hpp:68
Intrepid::Basis_HGRAD_HEX_C2_FEM::Basis_HGRAD_HEX_C2_FEM
Basis_HGRAD_HEX_C2_FEM()
Constructor.
Definition: Intrepid_HGRAD_HEX_C2_FEMDef.hpp:55
Intrepid::Basis_HGRAD_HEX_C2_FEM::getDofCoords
void getDofCoords(ArrayScalar &DofCoords) const
Returns spatial locations (coordinates) of degrees of freedom on a reference Quadrilateral.
Definition: Intrepid_HGRAD_HEX_C2_FEMDef.hpp:1063
Generated by   doxygen 1.8.5