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
51namespace 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
65template<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
108template<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
690template<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
int isValidOperator(const EOperator operatorType)
Verifies validity of an operator enum.
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.
int getDkCardinality(const EOperator operatorType, const int spaceDim)
Returns cardinality of Dk, i.e., the number of all derivatives of order k.
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.