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 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.