Intrepid2
Intrepid2_BasisDef.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Intrepid2 Package
4//
5// Copyright 2007 NTESS and the Intrepid2 contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
16#ifndef __INTREPID2_BASIS_DEF_HPP__
17#define __INTREPID2_BASIS_DEF_HPP__
18
19#include <stdexcept>
20
21namespace Intrepid2 {
22
23 //--------------------------------------------------------------------------------------------//
24 // //
25 // Helper functions of the Basis class //
26 // //
27 //--------------------------------------------------------------------------------------------//
28
29 //
30 // functions for orders, cardinality and etc. of differential operators
31 //
32
33 KOKKOS_INLINE_FUNCTION
34 ordinal_type getFieldRank(const EFunctionSpace spaceType) {
35 ordinal_type fieldRank = -1;
36
37 switch (spaceType) {
38
39 case FUNCTION_SPACE_HGRAD:
40 case FUNCTION_SPACE_HVOL:
41 fieldRank = 0;
42 break;
43
44 case FUNCTION_SPACE_HCURL:
45 case FUNCTION_SPACE_HDIV:
46 case FUNCTION_SPACE_VECTOR_HGRAD:
47 fieldRank = 1;
48 break;
49
50 case FUNCTION_SPACE_TENSOR_HGRAD:
51 fieldRank = 2;
52 break;
53
54 default:
55 INTREPID2_TEST_FOR_ABORT( !isValidFunctionSpace(spaceType),
56 ">>> ERROR (Intrepid2::getFieldRank): Invalid function space type");
57 }
58 return fieldRank;
59 }
60
61
62 KOKKOS_INLINE_FUNCTION
63 ordinal_type getOperatorRank(const EFunctionSpace spaceType,
64 const EOperator operatorType,
65 const ordinal_type spaceDim) {
66
67 const auto fieldRank = Intrepid2::getFieldRank(spaceType);
68#ifdef HAVE_INTREPID2_DEBUG
69 // Verify arguments: field rank can be 0,1, or 2, spaceDim can be 1,2, or 3.
70 INTREPID2_TEST_FOR_ABORT( !(0 <= fieldRank && fieldRank <= 2),
71 ">>> ERROR (Intrepid2::getOperatorRank): Invalid field rank");
72 INTREPID2_TEST_FOR_ABORT( !(1 <= spaceDim && spaceDim <= 3),
73 ">>> ERROR (Intrepid2::getOperatorRank): Invalid space dimension");
74#endif
75 ordinal_type operatorRank = -999;
76
77 // In 1D GRAD, CURL, and DIV default to d/dx; Dk defaults to d^k/dx^k, no casing needed.
78 if (spaceDim == 1) {
79 if (fieldRank == 0) {
80 // By default, in 1D any operator other than VALUE has rank 1
81 if (operatorType == OPERATOR_VALUE) {
82 operatorRank = 0;
83 } else {
84 operatorRank = 1;
85 }
86 }
87
88 // Only scalar fields are allowed in 1D
89 else {
90 INTREPID2_TEST_FOR_ABORT( fieldRank > 0,
91 ">>> ERROR (getOperatorRank): Only scalar fields are allowed in 1D");
92 } // fieldRank == 0
93 } // spaceDim == 1
94
95 // We are either in 2D or 3D
96 else {
97 switch (operatorType) {
98 case OPERATOR_VALUE:
99 operatorRank = 0;
100 break;
101
102 case OPERATOR_GRAD:
103 case OPERATOR_D1:
104 case OPERATOR_D2:
105 case OPERATOR_D3:
106 case OPERATOR_D4:
107 case OPERATOR_D5:
108 case OPERATOR_D6:
109 case OPERATOR_D7:
110 case OPERATOR_D8:
111 case OPERATOR_D9:
112 case OPERATOR_D10:
113 operatorRank = 1;
114 break;
115
116 case OPERATOR_CURL:
117
118 // operator rank for vector and tensor fields equals spaceDim - 3 (-1 in 2D and 0 in 3D)
119 if (fieldRank > 0) {
120 operatorRank = spaceDim - 3;
121 } else {
122
123 // CURL can be applied to scalar functions (rank = 0) in 2D and gives a vector (rank = 1)
124 if (spaceDim == 2) {
125 operatorRank = 3 - spaceDim;
126 }
127
128 // If we are here, fieldRank=0, spaceDim=3: CURL is undefined for 3D scalar functions
129 else {
130 INTREPID2_TEST_FOR_ABORT( ( (spaceDim == 3) && (fieldRank == 0) ),
131 ">>> ERROR (Intrepid2::getOperatorRank): CURL cannot be applied to scalar fields in 3D");
132 }
133 }
134 break;
135
136 case OPERATOR_DIV:
137
138 // DIV can be applied to vectors and tensors and has rank -1 in 2D and 3D
139 if (fieldRank > 0) {
140 operatorRank = -1;
141 }
142
143 // DIV cannot be applied to scalar fields except in 1D where it defaults to d/dx
144 else {
145 INTREPID2_TEST_FOR_ABORT( ( (spaceDim > 1) && (fieldRank == 0) ),
146 ">>> ERROR (Intrepid2::getOperatorRank): DIV cannot be applied to scalar fields in 2D and 3D");
147 }
148 break;
149
150 default:
151 INTREPID2_TEST_FOR_ABORT( !( isValidOperator(operatorType) ),
152 ">>> ERROR (Intrepid2::getOperatorRank): Invalid operator type");
153 } // switch
154 }// 2D and 3D
155
156 return operatorRank;
157 }
158
159
160 KOKKOS_INLINE_FUNCTION
161 ordinal_type getOperatorOrder(const EOperator operatorType) {
162 ordinal_type opOrder = -1;
163
164 switch (operatorType) {
165
166 case OPERATOR_VALUE:
167 opOrder = 0;
168 break;
169
170 case OPERATOR_GRAD:
171 case OPERATOR_CURL:
172 case OPERATOR_DIV:
173 case OPERATOR_D1:
174 opOrder = 1;
175 break;
176
177 case OPERATOR_D2:
178 case OPERATOR_D3:
179 case OPERATOR_D4:
180 case OPERATOR_D5:
181 case OPERATOR_D6:
182 case OPERATOR_D7:
183 case OPERATOR_D8:
184 case OPERATOR_D9:
185 case OPERATOR_D10:
186 opOrder = (ordinal_type)operatorType - (ordinal_type)OPERATOR_D1 + 1;
187 break;
188
189 default:
190 INTREPID2_TEST_FOR_ABORT( !( Intrepid2::isValidOperator(operatorType) ),
191 ">>> ERROR (Intrepid2::getOperatorOrder): Invalid operator type");
192 }
193 return opOrder;
194 }
195
196 template<EOperator operatorType>
197 KOKKOS_INLINE_FUNCTION
198 constexpr ordinal_type getOperatorOrder() {
199 return (operatorType == OPERATOR_VALUE) ? 0 :
200 ((operatorType == OPERATOR_GRAD) || (operatorType == OPERATOR_CURL) || (operatorType == OPERATOR_DIV) || (operatorType == OPERATOR_D1)) ? 1 :
201 (ordinal_type)operatorType - (ordinal_type)OPERATOR_D1 + 1;
202 }
203
204
205 template<ordinal_type spaceDim>
206 KOKKOS_INLINE_FUNCTION
207 ordinal_type getDkEnumeration(const ordinal_type /*xMult*/,
208 const ordinal_type yMult,
209 const ordinal_type zMult) {
210 switch(spaceDim) {
211
212 case 1: return 0;
213 case 2: return yMult;
214 case 3: return zMult + (yMult+zMult)*(yMult+zMult+1)/2;
215
216 default: {
217 INTREPID2_TEST_FOR_ABORT( !( (0 < spaceDim ) && (spaceDim < 4) ),
218 ">>> ERROR (Intrepid2::getDkEnumeration): Invalid space dimension");
219 return 0;
220 }
221 }
222 }
223
224 template<ordinal_type spaceDim>
225 KOKKOS_INLINE_FUNCTION
226 ordinal_type getPnEnumeration(const ordinal_type p,
227 const ordinal_type q /* = 0*/,
228 const ordinal_type r /* = 0*/) {
229 return (spaceDim==1) ? p :
230 (spaceDim==2) ? (p+q)*(p+q+1)/2+q :
231 (p+q+r)*(p+q+r+1)*(p+q+r+2)/6+(q+r)*(q+r+1)/2+r;
232
233 }
234
235 template<typename value_type>
236 KOKKOS_INLINE_FUNCTION
238 value_type &an,
239 value_type &bn,
240 value_type &cn,
241 const ordinal_type alpha,
242 const ordinal_type beta ,
243 const ordinal_type n) {
244 an = ( (2.0 * n + 1.0 + alpha + beta) * ( 2.0 * n + 2.0 + alpha + beta ) /
245 value_type(2.0 * ( n + 1 ) * ( n + 1 + alpha + beta ) ) );
246 bn = ( (alpha*alpha-beta*beta)*(2.0*n+1.0+alpha+beta) /
247 value_type(2.0*(n+1.0)*(2.0*n+alpha+beta)*(n+1.0+alpha+beta) ) );
248 cn = ( (n+alpha)*(n+beta)*(2.0*n+2.0+alpha+beta) /
249 value_type( (n+1.0)*(n+1.0+alpha+beta)*(2.0*n+alpha+beta) ) );
250 }
251
252
253// template<typename OrdinalArrayType>
254// KOKKOS_INLINE_FUNCTION
255// void getDkMultiplicities(OrdinalArrayType partialMult,
256// const ordinal_type derivativeEnum,
257// const EOperator operatorType,
258// const ordinal_type spaceDim) {
259
260// /* Hash table to convert enumeration of partial derivative to multiplicities of dx,dy,dz in 3D.
261// Multiplicities {mx,my,mz} are arranged lexicographically in bins numbered from 0 to 10.
262// The size of bins is an arithmetic progression, i.e., 1,2,3,4,5,...,11. Conversion formula is:
263// \verbatim
264// mx = derivativeOrder - binNumber
265// mz = derivativeEnum - binBegin
266// my = derivativeOrder - mx - mz = binNumber + binBegin - derivativeEnum
267// \endverbatim
268// where binBegin is the enumeration of the first element in the bin. Bin numbers and binBegin
269// values are stored in hash tables for quick access by derivative enumeration value.
270// */
271
272// // Returns the bin number for the specified derivative enumeration
273// const ordinal_type binNumber[66] = {
274// 0,
275// 1, 1,
276// 2, 2, 2,
277// 3, 3, 3, 3,
278// 4, 4, 4, 4, 4,
279// 5, 5, 5, 5, 5, 5,
280// 6, 6, 6, 6, 6, 6, 6,
281// 7, 7, 7, 7, 7, 7, 7, 7,
282// 8, 8, 8, 8, 8, 8, 8, 8, 8,
283// 9, 9, 9, 9, 9, 9, 9, 9, 9, 9,
284// 10,10,10,10,10,10,10,10,10,10,10
285// };
286
287// // Returns the binBegin value for the specified derivative enumeration
288// const ordinal_type binBegin[66] = {
289// 0,
290// 1, 1,
291// 3, 3 ,3,
292// 6, 6, 6, 6,
293// 10,10,10,10,10,
294// 15,15,15,15,15,15,
295// 21,21,21,21,21,21,21,
296// 28,28,28,28,28,28,28,28,
297// 36,36,36,36,36,36,36,36,36,
298// 45,45,45,45,45,45,45,45,45,45,
299// 55,55,55,55,55,55,55,55,55,55,55
300// };
301
302// #ifdef HAVE_INTREPID2_DEBUG
303// // Enumeration value must be between 0 and the cardinality of the derivative set
304// INTREPID2_TEST_FOR_ABORT( !( (0 <= derivativeEnum) && (derivativeEnum < getDkCardinality(operatorType,spaceDim) ) ),
305// ">>> ERROR (Intrepid2::getDkMultiplicities): Invalid derivative enumeration value for this order and space dimension");
306// #endif
307
308// // This method should only be called for Dk operators
309// ordinal_type derivativeOrder;
310// switch(operatorType) {
311
312// case OPERATOR_D1:
313// case OPERATOR_D2:
314// case OPERATOR_D3:
315// case OPERATOR_D4:
316// case OPERATOR_D5:
317// case OPERATOR_D6:
318// case OPERATOR_D7:
319// case OPERATOR_D8:
320// case OPERATOR_D9:
321// case OPERATOR_D10:
322// derivativeOrder = Intrepid2::getOperatorOrder(operatorType);
323// break;
324
325// default:
326// INTREPID2_TEST_FOR_ABORT(true,
327// ">>> ERROR (Intrepid2::getDkMultiplicities): operator type Dk required for this method");
328// }// switch
329
330// #ifdef HAVE_INTREPID2_DEBUG
331// INTREPID2_TEST_FOR_ABORT( partialMult.extent(0) != spaceDim,
332// ">>> ERROR (Intrepid2::getDkMultiplicities): partialMult must have the same dimension as spaceDim" );
333// #endif
334
335// switch(spaceDim) {
336
337// case 1:
338// // Multiplicity of dx equals derivativeOrder
339// partialMult(0) = derivativeOrder;
340// break;
341
342// case 2:
343// // Multiplicity of dy equals the enumeration of the derivative; of dx - the complement
344// partialMult(1) = derivativeEnum;
345// partialMult(0) = derivativeOrder - derivativeEnum;
346// break;
347
348// case 3:
349// // Recover multiplicities
350// partialMult(0) = derivativeOrder - binNumber[derivativeEnum];
351// partialMult(1) = binNumber[derivativeEnum] + binBegin[derivativeEnum] - derivativeEnum;
352// partialMult(2) = derivativeEnum - binBegin[derivativeEnum];
353// break;
354
355// default:
356// INTREPID2_TEST_FOR_ABORT( !( (0 < spaceDim ) && (spaceDim < 4) ),
357// ">>> ERROR (Intrepid2::getDkMultiplicities): Invalid space dimension");
358// }
359// }
360
361
362 KOKKOS_INLINE_FUNCTION
363 ordinal_type getDkCardinality(const EOperator operatorType,
364 const ordinal_type spaceDim) {
365
366#ifdef HAVE_INTREPID2_DEBUG
367 INTREPID2_TEST_FOR_ABORT( !( (0 < spaceDim ) && (spaceDim < 8) ),
368 ">>> ERROR (Intrepid2::getDkcardinality): Invalid space dimension");
369 switch (operatorType) {
370 case OPERATOR_VALUE:
371 case OPERATOR_GRAD:
372 case OPERATOR_D1:
373 case OPERATOR_D2:
374 case OPERATOR_D3:
375 case OPERATOR_D4:
376 case OPERATOR_D5:
377 case OPERATOR_D6:
378 case OPERATOR_D7:
379 case OPERATOR_D8:
380 case OPERATOR_D9:
381 case OPERATOR_D10:
382 break;
383 default:
384 INTREPID2_TEST_FOR_ABORT( true, ">>> ERROR (Intrepid2::getDkCardinality): Cannot be used for this operator ");
385 break;
386 }
387#endif
388
389 ordinal_type n = Intrepid2::getOperatorOrder(operatorType);
390 return (spaceDim==1) ? 1 :
391 (spaceDim==2) ? n + 1 :
392 (spaceDim==3) ? (n + 1) * (n + 2) / 2 :
393 (spaceDim==4) ? (n + 1) * (n + 2) * (n + 3) / 6 :
394 (spaceDim==5) ? (n + 1) * (n + 2) * (n + 3) * (n + 4) / 24 :
395 (spaceDim==6) ? (n + 1) * (n + 2) * (n + 3) * (n + 4) * (n + 5) / 120 :
396 (n + 1) * (n + 2) * (n + 3) * (n + 4) * (n + 5) * (n + 6) / 720;
397 }
398
399 template<EOperator operatorType, ordinal_type spaceDim>
400 KOKKOS_INLINE_FUNCTION
401 constexpr ordinal_type getDkCardinality() {
402 return getPnCardinality<spaceDim-1,Intrepid2::getOperatorOrder<operatorType>()>();
403 }
404
405 template<ordinal_type spaceDim>
406 KOKKOS_INLINE_FUNCTION
407 ordinal_type getPnCardinality (ordinal_type n) {
408
409#ifdef HAVE_INTREPID2_DEBUG
410 INTREPID2_TEST_FOR_ABORT( !( (0 <= spaceDim ) && (spaceDim < 4) ),
411 ">>> ERROR (Intrepid2::getPnCardinality): Invalid space dimension");
412#endif
413
414 return (spaceDim==0) ? 1 :
415 (spaceDim==1) ? n+1 :
416 (spaceDim==2) ? (n + 1) * (n + 2) / 2 :
417 (n + 1) * (n + 2) * (n + 3) / 6;
418 }
419
420
421 template<ordinal_type spaceDim, ordinal_type n>
422 KOKKOS_INLINE_FUNCTION
423 constexpr ordinal_type getPnCardinality () {
424
425 return (spaceDim==0) ? 1 :
426 (spaceDim==1) ? n+1 :
427 (spaceDim==2) ? (n + 1) * (n + 2) / 2 :
428 (n + 1) * (n + 2) * (n + 3) / 6;
429
430 }
431
432
433
434
435 //
436 // Argument check
437 //
438
439
440 template<typename outputValueViewType,
441 typename inputPointViewType>
442 void getValues_HGRAD_Args( const outputValueViewType outputValues,
443 const inputPointViewType inputPoints,
444 const EOperator operatorType,
445 const shards::CellTopology cellTopo,
446 const ordinal_type basisCard ) {
447 const auto spaceDim = cellTopo.getDimension();
448
449 // Verify inputPoints array
450 INTREPID2_TEST_FOR_EXCEPTION( !(rank(inputPoints) == 2), std::invalid_argument,
451 ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) rank = 2 required for inputPoints array");
452
453 INTREPID2_TEST_FOR_EXCEPTION( (inputPoints.extent(0) <= 0), std::invalid_argument,
454 ">>> ERROR (Intrepid2::getValues_HGRAD_Args): dim 0 (number of points) > 0 required for inputPoints array");
455
456 INTREPID2_TEST_FOR_EXCEPTION( !(inputPoints.extent(1) == spaceDim), std::invalid_argument,
457 ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) dim 1 (spatial dimension) of inputPoints array does not match cell dimension");
458
459
460 // Verify that all inputPoints are in the reference cell
461 /*
462 INTREPID2_TEST_FOR_EXCEPTION( !CellTools<Scalar>::checkPointSetInclusion(inputPoints, cellTopo), std::invalid_argument,
463 ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) One or more points are outside the "
464 << cellTopo <<" reference cell");
465 */
466
467
468 // Verify that operatorType is admissible for HGRAD fields
469 INTREPID2_TEST_FOR_EXCEPTION( ( (spaceDim == 2) && (operatorType == OPERATOR_DIV) ), std::invalid_argument,
470 ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) DIV is invalid operator for rank-0 (scalar) fields in 2D.");
471
472 INTREPID2_TEST_FOR_EXCEPTION( ( (spaceDim == 3) && ( (operatorType == OPERATOR_DIV) ||
473 (operatorType == OPERATOR_CURL) ) ), std::invalid_argument,
474 ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) DIV and CURL are invalid operators for rank-0 (scalar) fields in 3D.");
475
476
477 // Check rank of outputValues (all operators are admissible in 1D) and its dim 2 when operator is
478 // GRAD, CURL (only in 2D), or Dk.
479
480 if(spaceDim == 1) {
481 switch(operatorType){
482 case OPERATOR_VALUE:
483 INTREPID2_TEST_FOR_EXCEPTION( !(rank(outputValues) == 2), std::invalid_argument,
484 ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) rank = 2 required for outputValues when operator = VALUE.");
485 break;
486 case OPERATOR_GRAD:
487 case OPERATOR_CURL:
488 case OPERATOR_DIV:
489 case OPERATOR_D1:
490 case OPERATOR_D2:
491 case OPERATOR_D3:
492 case OPERATOR_D4:
493 case OPERATOR_D5:
494 case OPERATOR_D6:
495 case OPERATOR_D7:
496 case OPERATOR_D8:
497 case OPERATOR_D9:
498 case OPERATOR_D10:
499 INTREPID2_TEST_FOR_EXCEPTION( !(rank(outputValues) == 3), std::invalid_argument,
500 ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) rank = 3 required for outputValues in 1D when operator = GRAD, CURL, DIV, or Dk.");
501
502 INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.extent(2) == 1 ),
503 std::invalid_argument,
504 ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) dim 2 of outputValues must equal 1 when operator = GRAD, CURL, DIV, or Dk.");
505 break;
506 default:
507 INTREPID2_TEST_FOR_EXCEPTION( (true), std::invalid_argument, ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) Invalid operator");
508 }
509 }
510 else if(spaceDim > 1) {
511 switch(operatorType){
512 case OPERATOR_VALUE:
513 INTREPID2_TEST_FOR_EXCEPTION( !(rank(outputValues) == 2), std::invalid_argument,
514 ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) rank = 2 required for outputValues when operator = VALUE.");
515 break;
516 case OPERATOR_GRAD:
517 case OPERATOR_CURL:
518 case OPERATOR_D1:
519 INTREPID2_TEST_FOR_EXCEPTION( !(rank(outputValues) == 3), std::invalid_argument,
520 ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) rank = 3 required for outputValues in 2D and 3D when operator = GRAD, CURL (in 2D), or Dk.");
521
522 INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.extent(2) == spaceDim ),
523 std::invalid_argument,
524 ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) dim 2 of outputValues must equal cell dimension when operator = GRAD, CURL (in 2D), or D1.");
525 break;
526 case OPERATOR_D2:
527 case OPERATOR_D3:
528 case OPERATOR_D4:
529 case OPERATOR_D5:
530 case OPERATOR_D6:
531 case OPERATOR_D7:
532 case OPERATOR_D8:
533 case OPERATOR_D9:
534 case OPERATOR_D10:
535 INTREPID2_TEST_FOR_EXCEPTION( !(rank(outputValues) == 3), std::invalid_argument,
536 ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) rank = 3 required for outputValues in 2D and 3D when operator = GRAD, CURL (in 2D), or Dk.");
537
538 INTREPID2_TEST_FOR_EXCEPTION( !(static_cast<ordinal_type>(outputValues.extent(2)) == getDkCardinality(operatorType, spaceDim)),
539 std::invalid_argument,
540 ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) dim 2 of outputValues must equal cardinality of the Dk multiset.");
541 break;
542 default:
543 INTREPID2_TEST_FOR_EXCEPTION( (true), std::invalid_argument, ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) Invalid operator");
544 }
545 }
546
547
548 // Verify dim 0 and dim 1 of outputValues
549 INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.extent(1) == inputPoints.extent(0) ),
550 std::invalid_argument,
551 ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) dim 1 (number of points) of outputValues must equal dim 0 of inputPoints.");
552
553 INTREPID2_TEST_FOR_EXCEPTION( !(static_cast<ordinal_type>(outputValues.extent(0)) == basisCard ),
554 std::invalid_argument,
555 ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) dim 0 (number of basis functions) of outputValues must equal basis cardinality.");
556 }
557
558
559 template<typename outputValueViewType,
560 typename inputPointViewType>
561 void getValues_HCURL_Args( const outputValueViewType outputValues,
562 const inputPointViewType inputPoints,
563 const EOperator operatorType,
564 const shards::CellTopology cellTopo,
565 const ordinal_type basisCard ) {
566
567 const auto spaceDim = cellTopo.getDimension();
568
569 // Verify that cell is 2D or 3D (this is redundant for default bases where we use correct cells,
570 // but will force any user-defined basis for HCURL spaces to use 2D or 3D cells
571 INTREPID2_TEST_FOR_EXCEPTION( !( (spaceDim == 2) || (spaceDim == 3) ), std::invalid_argument,
572 ">>> ERROR: (Intrepid2::getValues_HCURL_Args) cell dimension = 2 or 3 required for HCURL basis");
573
574
575 // Verify inputPoints array
576 INTREPID2_TEST_FOR_EXCEPTION( !(rank(inputPoints) == 2), std::invalid_argument,
577 ">>> ERROR: (Intrepid2::getValues_HCURL_Args) rank = 2 required for inputPoints array");
578 INTREPID2_TEST_FOR_EXCEPTION( (inputPoints.extent(0) <= 0), std::invalid_argument,
579 ">>> ERROR (Intrepid2::getValues_HCURL_Args): dim 0 (number of points) > 0 required for inputPoints array");
580
581 INTREPID2_TEST_FOR_EXCEPTION( !(inputPoints.extent(1) == spaceDim), std::invalid_argument,
582 ">>> ERROR: (Intrepid2::getValues_HCURL_Args) dim 1 (spatial dimension) of inputPoints array does not match cell dimension");
583
584 // Verify that all inputPoints are in the reference cell
585 /*
586 INTREPID2_TEST_FOR_EXCEPTION( !CellTools<Scalar>::checkPointSetInclusion(inputPoints, cellTopo), std::invalid_argument,
587 ">>> ERROR: (Intrepid2::getValues_HCURL_Args) One or more points are outside the "
588 << cellTopo <<" reference cell");
589 */
590
591 // Verify that operatorType is admissible for HCURL fields
592 INTREPID2_TEST_FOR_EXCEPTION( !( (operatorType == OPERATOR_VALUE) || (operatorType == OPERATOR_CURL) ), std::invalid_argument,
593 ">>> ERROR: (Intrepid2::getValues_HCURL_Args) operator = VALUE or CURL required for HCURL fields.");
594
595
596 // Check rank of outputValues: for VALUE should always be rank-3 array with (F,P,D) layout
597 switch(operatorType) {
598
599 case OPERATOR_VALUE:
600 INTREPID2_TEST_FOR_EXCEPTION( !(rank(outputValues) == 3), std::invalid_argument,
601 ">>> ERROR: (Intrepid2::getValues_HCURL_Args) rank = 3 required for outputValues when operator is VALUE");
602 INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.extent(2) == spaceDim ),
603 std::invalid_argument,
604 ">>> ERROR: (Intrepid2::getValues_HCURL_Args) dim 2 of outputValues must equal cell dimension when operator is VALUE.");
605 break;
606
607 case OPERATOR_CURL:
608
609 // in 3D we need an (F,P,D) container because CURL gives a vector field:
610 if(spaceDim == 3) {
611 INTREPID2_TEST_FOR_EXCEPTION( !(rank(outputValues) == 3 ) ,
612 std::invalid_argument,
613 ">>> ERROR: (Intrepid2::getValues_HCURL_Args) rank = 3 required for outputValues in 3D when operator is CURL");
614 INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.extent(2) == spaceDim),
615 std::invalid_argument,
616 ">>> ERROR: (Intrepid2::getValues_HCURL_Args) dim 2 of outputValues must equal cell dimension in 3D when operator is CURL.");
617 }
618 // In 2D we need an (F,P) container because CURL gives a scalar field
619 else if(spaceDim == 2) {
620 INTREPID2_TEST_FOR_EXCEPTION( !(rank(outputValues) == 2 ) ,
621 std::invalid_argument,
622 ">>> ERROR: (Intrepid2::getValues_HCURL_Args) rank = 2 required for outputValues in 2D when operator is CURL");
623 }
624 break;
625
626 default:
627 INTREPID2_TEST_FOR_EXCEPTION( (true), std::invalid_argument, ">>> ERROR: (Intrepid2::getValues_HCURL_Args) Invalid operator");
628 }
629
630
631 // Verify dim 0 and dim 1 of outputValues
632 INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.extent(1) == inputPoints.extent(0) ),
633 std::invalid_argument,
634 ">>> ERROR: (Intrepid2::getValues_HCURL_Args) dim 1 (number of points) of outputValues must equal dim 0 of inputPoints.");
635
636 INTREPID2_TEST_FOR_EXCEPTION( !(static_cast<ordinal_type>(outputValues.extent(0)) == basisCard ),
637 std::invalid_argument,
638 ">>> ERROR: (Intrepid2::getValues_HCURL_Args) dim 0 (number of basis functions) of outputValues must equal basis cardinality.");
639
640 }
641
642
643
644 template<typename outputValueViewType,
645 typename inputPointViewType>
646 void getValues_HDIV_Args( const outputValueViewType outputValues,
647 const inputPointViewType inputPoints,
648 const EOperator operatorType,
649 const shards::CellTopology cellTopo,
650 const ordinal_type basisCard ) {
651
652 const auto spaceDim = cellTopo.getDimension();
653
654 // Verify inputPoints array
655 INTREPID2_TEST_FOR_EXCEPTION( !(rank(inputPoints) == 2), std::invalid_argument,
656 ">>> ERROR: (Intrepid2::getValues_HDIV_Args) rank = 2 required for inputPoints array");
657 INTREPID2_TEST_FOR_EXCEPTION( (inputPoints.extent(0) <= 0), std::invalid_argument,
658 ">>> ERROR (Intrepid2::getValues_HDIV_Args): dim 0 (number of points) > 0 required for inputPoints array");
659
660 INTREPID2_TEST_FOR_EXCEPTION( !(inputPoints.extent(1) == spaceDim), std::invalid_argument,
661 ">>> ERROR: (Intrepid2::getValues_HDIV_Args) dim 1 (spatial dimension) of inputPoints array does not match cell dimension");
662
663 // Verify that all inputPoints are in the reference cell
664 /*
665 INTREPID2_TEST_FOR_EXCEPTION( !CellTools<Scalar>::checkPointSetInclusion(inputPoints, cellTopo), std::invalid_argument,
666 ">>> ERROR: (Intrepid2::getValues_HDIV_Args) One or more points are outside the "
667 << cellTopo <<" reference cell");
668 */
669
670 // Verify that operatorType is admissible for HDIV fields
671 INTREPID2_TEST_FOR_EXCEPTION( !( (operatorType == OPERATOR_VALUE) || (operatorType == OPERATOR_DIV) ), std::invalid_argument,
672 ">>> ERROR: (Intrepid2::getValues_HDIV_Args) operator = VALUE or DIV required for HDIV fields.");
673
674
675 // Check rank of outputValues
676 switch(operatorType) {
677 case OPERATOR_VALUE:
678 INTREPID2_TEST_FOR_EXCEPTION( !(rank(outputValues) == 3), std::invalid_argument,
679 ">>> ERROR: (Intrepid2::getValues_HDIV_Args) rank = 3 required for outputValues when operator is VALUE.");
680
681 INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.extent(2) == spaceDim ),
682 std::invalid_argument,
683 ">>> ERROR: (Intrepid2::getValues_HDIV_Args) dim 2 of outputValues must equal cell dimension for operator VALUE.");
684 break;
685 case OPERATOR_DIV:
686 INTREPID2_TEST_FOR_EXCEPTION( !(rank(outputValues) == 2), std::invalid_argument,
687 ">>> ERROR: (Intrepid2::getValues_HDIV_Args) rank = 2 required for outputValues when operator is DIV.");
688 break;
689
690 default:
691 INTREPID2_TEST_FOR_EXCEPTION( (true), std::invalid_argument, ">>> ERROR: (Intrepid2::getValues_HDIV_Args) Invalid operator");
692 }
693
694
695 // Verify dim 0 and dim 1 of outputValues
696 INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.extent(1) == inputPoints.extent(0) ),
697 std::invalid_argument,
698 ">>> ERROR: (Intrepid2::getValues_HDIV_Args) dim 1 (number of points) of outputValues must equal dim 0 of inputPoints.");
699
700 INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.extent(0) == static_cast<size_type>(basisCard) ),
701 std::invalid_argument,
702 ">>> ERROR: (Intrepid2::getValues_HDIV_Args) dim 0 (number of basis functions) of outputValues must equal basis cardinality.");
703 }
704
705 template<typename outputValueViewType,
706 typename inputPointViewType>
707 void getValues_HVOL_Args( const outputValueViewType outputValues,
708 const inputPointViewType inputPoints,
709 const EOperator operatorType,
710 const shards::CellTopology cellTopo,
711 const ordinal_type basisCard ) {
712 const auto spaceDim = cellTopo.getDimension();
713
714 // Verify inputPoints array
715 INTREPID2_TEST_FOR_EXCEPTION( !(rank(inputPoints) == 2), std::invalid_argument,
716 ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) rank = 2 required for inputPoints array");
717
718 INTREPID2_TEST_FOR_EXCEPTION( (inputPoints.extent(0) <= 0), std::invalid_argument,
719 ">>> ERROR (Intrepid2::getValues_HGRAD_Args): dim 0 (number of points) > 0 required for inputPoints array");
720
721 INTREPID2_TEST_FOR_EXCEPTION( !(inputPoints.extent(1) == spaceDim), std::invalid_argument,
722 ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) dim 1 (spatial dimension) of inputPoints array does not match cell dimension");
723
724
725 // Verify that all inputPoints are in the reference cell
726 /*
727 INTREPID2_TEST_FOR_EXCEPTION( !CellTools<Scalar>::checkPointSetInclusion(inputPoints, cellTopo), std::invalid_argument,
728 ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) One or more points are outside the "
729 << cellTopo <<" reference cell");
730 */
731
732
733 // Verify that operatorType is admissible for HGRAD fields
734 INTREPID2_TEST_FOR_EXCEPTION( ( (spaceDim == 2) && (operatorType == OPERATOR_DIV) ), std::invalid_argument,
735 ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) DIV is invalid operator for rank-0 (scalar) fields in 2D.");
736
737 INTREPID2_TEST_FOR_EXCEPTION( ( (spaceDim == 3) && ( (operatorType == OPERATOR_DIV) ||
738 (operatorType == OPERATOR_CURL) ) ), std::invalid_argument,
739 ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) DIV and CURL are invalid operators for rank-0 (scalar) fields in 3D.");
740
741
742 // Check rank of outputValues (all operators are admissible in 1D) and its dim 2 when operator is
743 // GRAD, CURL (only in 2D), or Dk.
744
745 if(spaceDim == 1) {
746 switch(operatorType){
747 case OPERATOR_VALUE:
748 INTREPID2_TEST_FOR_EXCEPTION( !(rank(outputValues) == 2), std::invalid_argument,
749 ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) rank = 2 required for outputValues when operator = VALUE.");
750 break;
751 case OPERATOR_GRAD:
752 case OPERATOR_CURL:
753 case OPERATOR_DIV:
754 case OPERATOR_D1:
755 case OPERATOR_D2:
756 case OPERATOR_D3:
757 case OPERATOR_D4:
758 case OPERATOR_D5:
759 case OPERATOR_D6:
760 case OPERATOR_D7:
761 case OPERATOR_D8:
762 case OPERATOR_D9:
763 case OPERATOR_D10:
764 INTREPID2_TEST_FOR_EXCEPTION( !(rank(outputValues) == 3), std::invalid_argument,
765 ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) rank = 3 required for outputValues in 1D when operator = GRAD, CURL, DIV, or Dk.");
766
767 INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.extent(2) == 1 ),
768 std::invalid_argument,
769 ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) dim 2 of outputValues must equal 1 when operator = GRAD, CURL, DIV, or Dk.");
770 break;
771 default:
772 INTREPID2_TEST_FOR_EXCEPTION( (true), std::invalid_argument, ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) Invalid operator");
773 }
774 }
775 else if(spaceDim > 1) {
776 switch(operatorType){
777 case OPERATOR_VALUE:
778 INTREPID2_TEST_FOR_EXCEPTION( !(rank(outputValues) == 2), std::invalid_argument,
779 ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) rank = 2 required for outputValues when operator = VALUE.");
780 break;
781 case OPERATOR_GRAD:
782 case OPERATOR_CURL:
783 case OPERATOR_D1:
784 INTREPID2_TEST_FOR_EXCEPTION( !(rank(outputValues) == 3), std::invalid_argument,
785 ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) rank = 3 required for outputValues in 2D and 3D when operator = GRAD, CURL (in 2D), or Dk.");
786
787 INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.extent(2) == spaceDim ),
788 std::invalid_argument,
789 ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) dim 2 of outputValues must equal cell dimension when operator = GRAD, CURL (in 2D), or D1.");
790 break;
791 case OPERATOR_D2:
792 case OPERATOR_D3:
793 case OPERATOR_D4:
794 case OPERATOR_D5:
795 case OPERATOR_D6:
796 case OPERATOR_D7:
797 case OPERATOR_D8:
798 case OPERATOR_D9:
799 case OPERATOR_D10:
800 INTREPID2_TEST_FOR_EXCEPTION( !(rank(outputValues) == 3), std::invalid_argument,
801 ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) rank = 3 required for outputValues in 2D and 3D when operator = GRAD, CURL (in 2D), or Dk.");
802
803 INTREPID2_TEST_FOR_EXCEPTION( !(static_cast<ordinal_type>(outputValues.extent(2)) == getDkCardinality(operatorType, spaceDim)),
804 std::invalid_argument,
805 ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) dim 2 of outputValues must equal cardinality of the Dk multiset.");
806 break;
807 default:
808 INTREPID2_TEST_FOR_EXCEPTION( (true), std::invalid_argument, ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) Invalid operator");
809 }
810 }
811
812
813 // Verify dim 0 and dim 1 of outputValues
814 INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.extent(1) == inputPoints.extent(0) ),
815 std::invalid_argument,
816 ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) dim 1 (number of points) of outputValues must equal dim 0 of inputPoints.");
817
818 INTREPID2_TEST_FOR_EXCEPTION( !(static_cast<ordinal_type>(outputValues.extent(0)) == basisCard ),
819 std::invalid_argument,
820 ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) dim 0 (number of basis functions) of outputValues must equal basis cardinality.");
821 }
822
823 template<typename Device,
824 typename outputValueType,
825 typename pointValueType>
826 Kokkos::DynRankView<outputValueType,Device>
827 Basis<Device,outputValueType,pointValueType>::allocateOutputView( const int numPoints, const EOperator operatorType) const
828 {
829 const bool operatorIsDk = (operatorType >= OPERATOR_D1) && (operatorType <= OPERATOR_D10);
830 const bool operatorSupported = (operatorType == OPERATOR_VALUE) || (operatorType == OPERATOR_GRAD) || (operatorType == OPERATOR_CURL) || (operatorType == OPERATOR_DIV) || operatorIsDk;
831 INTREPID2_TEST_FOR_EXCEPTION(!operatorSupported, std::invalid_argument, "operator is not supported by allocateOutputView()");
832
833 const int numFields = this->getCardinality();
834 const int spaceDim = this->getBaseCellTopology().getDimension() + this->getNumTensorialExtrusions();
835
836 using OutputViewAllocatable = Kokkos::DynRankView<outputValueType,DeviceType>;
837
838 switch (functionSpace_)
839 {
840 case FUNCTION_SPACE_HGRAD:
841 if (operatorType == OPERATOR_VALUE)
842 {
843 // scalar-valued container
844 OutputViewAllocatable dataView("BasisValues HGRAD VALUE data", numFields, numPoints);
845 return dataView;
846 }
847 else if (operatorType == OPERATOR_GRAD)
848 {
849 OutputViewAllocatable dataView("BasisValues HGRAD GRAD data", numFields, numPoints, spaceDim);
850 return dataView;
851 }
852 else if (operatorIsDk)
853 {
854 ordinal_type dkCardinality = getDkCardinality(operatorType, spaceDim);
855 OutputViewAllocatable dataView("BasisValues HGRAD Dk data", numFields, numPoints, dkCardinality);
856 return dataView;
857 }
858 else
859 {
860 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "operator/space combination not supported by allocateOutputView()");
861 }
862 case FUNCTION_SPACE_HDIV:
863 if (operatorType == OPERATOR_VALUE)
864 {
865 // vector-valued container
866 OutputViewAllocatable dataView("BasisValues HDIV VALUE data", numFields, numPoints, spaceDim);
867 return dataView;
868 }
869 else if (operatorType == OPERATOR_DIV)
870 {
871 // scalar-valued curl
872 OutputViewAllocatable dataView("BasisValues HDIV DIV data", numFields, numPoints);
873 return dataView;
874 }
875 else
876 {
877 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "operator/space combination not supported by allocateOutputView()");
878 }
879 case FUNCTION_SPACE_HCURL:
880 if (operatorType == OPERATOR_VALUE)
881 {
882 OutputViewAllocatable dataView("BasisValues HCURL VALUE data", numFields, numPoints, spaceDim);
883 return dataView;
884 }
885 else if (operatorType == OPERATOR_CURL)
886 {
887 if (spaceDim != 2)
888 {
889 // vector-valued curl
890 OutputViewAllocatable dataView("BasisValues HCURL CURL data", numFields, numPoints, spaceDim);
891 return dataView;
892 }
893 else
894 {
895 // scalar-valued curl
896 OutputViewAllocatable dataView("BasisValues HCURL CURL data (scalar)", numFields, numPoints);
897 return dataView;
898 }
899 }
900 else
901 {
902 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "operator/space combination not supported by allocateOutputView()");
903 }
904 case FUNCTION_SPACE_HVOL:
905 if (operatorType == OPERATOR_VALUE)
906 {
907 // vector-valued container
908 OutputViewAllocatable dataView("BasisValues HVOL VALUE data", numFields, numPoints);
909 return dataView;
910 }
911 else if (operatorIsDk || (operatorType == OPERATOR_GRAD))
912 {
913 ordinal_type dkCardinality = getDkCardinality(operatorType, spaceDim);
914 OutputViewAllocatable dataView("BasisValues HVOL Dk data", numFields, numPoints, dkCardinality);
915 return dataView;
916 }
917 else
918 {
919 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "operator/space combination not supported by allocateOutputView()");
920 }
921 default:
922 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "operator/space combination not supported by allocateOutputView()");
923 }
924 }
925}
926
927#endif
KOKKOS_INLINE_FUNCTION void getJacobyRecurrenceCoeffs(value_type &an, value_type &bn, value_type &cn, const ordinal_type alpha, const ordinal_type beta, const ordinal_type n)
function for computing the Jacobi recurrence coefficients so that
KOKKOS_INLINE_FUNCTION ordinal_type getOperatorRank(const EFunctionSpace spaceType, const EOperator operatorType, const ordinal_type spaceDim)
Returns rank of an operator.
KOKKOS_INLINE_FUNCTION ordinal_type getPnEnumeration(const ordinal_type p, const ordinal_type q=0, const ordinal_type r=0)
Returns the index of the term x^p y^q z^r of a polynomial of degree n (p+q+r <= n)....
KOKKOS_INLINE_FUNCTION ordinal_type getFieldRank(const EFunctionSpace spaceType)
Returns the rank of fields in a function space of the specified type.
KOKKOS_INLINE_FUNCTION ordinal_type getOperatorOrder(const EOperator operatorType)
Returns order of an operator.
void getValues_HGRAD_Args(const outputValueViewType outputValues, const inputPointViewType inputPoints, const EOperator operatorType, const shards::CellTopology cellTopo, const ordinal_type basisCard)
Runtime check of the arguments for the getValues method in an HGRAD-conforming FEM basis....
KOKKOS_INLINE_FUNCTION ordinal_type getPnCardinality(ordinal_type n)
Returns cardinality of Polynomials of order n (P^n).
KOKKOS_INLINE_FUNCTION ordinal_type getDkCardinality(const EOperator operatorType, const ordinal_type spaceDim)
Returns cardinality of Dk, i.e., the number of all derivatives of order k.
void getValues_HVOL_Args(const outputValueViewType outputValues, const inputPointViewType inputPoints, const EOperator operatorType, const shards::CellTopology cellTopo, const ordinal_type basisCard)
Runtime check of the arguments for the getValues method in an HVOL-conforming FEM basis....
void getValues_HDIV_Args(const outputValueViewType outputValues, const inputPointViewType inputPoints, const EOperator operatorType, const shards::CellTopology cellTopo, const ordinal_type basisCard)
Runtime check of the arguments for the getValues method in an HDIV-conforming FEM basis....
void getValues_HCURL_Args(const outputValueViewType outputValues, const inputPointViewType inputPoints, const EOperator operatorType, const shards::CellTopology cellTopo, const ordinal_type basisCard)
Runtime check of the arguments for the getValues method in an HCURL-conforming FEM basis....
KOKKOS_INLINE_FUNCTION ordinal_type getDkEnumeration(const ordinal_type xMult, const ordinal_type yMult=-1, const ordinal_type zMult=-1)
Returns the ordinal of a partial derivative of order k based on the multiplicities of the partials dx...
KOKKOS_FORCEINLINE_FUNCTION bool isValidFunctionSpace(const EFunctionSpace spaceType)
Verifies validity of a function space enum.
KOKKOS_FORCEINLINE_FUNCTION bool isValidOperator(const EOperator operatorType)
Verifies validity of an operator enum.
Kokkos::DynRankView< OutputValueType, DeviceType > allocateOutputView(const int numPoints, const EOperator operatorType=OPERATOR_VALUE) const
Allocate a View container suitable for passing to the getValues() variant that accepts Kokkos DynRank...