Intrepid
Intrepid_HGRAD_PYR_I2_FEMDef.hpp
Go to the documentation of this file.
1 #ifndef INTREPID_HGRAD_PYR_I2_FEMDEF_HPP
2 #define INTREPID_HGRAD_PYR_I2_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  template<class Scalar, class ArrayScalar>
55  {
56  this -> basisCardinality_ = 13;
57  this -> basisDegree_ = 2;
58  this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Pyramid<5> >() );
59  this -> basisType_ = BASIS_FEM_DEFAULT;
60  this -> basisCoordinates_ = COORDINATES_CARTESIAN;
61  this -> basisTagsAreSet_ = false;
62  }
63 
64 
65 template<class Scalar, class ArrayScalar>
67 
68  // Basis-dependent intializations
69  int tagSize = 4; // size of DoF tag
70  int posScDim = 0; // position in the tag, counting from 0, of the subcell dim
71  int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
72  int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
73 
74  // An array with local DoF tags assigned to basis functions, in the order of their local enumeration
75  int tags[] = { 0, 0, 0, 1,
76  0, 1, 0, 1,
77  0, 2, 0, 1,
78  0, 3, 0, 1,
79  0, 4, 0, 1,
80  1, 0, 0, 1,
81  1, 1, 0, 1,
82  1, 2, 0, 1,
83  1, 3, 0, 1,
84  1, 4, 0, 1,
85  1, 5, 0, 1,
86  1, 6, 0, 1,
87  1, 7, 0, 1,
88  };
89 
90  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
91  Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
92  this -> ordinalToTag_,
93  tags,
94  this -> basisCardinality_,
95  tagSize,
96  posScDim,
97  posScOrd,
98  posDfOrd);
99 }
100 
101 
102 
103 template<class Scalar, class ArrayScalar>
105  const ArrayScalar & inputPoints,
106  const EOperator operatorType) const {
107 
108  // Verify arguments
109 #ifdef HAVE_INTREPID_DEBUG
110  Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
111  inputPoints,
112  operatorType,
113  this -> getBaseCellTopology(),
114  this -> getCardinality() );
115 #endif
116 
117  // Number of evaluation points = dim 0 of inputPoints
118  int dim0 = inputPoints.dimension(0);
119 
120  // Temporaries: (x,y,z) coordinates of the evaluation point
121  Scalar x = 0.0;
122  Scalar y = 0.0;
123  Scalar z = 0.0;
124  const Scalar eps = std::numeric_limits<Scalar>::epsilon( );
125 
126  switch (operatorType) {
127 
128  case OPERATOR_VALUE:
129  for (int i0 = 0; i0 < dim0; i0++) {
130  x = inputPoints(i0, 0);
131  y = inputPoints(i0, 1);
132  z = inputPoints(i0, 2);
133 
134  //be sure that the basis functions are defined when z is very close to 1.
135 
136  if(fabs(z-1.0) < eps) {
137  if(z <= 1.0) z = 1.0-eps;
138  else z = 1.0+eps;
139  //z = 1.0-eps;
140  }
141 
142  Scalar w = 1.0/(1.0 - z);
143 
144  // outputValues is a rank-2 array with dimensions (basisCardinality_, dim0)
145  outputValues(0, i0) = 0.25 * (-x - y - 1.0)*((1.0-x)*(1.0-y) - z + x*y*z*w);
146  outputValues(1, i0) = 0.25 * ( x - y - 1.0)*((1.0+x)*(1.0-y) - z - x*y*z*w);
147  outputValues(2, i0) = 0.25 * ( x + y - 1.0)*((1.0+x)*(1.0+y) - z + x*y*z*w);
148  outputValues(3, i0) = 0.25 * (-x + y - 1.0)*((1.0-x)*(1.0+y) - z - x*y*z*w);
149 
150  outputValues(4, i0) = z * (2.0*z - 1.0);
151 
152  outputValues(5, i0) = 0.5 * (1.0 + x - z)*(1.0 - x - z)*(1.0 - y - z)*w;
153  outputValues(6, i0) = 0.5 * (1.0 + y - z)*(1.0 - y - z)*(1.0 + x - z)*w;
154  outputValues(7, i0) = 0.5 * (1.0 + x - z)*(1.0 - x - z)*(1.0 + y - z)*w;
155  outputValues(8, i0) = 0.5 * (1.0 + y - z)*(1.0 - y - z)*(1.0 - x - z)*w;
156 
157  outputValues(9, i0) = z*(1.0 - x - z)*(1.0 - y - z)*w;
158  outputValues(10,i0) = z*(1.0 + x - z)*(1.0 - y - z)*w;
159  outputValues(11,i0) = z*(1.0 + x - z)*(1.0 + y - z)*w;
160  outputValues(12,i0) = z*(1.0 - x - z)*(1.0 + y - z)*w;
161  }
162  break;
163 
164  case OPERATOR_GRAD:
165  case OPERATOR_D1:
166  for (int i0 = 0; i0 < dim0; i0++) {
167  x = inputPoints(i0,0);
168  y = inputPoints(i0,1);
169  z = inputPoints(i0,2);
170 
171  //be sure that the basis functions are defined when z is very close to 1.
172 
173  if(fabs(z-1.0) < eps) {
174  if(z <= 1.0) z = 1.0-eps;
175  else z = 1.0+eps;
176  }
177 
178  Scalar w = 1.0/(1.0 - z);
179 
180  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
181  outputValues(0, i0, 0) = 0.25*(-1.0-x-y)*(-1.0+y + y*z*w) - 0.25*((1.0-x)*(1.0-y)-z + x*y*z*w);
182  outputValues(0, i0, 1) = 0.25*(-1.0-x-y)*(-1.0+x + x*z*w) - 0.25*((1.0-x)*(1.0-y)-z + x*y*z*w);
183  outputValues(0, i0, 2) = 0.25*(-1.0-x-y)*(-1.0 + x*y*w + x*y*z*w*w);
184 
185  outputValues(1, i0, 0) = 0.25*(-1.0+x-y)*( 1.0-y - y*z*w) + 0.25*((1.0+x)*(1.0-y)-z - x*y*z*w);
186  outputValues(1, i0, 1) = 0.25*(-1.0+x-y)*(-1.0-x - x*z*w) - 0.25*((1.0+x)*(1.0-y)-z - x*y*z*w);
187  outputValues(1, i0, 2) = 0.25*(-1.0+x-y)*(-1.0 - x*y*w - x*y*z*w*w);
188 
189  outputValues(2, i0, 0) = 0.25*(-1.0+x+y)*(1.0+y + y*z*w) + 0.25*((1.0+x)*(1.0+y)-z + x*y*z*w);
190  outputValues(2, i0, 1) = 0.25*(-1.0+x+y)*(1.0+x + x*z*w) + 0.25*((1.0+x)*(1.0+y)-z + x*y*z*w);
191  outputValues(2, i0, 2) = 0.25*(-1.0+x+y)*(-1.0 + x*y*w + x*y*z*w*w);
192 
193  outputValues(3, i0, 0) = 0.25*(-1.0-x+y)*(-1.0-y - y*z*w) - 0.25*((1.0-x)*(1.0+y)-z - x*y*z*w);
194  outputValues(3, i0, 1) = 0.25*(-1.0-x+y)*( 1.0-x - x*z*w) + 0.25*((1.0-x)*(1.0+y)-z - x*y*z*w);
195  outputValues(3, i0, 2) = 0.25*(-1.0-x+y)*(-1.0 - x*y*w - x*y*z*w*w);
196 
197  outputValues(4, i0, 0) = 0.0;
198  outputValues(4, i0, 1) = 0.0;
199  outputValues(4, i0, 2) = -1.0 + 4.0*z;
200 
201  outputValues(5, i0, 0) = -x*w*(1.0-y-z);
202  outputValues(5, i0, 1) = -0.5*(1.0-x-z)*(1.0+x-z)*w;
203  outputValues(5, i0, 2) = 0.5*y*x*x*w*w + 0.5*y - 1.0+z;
204 
205  outputValues(6, i0, 0) = 0.5*(1.0-y-z)*(1.0+y-z)*w;
206  outputValues(6, i0, 1) = -y*w*(1.0+x-z);
207  outputValues(6, i0, 2) = -0.5*x*y*y*w*w - 0.5*x - 1.0+z;
208 
209  outputValues(7, i0, 0) = -x*w*(1.0+y-z);
210  outputValues(7, i0, 1) = 0.5*(1.0-x-z)*(1.0+x-z)*w;
211  outputValues(7, i0, 2) = -0.5*y*x*x*w*w - 0.5*y - 1.0+z;
212 
213  outputValues(8, i0, 0) = -0.5*(1.0-y-z)*(1.0+y-z)*w;
214  outputValues(8, i0, 1) = -y*w*(1.0-x-z);
215  outputValues(8, i0, 2) = 0.5*x*y*y*w*w + 0.5*x - 1.0+z;
216 
217  outputValues(9, i0, 0) = -(1.0-y-z)*z*w;
218  outputValues(9, i0, 1) = -(1.0-x-z)*z*w;
219  outputValues(9, i0, 2) = x*y*w*w + 1.0 - x - y - 2.0*z;
220 
221  outputValues(10,i0, 0) = (1.0-y-z)*z*w;
222  outputValues(10,i0, 1) = -(1.0+x-z)*z*w;
223  outputValues(10,i0, 2) = -x*y*w*w + 1.0 + x - y - 2.0*z;
224 
225  outputValues(11,i0, 0) = (1.0+y-z)*z*w;
226  outputValues(11,i0, 1) = (1.0+x-z)*z*w;
227  outputValues(11,i0, 2) = x*y*w*w + 1.0 + x + y - 2.0*z;
228 
229  outputValues(12,i0, 0) = -(1.0+y-z)*z*w;
230  outputValues(12,i0, 1) = (1.0-x-z)*z*w;
231  outputValues(12,i0, 2) = -x*y*w*w + 1.0 - x + y - 2.0*z;
232 
233  }
234  break;
235 
236  case OPERATOR_CURL:
237  TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_CURL), std::invalid_argument,
238  ">>> ERROR (Basis_HGRAD_PYR_I2_FEM): CURL is invalid operator for rank-0 (scalar) functions in 3D");
239  break;
240 
241  case OPERATOR_DIV:
242  TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
243  ">>> ERROR (Basis_HGRAD_PYR_I2_FEM): DIV is invalid operator for rank-0 (scalar) functions in 3D");
244  break;
245 
246  case OPERATOR_D2:
247  for (int i0 = 0; i0 < dim0; i0++) {
248  x = inputPoints(i0,0);
249  y = inputPoints(i0,1);
250  z = inputPoints(i0,2);
251 
252  if(fabs(z-1.0) < eps) {
253  if(z <= 1.0) z = 1.0-eps;
254  else z = 1.0+eps;
255  }
256  Scalar w = 1.0/(1.0 - z);
257 
258  outputValues(0, i0, 0) = -0.5*(-1.0+y+y*z*w);
259  outputValues(0, i0, 1) = -(-0.25 + 0.5*x + 0.5*y + 0.5*z)*w;
260  outputValues(0, i0, 2) = 0.25 + (-0.25*y-0.5*x*y-0.25*y*y)*w*w;
261  outputValues(0, i0, 3) = -0.5*(-1.0+x+x*z*w);
262  outputValues(0, i0, 4) = 0.25 + (-0.25*x*x-0.25*x-0.5*x*y)*w*w;
263  outputValues(0, i0, 5) = 0.5*x*y*(-1.0-x-y)*w*w*w;
264 
265  outputValues(1, i0, 0) = 0.5*(1.0-y-y*z*w);
266  outputValues(1, i0, 1) =-(0.25 + 0.5*x - 0.5*y - 0.5*z)*w;
267  outputValues(1, i0, 2) = -0.25 + (0.25*y-0.5*x*y+0.25*y*y)*w*w;
268  outputValues(1, i0, 3) = -0.5*(-1.0-x-x*z*w);
269  outputValues(1, i0, 4) = 0.25 + (-0.25*x*x + 0.25*x + 0.5*x*y)*w*w;
270  outputValues(1, i0, 5) = -0.5*x*y*(-1.0+x-y)*w*w*w;
271 
272  outputValues(2, i0, 0) = 0.5*(1.0+y+y*z*w);
273  outputValues(2, i0, 1) =-(-0.25 - 0.5*x - 0.5*y + 0.5*z)*w;
274  outputValues(2, i0, 2) = -0.25 + (-0.25*y+0.5*x*y+0.25*y*y)*w*w;
275  outputValues(2, i0, 3) = 0.5*(1.0+x+x*z*w);
276  outputValues(2, i0, 4) = -0.25 + (0.25*x*x -0.25*x + 0.5*x*y)*w*w;
277  outputValues(2, i0, 5) = 0.5*x*y*(-1.0+x+y)*w*w*w;
278 
279  outputValues(3, i0, 0) = -0.5*(-1.0-y-y*z*w);
280  outputValues(3, i0, 1) =-(0.25 - 0.5*x + 0.5*y - 0.5*z)*w;
281  outputValues(3, i0, 2) = 0.25 + (0.25*y+0.5*x*y-0.25*y*y)*w*w;
282  outputValues(3, i0, 3) = 0.5*(1.0-x-x*z*w);
283  outputValues(3, i0, 4) = -0.25 + (0.25*x + 0.25*x*x - 0.5*x*y)*w*w;
284  outputValues(3, i0, 5) = -0.5*x*y*(-1.0-x+y)*w*w*w;
285 
286  outputValues(4, i0, 0) = 0.0;
287  outputValues(4, i0, 1) = 0.0;
288  outputValues(4, i0, 2) = 0.0;
289  outputValues(4, i0, 3) = 0.0;
290  outputValues(4, i0, 4) = 0.0;
291  outputValues(4, i0, 5) = 4.0;
292 
293  outputValues(5, i0, 0) = -(1.0-y-z)*w;
294  outputValues(5, i0, 1) = x*w;
295  outputValues(5, i0, 2) = x*y*w*w;
296  outputValues(5, i0, 3) = 0.0;
297  outputValues(5, i0, 4) = 0.5*x*x*w*w + 0.5;
298  outputValues(5, i0, 5) = x*x*y*w*w*w + 1.0;
299 
300  outputValues(6, i0, 0) = 0.0;
301  outputValues(6, i0, 1) = -y*w;
302  outputValues(6, i0, 2) = -0.5*y*y*w*w - 0.5;
303  outputValues(6, i0, 3) =-(1.0+x-z)*w;
304  outputValues(6, i0, 4) = -x*y*w*w;
305  outputValues(6, i0, 5) = -x*y*y*w*w*w + 1.0;
306 
307  outputValues(7, i0, 0) = -(1.0+y-z)*w;
308  outputValues(7, i0, 1) = -x*w;
309  outputValues(7, i0, 2) = -x*y*w*w;
310  outputValues(7, i0, 3) = 0.0;
311  outputValues(7, i0, 4) = -0.5*x*x*w*w - 0.5;
312  outputValues(7, i0, 5) = -x*x*y*w*w*w + 1.0;
313 
314  outputValues(8, i0, 0) = 0.0;
315  outputValues(8, i0, 1) = y*w;
316  outputValues(8, i0, 2) = 0.5*y*y*w*w + 0.5;
317  outputValues(8, i0, 3) = -(1.0-x-z)*w;
318  outputValues(8, i0, 4) = x*y*w*w;
319  outputValues(8, i0, 5) = x*y*y*w*w*w + 1.0;
320 
321  outputValues(9, i0, 0) = 0.0;
322  outputValues(9, i0, 1) = z*w;
323  outputValues(9, i0, 2) = y*w*w - 1.0;
324  outputValues(9, i0, 3) = 0.0;
325  outputValues(9, i0, 4) = x*w*w - 1.0;
326  outputValues(9, i0, 5) = 2.0*x*y*w*w*w - 2.0;
327 
328  outputValues(10,i0, 0) = 0.0;
329  outputValues(10,i0, 1) = -z*w;
330  outputValues(10,i0, 2) = -y*w*w + 1.0;
331  outputValues(10,i0, 3) = 0.0;
332  outputValues(10,i0, 4) = -x*w*w - 1.0;
333  outputValues(10,i0, 5) = -2.0*x*y*w*w*w - 2.0;
334 
335  outputValues(11,i0, 0) = 0.0;
336  outputValues(11,i0, 1) = z*w;
337  outputValues(11,i0, 2) = y*w*w + 1.0;
338  outputValues(11,i0, 3) = 0.0;
339  outputValues(11,i0, 4) = x*w*w + 1.0;
340  outputValues(11,i0, 5) = 2.0*x*y*w*w*w - 2.0;
341 
342  outputValues(12,i0, 0) = 0.0;
343  outputValues(12,i0, 1) = -z*w;
344  outputValues(12,i0, 2) = -y*w*w - 1.0;
345  outputValues(12,i0, 3) = 0.0;
346  outputValues(12,i0, 4) = -x*w*w + 1.0;
347  outputValues(12,i0, 5) = -2.0*x*y*w*w*w - 2.0;
348 
349  }
350  break;
351 
352  case OPERATOR_D3:
353  for (int i0 = 0; i0 < dim0; i0++) {
354  x = inputPoints(i0,0);
355  y = inputPoints(i0,1);
356  z = inputPoints(i0,2);
357 
358  if(fabs(z-1.0) < eps) {
359  if(z <= 1.0) z = 1.0-eps;
360  else z = 1.0+eps;
361  }
362 
363  Scalar w = 1.0/(1.0 - z);
364 
365  outputValues(0, i0, 0) = 0.0;
366  outputValues(0, i0, 1) = -0.5*w;
367  outputValues(0, i0, 2) = -0.5*y*w*w;
368  outputValues(0, i0, 3) = -0.5*w;
369  outputValues(0, i0, 4) =(-0.25 - 0.5*x - 0.5*y)*w*w;
370  outputValues(0, i0, 5) =-(0.5*y + x*y + 0.5*y*y)*w*w*w;
371  outputValues(0, i0, 6) = 0.0;
372  outputValues(0, i0, 7) = -0.5*x*w*w;
373  outputValues(0, i0, 8) =-(0.5*x + x*y + 0.5*x*x)*w*w*w;
374  outputValues(0, i0, 9) = 1.5*x*y*(-1.0 - x - y)*w*w*w*w;
375 
376  outputValues(1, i0, 0) = 0.0;
377  outputValues(1, i0, 1) = -0.5*w;
378  outputValues(1, i0, 2) = -0.5*y*w*w;
379  outputValues(1, i0, 3) = 0.5*w;
380  outputValues(1, i0, 4) = (0.25 - 0.5*x - 0.5*y)*w*w;
381  outputValues(1, i0, 5) =-(-0.5*y + x*y - 0.5*y*y)*w*w*w;
382  outputValues(1, i0, 6) = 0.0;
383  outputValues(1, i0, 7) = 0.5*x*w*w;
384  outputValues(1, i0, 8) =-(-0.5*x - x*y + 0.5*x*x)*w*w*w;
385  outputValues(1, i0, 9) = -1.5*x*y*(-1.0 + x - y)*w*w*w*w;
386 
387  outputValues(2, i0, 0) = 0.0;
388  outputValues(2, i0, 1) = 0.5*w;
389  outputValues(2, i0, 2) = 0.5*y*w*w;
390  outputValues(2, i0, 3) = 0.5*w;
391  outputValues(2, i0, 4) = (-0.25 + 0.5*x + 0.5*y)*w*w;
392  outputValues(2, i0, 5) =-(0.5*y - x*y - 0.5*y*y)*w*w*w;
393  outputValues(2, i0, 6) = 0.0;
394  outputValues(2, i0, 7) = 0.5*x*w*w;
395  outputValues(2, i0, 8) =-(0.5*x - x*y - 0.5*x*x)*w*w*w;
396  outputValues(2, i0, 9) = 1.5*x*y*(-1.0 + x + y)*w*w*w*w;
397 
398  outputValues(3, i0, 0) = 0.0;
399  outputValues(3, i0, 1) = 0.5*w;
400  outputValues(3, i0, 2) = 0.5*y*w*w;
401  outputValues(3, i0, 3) = -0.5*w;
402  outputValues(3, i0, 4) = (0.25 + 0.5*x - 0.5*y)*w*w;
403  outputValues(3, i0, 5) =-(-0.5*y - x*y + 0.5*y*y)*w*w*w;
404  outputValues(3, i0, 6) = 0.0;
405  outputValues(3, i0, 7) = -0.5*x*w*w;
406  outputValues(3, i0, 8) =-(-0.5*x + x*y - 0.5*x*x)*w*w*w;
407  outputValues(3, i0, 9) = -1.5*x*y*(-1.0 - x + y)*w*w*w*w;
408 
409  outputValues(4, i0, 0) = 0.0;
410  outputValues(4, i0, 1) = 0.0;
411  outputValues(4, i0, 2) = 0.0;
412  outputValues(4, i0, 3) = 0.0;
413  outputValues(4, i0, 4) = 0.0;
414  outputValues(4, i0, 5) = 0.0;
415  outputValues(4, i0, 6) = 0.0;
416  outputValues(4, i0, 7) = 0.0;
417  outputValues(4, i0, 8) = 0.0;
418  outputValues(4, i0, 9) = 0.0;
419 
420  outputValues(5, i0, 0) = 0.0;
421  outputValues(5, i0, 1) = w;
422  outputValues(5, i0, 2) = y*w*w;
423  outputValues(5, i0, 3) = 0.0;
424  outputValues(5, i0, 4) = x*w*w;
425  outputValues(5, i0, 5) = 2.0*y*x*w*w*w;
426  outputValues(5, i0, 6) = 0.0;
427  outputValues(5, i0, 7) = 0.0;
428  outputValues(5, i0, 8) = x*x*w*w*w;
429  outputValues(5, i0, 9) = 3.0*x*x*y*w*w*w*w;
430 
431  outputValues(6, i0, 0) = 0.0;
432  outputValues(6, i0, 1) = 0.0;
433  outputValues(6, i0, 2) = 0.0;
434  outputValues(6, i0, 3) = -w;
435  outputValues(6, i0, 4) = -y*w*w;
436  outputValues(6, i0, 5) = -y*y*w*w*w;
437  outputValues(6, i0, 6) = 0.0;
438  outputValues(6, i0, 7) = -x*w*w ;
439  outputValues(6, i0, 8) = -2.0*x*y*w*w*w;
440  outputValues(6, i0, 9) = -3.0*x*y*y*w*w*w*w;
441 
442  outputValues(7, i0, 0) = 0.0;
443  outputValues(7, i0, 1) = -w;
444  outputValues(7, i0, 2) = -y*w*w;
445  outputValues(7, i0, 3) = 0.0;
446  outputValues(7, i0, 4) = -x*w*w;
447  outputValues(7, i0, 5) = -2.0*x*y*w*w*w;
448  outputValues(7, i0, 6) = 0.0;
449  outputValues(7, i0, 7) = 0.0;
450  outputValues(7, i0, 8) = -x*x*w*w*w;
451  outputValues(7, i0, 9) = -3.0*x*x*y*w*w*w*w;
452 
453  outputValues(8, i0, 0) = 0.0;
454  outputValues(8, i0, 1) = 0.0;
455  outputValues(8, i0, 2) = 0.0;
456  outputValues(8, i0, 3) = w;
457  outputValues(8, i0, 4) = y*w*w;
458  outputValues(8, i0, 5) = y*y*w*w*w;
459  outputValues(8, i0, 6) = 0.0;
460  outputValues(8, i0, 7) = x*w*w;
461  outputValues(8, i0, 8) = 2.0*x*y*w*w*w;
462  outputValues(8, i0, 9) = 3.0*x*y*y*w*w*w*w ;
463 
464  outputValues(9, i0, 0) = 0.0;
465  outputValues(9, i0, 1) = 0.0;
466  outputValues(9, i0, 2) = 0.0;
467  outputValues(9, i0, 3) = 0.0;
468  outputValues(9, i0, 4) = w*w;
469  outputValues(9, i0, 5) = 2.0*y*w*w*w;
470  outputValues(9, i0, 6) = 0.0;
471  outputValues(9, i0, 7) = 0.0;
472  outputValues(9, i0, 8) = 2.0*x*w*w*w;
473  outputValues(9, i0, 9) = 6.0*x*y*w*w*w*w;
474 
475  outputValues(10,i0, 0) = 0.0;
476  outputValues(10,i0, 1) = 0.0;
477  outputValues(10,i0, 2) = 0.0;
478  outputValues(10,i0, 3) = 0.0;
479  outputValues(10,i0, 4) = -w*w;
480  outputValues(10,i0, 5) = -2.0*y*w*w*w;
481  outputValues(10,i0, 6) = 0.0;
482  outputValues(10,i0, 7) = 0.0;
483  outputValues(10,i0, 8) = -2.0*x*w*w*w;
484  outputValues(10,i0, 9) = -6.0*x*y*w*w*w*w;
485 
486  outputValues(11,i0, 0) = 0.0;
487  outputValues(11,i0, 1) = 0.0;
488  outputValues(11,i0, 2) = 0.0;
489  outputValues(11,i0, 3) = 0.0;
490  outputValues(11,i0, 4) = w*w;
491  outputValues(11,i0, 5) = 2.0*y*w*w*w;
492  outputValues(11,i0, 6) = 0.0;
493  outputValues(11,i0, 7) = 0.0;
494  outputValues(11,i0, 8) = 2.0*x*w*w*w;
495  outputValues(11,i0, 9) = 6.0*x*y*w*w*w*w;
496 
497  outputValues(12,i0, 0) = 0.0;
498  outputValues(12,i0, 1) = 0.0;
499  outputValues(12,i0, 2) = 0.0;
500  outputValues(12,i0, 3) = 0.0;
501  outputValues(12,i0, 4) = -w*w;
502  outputValues(12,i0, 5) = -2.0*y*w*w*w;
503  outputValues(12,i0, 6) = 0.0;
504  outputValues(12,i0, 7) = 0.0;
505  outputValues(12,i0, 8) = -2.0*x*w*w*w;
506  outputValues(12,i0, 9) = -6.0*x*y*w*w*w*w;
507 
508  }
509  break;
510 
511  case OPERATOR_D4:
512  case OPERATOR_D5:
513  case OPERATOR_D6:
514  case OPERATOR_D7:
515  case OPERATOR_D8:
516  case OPERATOR_D9:
517  case OPERATOR_D10:
518  {
519  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, DkCardinality)
520  int DkCardinality = Intrepid::getDkCardinality(operatorType,
521  this -> basisCellTopology_.getDimension() );
522  for(int dofOrd = 0; dofOrd < this -> basisCardinality_; dofOrd++) {
523  for (int i0 = 0; i0 < dim0; i0++) {
524  for(int dkOrd = 0; dkOrd < DkCardinality; dkOrd++){
525  outputValues(dofOrd, i0, dkOrd) = 0.0;
526  }
527  }
528  }
529  }
530  break;
531 
532  default:
533  TEUCHOS_TEST_FOR_EXCEPTION( !( Intrepid::isValidOperator(operatorType) ), std::invalid_argument,
534  ">>> ERROR (Basis_HGRAD_PYR_I2_FEM): Invalid operator type");
535  }
536 }
537 
538 
539 
540 template<class Scalar, class ArrayScalar>
542  const ArrayScalar & inputPoints,
543  const ArrayScalar & cellVertices,
544  const EOperator operatorType) const {
545  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
546  ">>> ERROR (Basis_HGRAD_PYR_I2_FEM): FEM Basis calling an FVD member function");
547 }
548 }// namespace Intrepid
549 #endif
void setOrdinalTagData(std::vector< std::vector< std::vector< int > > > &tagToOrdinal, std::vector< std::vector< int > > &ordinalToTag, const int *tags, const int basisCard, const int tagSize, const int posScDim, const int posScOrd, const int posDfOrd)
Fills ordinalToTag_ and tagToOrdinal_ by basis-specific tag data.
EOperator
Enumeration of primitive operators available in Intrepid. Primitive operators act on reconstructed fu...
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
int getDkCardinality(const EOperator operatorType, const int spaceDim)
Returns cardinality of Dk, i.e., the number of all derivatives of order k.
int isValidOperator(const EOperator operatorType)
Verifies validity of an operator enum.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
FEM basis evaluation on a reference Pyramid cell.