Intrepid
Intrepid_HGRAD_WEDGE_C2_FEMDef.hpp
Go to the documentation of this file.
1 #ifndef INTREPID_HGRAD_WEDGE_C2_FEMDEF_HPP
2 #define INTREPID_HGRAD_WEDGE_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  template<class Scalar, class ArrayScalar>
55  {
56  this -> basisCardinality_ = 18;
57  this -> basisDegree_ = 2;
58  this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Wedge<6> >() );
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  0, 5, 0, 1,
81  1, 0, 0, 1,
82  1, 1, 0, 1,
83  1, 2, 0, 1,
84  1, 6, 0, 1,
85  1, 7, 0, 1,
86  1, 8, 0, 1,
87  1, 3, 0, 1,
88  1, 4, 0, 1,
89  1, 5, 0, 1,
90  2, 0, 0, 1,
91  2, 1, 0, 1,
92  2, 2, 0, 1
93  };
94 
95  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
96  Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
97  this -> ordinalToTag_,
98  tags,
99  this -> basisCardinality_,
100  tagSize,
101  posScDim,
102  posScOrd,
103  posDfOrd);
104 }
105 
106 
107 
108 template<class Scalar, class ArrayScalar>
110  const ArrayScalar & inputPoints,
111  const EOperator operatorType) const {
112 
113  // Verify arguments
114 #ifdef HAVE_INTREPID_DEBUG
115  Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
116  inputPoints,
117  operatorType,
118  this -> getBaseCellTopology(),
119  this -> getCardinality() );
120 #endif
121 
122  // Number of evaluation points = dim 0 of inputPoints
123  int dim0 = inputPoints.dimension(0);
124 
125  // Temporaries: (x,y,z) coordinates of the evaluation point
126  Scalar x = 0.0;
127  Scalar y = 0.0;
128  Scalar z = 0.0;
129 
130  switch (operatorType) {
131 
132  case OPERATOR_VALUE:
133  for (int i0 = 0; i0 < dim0; i0++) {
134  x = inputPoints(i0, 0);
135  y = inputPoints(i0, 1);
136  z = inputPoints(i0, 2);
137 
138  // outputValues is a rank-2 array with dimensions (basisCardinality_, dim0)
139  outputValues(0, i0) = ((-1. + x + y)*(-1. + 2.*x + 2.*y)*(-1. + z)*z)/2.;
140  outputValues(1, i0) = (x*(-1. + 2.*x)*(-1. + z)*z)/2.;
141  outputValues(2, i0) = (y*(-1. + 2.*y)*(-1. + z)*z)/2.;
142  outputValues(3, i0) = ((-1. + x + y)*(-1. + 2.*x + 2.*y)*z*(1. + z))/2.;
143  outputValues(4, i0) = (x*(-1. + 2.*x)*z*(1. + z))/2.;
144  outputValues(5, i0) = (y*(-1. + 2.*y)*z*(1. + z))/2.;
145 
146  outputValues(6, i0) = -2.*x*(-1. + x + y)*(-1. + z)*z;
147  outputValues(7, i0) = 2.*x*y*(-1. + z)*z;
148  outputValues(8, i0) = -2.*y*(-1. + x + y)*(-1. + z)*z;
149  outputValues(9, i0) = -((-1. + x + y)*(-1. + 2.*x + 2.*y)*(-1. + z)*(1. + z));
150  outputValues(10,i0) = -(x*(-1. + 2.*x)*(-1. + z)*(1. + z));
151  outputValues(11,i0) = -(y*(-1. + 2.*y)*(-1. + z)*(1. + z));
152  outputValues(12,i0) = -2.*x*(-1. + x + y)*z*(1. + z);
153  outputValues(13,i0) = 2.*x*y*z*(1. + z);
154  outputValues(14,i0) = -2.*y*(-1. + x + y)*z*(1. + z);
155  outputValues(15,i0) = 4.*x*(-1. + x + y)*(-1. + z)*(1. + z);
156  outputValues(16,i0) = -4.*x*y*(-1. + z)*(1. + z);
157  outputValues(17,i0) = 4.*y*(-1. + x + y)*(-1. + z)*(1. + z);
158  }
159  break;
160 
161  case OPERATOR_GRAD:
162  case OPERATOR_D1:
163  for (int i0 = 0; i0 < dim0; i0++) {
164  x = inputPoints(i0,0);
165  y = inputPoints(i0,1);
166  z = inputPoints(i0,2);
167 
168  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
169  outputValues(0, i0, 0) = ((-3 + 4*x + 4*y)*(-1 + z)*z)/2.;
170  outputValues(0, i0, 1) = ((-3 + 4*x + 4*y)*(-1 + z)*z)/2.;
171  outputValues(0, i0, 2) = ((-1 + x + y)*(-1 + 2*x + 2*y)*(-1 + 2*z))/2.;
172 
173  outputValues(1, i0, 0) = ((-1 + 4*x)*(-1 + z)*z)/2.;
174  outputValues(1, i0, 1) = 0.;
175  outputValues(1, i0, 2) = (x*(-1 + 2*x)*(-1 + 2*z))/2.;
176 
177  outputValues(2, i0, 0) = 0.;
178  outputValues(2, i0, 1) = ((-1 + 4*y)*(-1 + z)*z)/2.;
179  outputValues(2, i0, 2) = (y*(-1 + 2*y)*(-1 + 2*z))/2.;
180 
181  outputValues(3, i0, 0) = ((-3 + 4*x + 4*y)*z*(1 + z))/2.;
182  outputValues(3, i0, 1) = ((-3 + 4*x + 4*y)*z*(1 + z))/2.;
183  outputValues(3, i0, 2) = ((-1 + x + y)*(-1 + 2*x + 2*y)*(1 + 2*z))/2.;
184 
185  outputValues(4, i0, 0) = ((-1 + 4*x)*z*(1 + z))/2.;
186  outputValues(4, i0, 1) = 0.;
187  outputValues(4, i0, 2) = (x*(-1 + 2*x)*(1 + 2*z))/2.;
188 
189  outputValues(5, i0, 0) = 0.;
190  outputValues(5, i0, 1) = ((-1 + 4*y)*z*(1 + z))/2.;
191  outputValues(5, i0, 2) = (y*(-1 + 2*y)*(1 + 2*z))/2.;
192 
193  outputValues(6, i0, 0) = -2*(-1 + 2*x + y)*(-1 + z)*z;
194  outputValues(6, i0, 1) = -2*x*(-1 + z)*z;
195  outputValues(6, i0, 2) = 2*x*(-1 + x + y)*(1 - 2*z);
196 
197  outputValues(7, i0, 0) = 2*y*(-1 + z)*z;
198  outputValues(7, i0, 1) = 2*x*(-1 + z)*z;
199  outputValues(7, i0, 2) = 2*x*y*(-1 + 2*z);
200 
201  outputValues(8, i0, 0) = -2*y*(-1 + z)*z;
202  outputValues(8, i0, 1) = -2*(-1 + x + 2*y)*(-1 + z)*z;
203  outputValues(8, i0, 2) = 2*y*(-1 + x + y)*(1 - 2*z);
204 
205  outputValues(9, i0, 0) = -(-3 + 4*x + 4*y)*(-1 + z*z);
206  outputValues(9, i0, 1) = -(-3 + 4*x + 4*y)*(-1 + z*z);
207  outputValues(9, i0, 2) = -2*(1 + 2*x*x - 3*y + 2*y*y + x*(-3 + 4*y))*z;
208 
209  outputValues(10,i0, 0) = -(-1 + 4*x)*(-1 + z*z);
210  outputValues(10,i0, 1) = 0;
211  outputValues(10,i0, 2) = 2*(1 - 2*x)*x*z;
212 
213  outputValues(11,i0, 0) = 0;
214  outputValues(11,i0, 1) = -(-1 + 4*y)*(-1 + z*z);
215  outputValues(11,i0, 2) = 2*(1 - 2*y)*y*z;
216 
217  outputValues(12,i0, 0) = -2*(-1 + 2*x + y)*z*(1 + z);
218  outputValues(12,i0, 1) = -2*x*z*(1 + z);
219  outputValues(12,i0, 2) = -2*x*(-1 + x + y)*(1 + 2*z);
220 
221  outputValues(13,i0, 0) = 2*y*z*(1 + z);
222  outputValues(13,i0, 1) = 2*x*z*(1 + z);
223  outputValues(13,i0, 2) = 2*x*y*(1 + 2*z);
224 
225  outputValues(14,i0, 0) = -2*y*z*(1 + z);
226  outputValues(14,i0, 1) = -2*(-1 + x + 2*y)*z*(1 + z);
227  outputValues(14,i0, 2) = -2*y*(-1 + x + y)*(1 + 2*z);
228 
229  outputValues(15,i0, 0) = 4*(-1 + 2*x + y)*(-1 + z*z);
230  outputValues(15,i0, 1) = 4*x*(-1 + z)*(1 + z);
231  outputValues(15,i0, 2) = 8*x*(-1 + x + y)*z;
232 
233  outputValues(16,i0, 0) = -4*y*(-1 + z)*(1 + z);
234  outputValues(16,i0, 1) = -4*x*(-1 + z)*(1 + z);
235  outputValues(16,i0, 2) = -8*x*y*z;
236 
237  outputValues(17,i0, 0) = 4*y*(-1 + z)*(1 + z);
238  outputValues(17,i0, 1) = 4*(-1 + x + 2*y)*(-1 + z*z);
239  outputValues(17,i0, 2) = 8*y*(-1 + x + y)*z;
240 
241  }
242  break;
243 
244  case OPERATOR_CURL:
245  TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_CURL), std::invalid_argument,
246  ">>> ERROR (Basis_HGRAD_WEDGE_C2_FEM): CURL is invalid operator for rank-0 (scalar) functions in 3D");
247  break;
248 
249  case OPERATOR_DIV:
250  TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
251  ">>> ERROR (Basis_HGRAD_WEDGE_C2_FEM): DIV is invalid operator for rank-0 (scalar) functions in 3D");
252  break;
253 
254  case OPERATOR_D2:
255  for (int i0 = 0; i0 < dim0; i0++) {
256  x = inputPoints(i0,0);
257  y = inputPoints(i0,1);
258  z = inputPoints(i0,2);
259 
260  outputValues(0, i0, 0) = 2.*(-1. + z)*z;
261  outputValues(0, i0, 1) = 2.*(-1. + z)*z;
262  outputValues(0, i0, 2) = ((-3. + 4.*x + 4.*y)*(-1. + 2.*z))/2.;
263  outputValues(0, i0, 3) = 2.*(-1. + z)*z;
264  outputValues(0, i0, 4) = ((-3. + 4.*x + 4.*y)*(-1. + 2.*z))/2.;
265  outputValues(0, i0, 5) = (-1. + x + y)*(-1. + 2.*x + 2.*y);
266 
267  outputValues(1, i0, 0) = 2.*(-1. + z)*z;
268  outputValues(1, i0, 1) = 0.;
269  outputValues(1, i0, 2) = ((-1. + 4.*x)*(-1. + 2.*z))/2.;
270  outputValues(1, i0, 3) = 0.;
271  outputValues(1, i0, 4) = 0.;
272  outputValues(1, i0, 5) = x*(-1. + 2.*x);
273 
274  outputValues(2, i0, 0) = 0.;
275  outputValues(2, i0, 1) = 0.;
276  outputValues(2, i0, 2) = 0.;
277  outputValues(2, i0, 3) = 2.*(-1. + z)*z;
278  outputValues(2, i0, 4) = ((-1. + 4.*y)*(-1. + 2.*z))/2.;
279  outputValues(2, i0, 5) = y*(-1. + 2.*y);
280 
281  outputValues(3, i0, 0) = 2.*z*(1. + z);
282  outputValues(3, i0, 1) = 2.*z*(1. + z);
283  outputValues(3, i0, 2) = ((-3. + 4.*x + 4.*y)*(1. + 2.*z))/2.;
284  outputValues(3, i0, 3) = 2.*z*(1. + z);
285  outputValues(3, i0, 4) = ((-3. + 4.*x + 4.*y)*(1. + 2.*z))/2.;
286  outputValues(3, i0, 5) = (-1. + x + y)*(-1. + 2.*x + 2.*y);
287 
288  outputValues(4, i0, 0) = 2.*z*(1. + z);
289  outputValues(4, i0, 1) = 0.;
290  outputValues(4, i0, 2) = ((-1. + 4.*x)*(1. + 2.*z))/2.;;
291  outputValues(4, i0, 3) = 0.;
292  outputValues(4, i0, 4) = 0.;
293  outputValues(4, i0, 5) = x*(-1. + 2.*x);
294 
295  outputValues(5, i0, 0) = 0.;
296  outputValues(5, i0, 1) = 0.;
297  outputValues(5, i0, 2) = 0.;
298  outputValues(5, i0, 3) = 2.*z*(1. + z);
299  outputValues(5, i0, 4) = ((-1. + 4.*y)*(1. + 2.*z))/2.;
300  outputValues(5, i0, 5) = y*(-1. + 2.*y);
301 
302  outputValues(6, i0, 0) = -4.*(-1. + z)*z;
303  outputValues(6, i0, 1) = -2.*(-1. + z)*z;
304  outputValues(6, i0, 2) = -2.*(-1. + 2.*x + y)*(-1. + 2.*z);
305  outputValues(6, i0, 3) = 0.;
306  outputValues(6, i0, 4) = x*(2. - 4.*z);
307  outputValues(6, i0, 5) = -4.*x*(-1. + x + y);
308 
309  outputValues(7, i0, 0) = 0.;
310  outputValues(7, i0, 1) = 2.*(-1. + z)*z;
311  outputValues(7, i0, 2) = 2.*y*(-1. + 2.*z);
312  outputValues(7, i0, 3) = 0.;
313  outputValues(7, i0, 4) = 2.*x*(-1. + 2.*z);
314  outputValues(7, i0, 5) = 4.*x*y;
315 
316  outputValues(8, i0, 0) = 0.;
317  outputValues(8, i0, 1) = -2.*(-1. + z)*z;
318  outputValues(8, i0, 2) = y*(2. - 4.*z);
319  outputValues(8, i0, 3) = -4.*(-1. + z)*z;
320  outputValues(8, i0, 4) = -2.*(-1. + x + 2.*y)*(-1. + 2.*z);
321  outputValues(8, i0, 5) = -4.*y*(-1. + x + y);
322 
323  outputValues(9, i0, 0) = 4. - 4.*z*z;
324  outputValues(9, i0, 1) = 4. - 4.*z*z;
325  outputValues(9, i0, 2) = -2.*(-3. + 4.*x + 4.*y)*z;
326  outputValues(9, i0, 3) = 4. - 4.*z*z;
327  outputValues(9, i0, 4) = -2.*(-3. + 4.*x + 4.*y)*z;
328  outputValues(9, i0, 5) = -2.*(-1. + x + y)*(-1. + 2.*x + 2.*y);
329 
330  outputValues(10,i0, 0) = 4. - 4.*z*z;
331  outputValues(10,i0, 1) = 0.;
332  outputValues(10,i0, 2) = (2. - 8.*x)*z;
333  outputValues(10,i0, 3) = 0.;
334  outputValues(10,i0, 4) = 0.;
335  outputValues(10,i0, 5) = -2.*x*(-1. + 2.*x);
336 
337  outputValues(11,i0, 0) = 0.;
338  outputValues(11,i0, 1) = 0.;
339  outputValues(11,i0, 2) = 0.;
340  outputValues(11,i0, 3) = 4. - 4.*z*z;
341  outputValues(11,i0, 4) = (2. - 8.*y)*z;
342  outputValues(11,i0, 5) = -2.*y*(-1. + 2.*y);
343 
344  outputValues(12,i0, 0) = -4.*z*(1. + z);
345  outputValues(12,i0, 1) = -2.*z*(1. + z);
346  outputValues(12,i0, 2) = -2.*(-1. + 2.*x + y)*(1. + 2.*z);
347  outputValues(12,i0, 3) = 0.;
348  outputValues(12,i0, 4) = -2.*(x + 2.*x*z);
349  outputValues(12,i0, 5) = -4.*x*(-1. + x + y);
350 
351  outputValues(13,i0, 0) = 0.;
352  outputValues(13,i0, 1) = 2.*z*(1. + z);
353  outputValues(13,i0, 2) = 2.*(y + 2.*y*z);
354  outputValues(13,i0, 3) = 0.;
355  outputValues(13,i0, 4) = 2.*(x + 2.*x*z);
356  outputValues(13,i0, 5) = 4.*x*y;
357 
358  outputValues(14,i0, 0) = 0.;
359  outputValues(14,i0, 1) = -2.*z*(1. + z);
360  outputValues(14,i0, 2) = -2.*(y + 2.*y*z);
361  outputValues(14,i0, 3) = -4.*z*(1. + z);
362  outputValues(14,i0, 4) = -2.*(-1. + x + 2.*y)*(1. + 2.*z);
363  outputValues(14,i0, 5) = -4.*y*(-1. + x + y);
364 
365  outputValues(15,i0, 0) = 8.*(-1. + z*z);
366  outputValues(15,i0, 1) = 4.*(-1. + z*z);
367  outputValues(15,i0, 2) = 8.*(-1. + 2.*x + y)*z;
368  outputValues(15,i0, 3) = 0.;
369  outputValues(15,i0, 4) = 8.*x*z;
370  outputValues(15,i0, 5) = 8.*x*(-1. + x + y);
371 
372  outputValues(16,i0, 0) = 0.;
373  outputValues(16,i0, 1) = 4. - 4.*z*z;
374  outputValues(16,i0, 2) = -8.*y*z;
375  outputValues(16,i0, 3) = 0.;
376  outputValues(16,i0, 4) = -8.*x*z;
377  outputValues(16,i0, 5) = -8.*x*y;
378 
379 
380  outputValues(17,i0, 0) = 0.;
381  outputValues(17,i0, 1) = 4.*(-1. + z*z);
382  outputValues(17,i0, 2) = 8.*y*z;
383  outputValues(17,i0, 3) = 8.*(-1. + z*z);
384  outputValues(17,i0, 4) = 8.*(-1. + x + 2.*y)*z;
385  outputValues(17,i0, 5) = 8.*y*(-1. + x + y);
386  }
387  break;
388 
389  case OPERATOR_D3:
390  for (int i0 = 0; i0 < dim0; i0++) {
391  x = inputPoints(i0,0);
392  y = inputPoints(i0,1);
393  z = inputPoints(i0,2);
394 
395  outputValues(0, i0, 0) = 0.;
396  outputValues(0, i0, 1) = 0.;
397  outputValues(0, i0, 2) = -2. + 4.*z;
398  outputValues(0, i0, 3) = 0.;
399  outputValues(0, i0, 4) = -2. + 4.*z;
400  outputValues(0, i0, 5) = -3. + 4.*x + 4.*y;
401  outputValues(0, i0, 6) = 0.;
402  outputValues(0, i0, 7) = -2. + 4.*z;
403  outputValues(0, i0, 8) = -3. + 4.*x + 4.*y;
404  outputValues(0, i0, 9) = 0.;
405 
406  outputValues(1, i0, 0) = 0.;
407  outputValues(1, i0, 1) = 0.;
408  outputValues(1, i0, 2) = -2. + 4.*z;
409  outputValues(1, i0, 3) = 0.;
410  outputValues(1, i0, 4) = 0.;
411  outputValues(1, i0, 5) = -1 + 4.*x;
412  outputValues(1, i0, 6) = 0.;
413  outputValues(1, i0, 7) = 0.;
414  outputValues(1, i0, 8) = 0.;
415  outputValues(1, i0, 9) = 0.;
416 
417  outputValues(2, i0, 0) = 0.;
418  outputValues(2, i0, 1) = 0.;
419  outputValues(2, i0, 2) = 0.;
420  outputValues(2, i0, 3) = 0.;
421  outputValues(2, i0, 4) = 0.;
422  outputValues(2, i0, 5) = 0.;
423  outputValues(2, i0, 6) = 0.;
424  outputValues(2, i0, 7) = -2. + 4.*z;
425  outputValues(2, i0, 8) = -1 + 4.*y;
426  outputValues(2, i0, 9) = 0.;
427 
428  outputValues(3, i0, 0) = 0.;
429  outputValues(3, i0, 1) = 0.;
430  outputValues(3, i0, 2) = 2. + 4.*z;
431  outputValues(3, i0, 3) = 0.;
432  outputValues(3, i0, 4) = 2. + 4.*z;
433  outputValues(3, i0, 5) = -3. + 4.*x + 4.*y;
434  outputValues(3, i0, 6) = 0.;
435  outputValues(3, i0, 7) = 2. + 4.*z;
436  outputValues(3, i0, 8) = -3. + 4.*x + 4.*y;
437  outputValues(3, i0, 9) = 0.;
438 
439  outputValues(4, i0, 0) = 0.;
440  outputValues(4, i0, 1) = 0.;
441  outputValues(4, i0, 2) = 2. + 4.*z;
442  outputValues(4, i0, 3) = 0.;
443  outputValues(4, i0, 4) = 0.;
444  outputValues(4, i0, 5) = -1 + 4.*x;
445  outputValues(4, i0, 6) = 0.;
446  outputValues(4, i0, 7) = 0.;
447  outputValues(4, i0, 8) = 0.;
448  outputValues(4, i0, 9) = 0.;
449 
450  outputValues(5, i0, 0) = 0.;
451  outputValues(5, i0, 1) = 0.;
452  outputValues(5, i0, 2) = 0.;
453  outputValues(5, i0, 3) = 0.;
454  outputValues(5, i0, 4) = 0.;
455  outputValues(5, i0, 5) = 0.;
456  outputValues(5, i0, 6) = 0.;
457  outputValues(5, i0, 7) = 2. + 4.*z;
458  outputValues(5, i0, 8) = -1 + 4.*y;
459  outputValues(5, i0, 9) = 0.;
460 
461  outputValues(6, i0, 0) = 0.;
462  outputValues(6, i0, 1) = 0.;
463  outputValues(6, i0, 2) = 4. - 8.*z;
464  outputValues(6, i0, 3) = 0.;
465  outputValues(6, i0, 4) = 2. - 4.*z;
466  outputValues(6, i0, 5) = -4.*(-1 + 2*x + y);
467  outputValues(6, i0, 6) = 0.;
468  outputValues(6, i0, 7) = 0.;
469  outputValues(6, i0, 8) = -4.*x;
470  outputValues(6, i0, 9) = 0.;
471 
472  outputValues(7, i0, 0) = 0.;
473  outputValues(7, i0, 1) = 0.;
474  outputValues(7, i0, 2) = 0.;
475  outputValues(7, i0, 3) = 0.;
476  outputValues(7, i0, 4) = -2. + 4.*z;
477  outputValues(7, i0, 5) = 4.*y;
478  outputValues(7, i0, 6) = 0.;
479  outputValues(7, i0, 7) = 0.;
480  outputValues(7, i0, 8) = 4.*x;
481  outputValues(7, i0, 9) = 0.;
482 
483  outputValues(8, i0, 0) = 0.;
484  outputValues(8, i0, 1) = 0.;
485  outputValues(8, i0, 2) = 0.;
486  outputValues(8, i0, 3) = 0.;
487  outputValues(8, i0, 4) = 2. - 4.*z;
488  outputValues(8, i0, 5) = -4.*y;
489  outputValues(8, i0, 6) = 0.;
490  outputValues(8, i0, 7) = 4. - 8.*z;
491  outputValues(8, i0, 8) = -4.*(-1 + x + 2*y);
492  outputValues(8, i0, 9) = 0.;
493 
494  outputValues(9, i0, 0) = 0.;
495  outputValues(9, i0, 1) = 0.;
496  outputValues(9, i0, 2) = -8.*z;
497  outputValues(9, i0, 3) = 0.;
498  outputValues(9, i0, 4) = -8.*z;
499  outputValues(9, i0, 5) = 6. - 8.*x - 8.*y;
500  outputValues(9, i0, 6) = 0.;
501  outputValues(9, i0, 7) = -8.*z;
502  outputValues(9, i0, 8) = 6. - 8.*x - 8.*y;
503  outputValues(9, i0, 9) = 0.;
504 
505  outputValues(10,i0, 0) = 0.;
506  outputValues(10,i0, 1) = 0.;
507  outputValues(10,i0, 2) = -8.*z;
508  outputValues(10,i0, 3) = 0.;
509  outputValues(10,i0, 4) = 0.;
510  outputValues(10,i0, 5) = 2. - 8.*x;
511  outputValues(10,i0, 6) = 0.;
512  outputValues(10,i0, 7) = 0.;
513  outputValues(10,i0, 8) = 0.;
514  outputValues(10,i0, 9) = 0.;
515 
516  outputValues(11,i0, 0) = 0.;
517  outputValues(11,i0, 1) = 0.;
518  outputValues(11,i0, 2) = 0.;
519  outputValues(11,i0, 3) = 0.;
520  outputValues(11,i0, 4) = 0.;
521  outputValues(11,i0, 5) = 0.;
522  outputValues(11,i0, 6) = 0.;
523  outputValues(11,i0, 7) = -8.*z;
524  outputValues(11,i0, 8) = 2. - 8.*y;
525  outputValues(11,i0, 9) = 0.;
526 
527  outputValues(12,i0, 0) = 0.;
528  outputValues(12,i0, 1) = 0.;
529  outputValues(12,i0, 2) = -4. - 8.*z;
530  outputValues(12,i0, 3) = 0.;
531  outputValues(12,i0, 4) = -2. - 4.*z;
532  outputValues(12,i0, 5) = -4.*(-1 + 2*x + y);
533  outputValues(12,i0, 6) = 0.;
534  outputValues(12,i0, 7) = 0.;
535  outputValues(12,i0, 8) = -4.*x;
536  outputValues(12,i0, 9) = 0.;
537 
538  outputValues(13,i0, 0) = 0.;
539  outputValues(13,i0, 1) = 0.;
540  outputValues(13,i0, 2) = 0.;
541  outputValues(13,i0, 3) = 0.;
542  outputValues(13,i0, 4) = 2. + 4.*z;
543  outputValues(13,i0, 5) = 4.*y;
544  outputValues(13,i0, 6) = 0.;
545  outputValues(13,i0, 7) = 0.;
546  outputValues(13,i0, 8) = 4.*x;
547  outputValues(13,i0, 9) = 0.;
548 
549  outputValues(14,i0, 0) = 0.;
550  outputValues(14,i0, 1) = 0.;
551  outputValues(14,i0, 2) = 0.;
552  outputValues(14,i0, 3) = 0.;
553  outputValues(14,i0, 4) = -2. - 4.*z;
554  outputValues(14,i0, 5) = -4.*y;
555  outputValues(14,i0, 6) = 0.;
556  outputValues(14,i0, 7) = -4. - 8.*z;
557  outputValues(14,i0, 8) = -4.*(-1 + x + 2*y);
558  outputValues(14,i0, 9) = 0.;
559 
560  outputValues(15,i0, 0) = 0.;
561  outputValues(15,i0, 1) = 0.;
562  outputValues(15,i0, 2) = 16.*z;
563  outputValues(15,i0, 3) = 0.;
564  outputValues(15,i0, 4) = 8.*z;
565  outputValues(15,i0, 5) = 8.*(-1 + 2*x + y);
566  outputValues(15,i0, 6) = 0.;
567  outputValues(15,i0, 7) = 0.;
568  outputValues(15,i0, 8) = 8.*x;
569  outputValues(15,i0, 9) = 0.;
570 
571  outputValues(16,i0, 0) = 0.;
572  outputValues(16,i0, 1) = 0.;
573  outputValues(16,i0, 2) = 0.;
574  outputValues(16,i0, 3) = 0.;
575  outputValues(16,i0, 4) = -8.*z;
576  outputValues(16,i0, 5) = -8.*y;
577  outputValues(16,i0, 6) = 0.;
578  outputValues(16,i0, 7) = 0.;
579  outputValues(16,i0, 8) = -8.*x;
580  outputValues(16,i0, 9) = 0.;
581 
582  outputValues(17,i0, 0) = 0.;
583  outputValues(17,i0, 1) = 0.;
584  outputValues(17,i0, 2) = 0.;
585  outputValues(17,i0, 3) = 0.;
586  outputValues(17,i0, 4) = 8.*z;
587  outputValues(17,i0, 5) = 8.*y;
588  outputValues(17,i0, 6) = 0.;
589  outputValues(17,i0, 7) = 16.*z;
590  outputValues(17,i0, 8) = 8.*(-1 + x + 2*y);
591  outputValues(17,i0, 9) = 0.;
592 
593  }
594  break;
595 
596  case OPERATOR_D4:
597  {
598  // There are only few constant non-zero entries. Initialize by zero and then assign non-zero entries.
599  int DkCardinality = Intrepid::getDkCardinality(operatorType, this -> basisCellTopology_.getDimension() );
600  for(int dofOrd = 0; dofOrd < this -> basisCardinality_; dofOrd++) {
601  for (int i0 = 0; i0 < dim0; i0++) {
602  for(int dkOrd = 0; dkOrd < DkCardinality; dkOrd++){
603  outputValues(dofOrd, i0, dkOrd) = 0.0;
604  }
605  }
606  }
607 
608  for (int i0 = 0; i0 < dim0; i0++) {
609 
610  outputValues(0, i0, 5) = 4.;
611  outputValues(0, i0, 8) = 4.;
612  outputValues(0, i0,12) = 4.;
613 
614  outputValues(1, i0, 5) = 4.;
615 
616  outputValues(2, i0,12) = 4.;
617 
618  outputValues(3, i0, 5) = 4.;
619  outputValues(3, i0, 8) = 4.;
620  outputValues(3, i0,12) = 4.;
621 
622  outputValues(4, i0, 5) = 4.0;
623 
624  outputValues(5, i0,12) = 4.0;
625 
626  outputValues(6, i0, 5) =-8.;
627  outputValues(6, i0, 8) =-4.;
628 
629  outputValues(7, i0, 8) = 4.;
630 
631  outputValues(8, i0, 8) =-4.;
632  outputValues(8, i0,12) =-8.;
633 
634  outputValues(9, i0, 5) =-8.;
635  outputValues(9, i0, 8) =-8.;
636  outputValues(9, i0,12) =-8.;
637 
638  outputValues(10,i0, 5) =-8.;
639 
640  outputValues(11,i0,12) =-8.;
641 
642  outputValues(12,i0, 5) =-8.;
643  outputValues(12,i0, 8) =-4.;
644 
645  outputValues(13,i0, 8) = 4.;
646 
647  outputValues(14,i0, 8) =-4;
648  outputValues(14,i0,12) =-8.;
649 
650  outputValues(15,i0, 5) =16.;
651  outputValues(15,i0, 8) = 8.;
652 
653  outputValues(16,i0, 8) =-8.;
654 
655 
656  outputValues(17,i0, 8) = 8.;
657  outputValues(17,i0,12) =16.;
658  }
659  }
660  break;
661 
662  case OPERATOR_D5:
663  case OPERATOR_D6:
664  case OPERATOR_D7:
665  case OPERATOR_D8:
666  case OPERATOR_D9:
667  case OPERATOR_D10:
668  {
669  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, DkCardinality)
670  int DkCardinality = Intrepid::getDkCardinality(operatorType,
671  this -> basisCellTopology_.getDimension() );
672  for(int dofOrd = 0; dofOrd < this -> basisCardinality_; dofOrd++) {
673  for (int i0 = 0; i0 < dim0; i0++) {
674  for(int dkOrd = 0; dkOrd < DkCardinality; dkOrd++){
675  outputValues(dofOrd, i0, dkOrd) = 0.0;
676  }
677  }
678  }
679  }
680  break;
681 
682  default:
683  TEUCHOS_TEST_FOR_EXCEPTION( !( Intrepid::isValidOperator(operatorType) ), std::invalid_argument,
684  ">>> ERROR (Basis_HGRAD_WEDGE_C2_FEM): Invalid operator type");
685  }
686 }
687 
688 
689 
690 template<class Scalar, class ArrayScalar>
692  const ArrayScalar & inputPoints,
693  const ArrayScalar & cellVertices,
694  const EOperator operatorType) const {
695  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
696  ">>> ERROR (Basis_HGRAD_WEDGE_C2_FEM): FEM Basis calling an FVD member function");
697 }
698 }// namespace Intrepid
699 #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.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
FEM basis evaluation on a reference Wedge cell.
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.