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 
21 namespace 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
237  void getJacobyRecurrenceCoeffs (
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
Intrepid2::Basis::allocateOutputView
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...
Definition: Intrepid2_BasisDef.hpp:827
Generated by   doxygen 1.8.5