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