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:
246  return 0;
247 #ifndef __NVCC__ //prevent nvcc warning
248  break;
249 #endif
250 
251  case 2:
252  return yMult;
253 #ifndef __NVCC__ //prevent nvcc warning
254  break;
255 #endif
256 
257  case 3:
258  return zMult + (yMult+zMult)*(yMult+zMult+1)/2;
259 #ifndef __NVCC__ //prevent nvcc warning
260  break;
261 #endif
262 
263  default: {
264  INTREPID2_TEST_FOR_ABORT( !( (0 < spaceDim ) && (spaceDim < 4) ),
265  ">>> ERROR (Intrepid2::getDkEnumeration): Invalid space dimension");
266  return 0;
267  }
268  }
269  }
270 
271  template<ordinal_type spaceDim>
272  KOKKOS_INLINE_FUNCTION
273  ordinal_type getPnEnumeration(const ordinal_type p,
274  const ordinal_type q /* = 0*/,
275  const ordinal_type r /* = 0*/) {
276  return (spaceDim==1) ? p :
277  (spaceDim==2) ? (p+q)*(p+q+1)/2+q :
278  (p+q+r)*(p+q+r+1)*(p+q+r+2)/6+(q+r)*(q+r+1)/2+r;
279 
280  }
281 
282  template<typename value_type>
283  KOKKOS_INLINE_FUNCTION
284  void getJacobyRecurrenceCoeffs (
285  value_type &an,
286  value_type &bn,
287  value_type &cn,
288  const ordinal_type alpha,
289  const ordinal_type beta ,
290  const ordinal_type n) {
291  an = ( (2.0 * n + 1.0 + alpha + beta) * ( 2.0 * n + 2.0 + alpha + beta ) /
292  value_type(2.0 * ( n + 1 ) * ( n + 1 + alpha + beta ) ) );
293  bn = ( (alpha*alpha-beta*beta)*(2.0*n+1.0+alpha+beta) /
294  value_type(2.0*(n+1.0)*(2.0*n+alpha+beta)*(n+1.0+alpha+beta) ) );
295  cn = ( (n+alpha)*(n+beta)*(2.0*n+2.0+alpha+beta) /
296  value_type( (n+1.0)*(n+1.0+alpha+beta)*(2.0*n+alpha+beta) ) );
297  }
298 
299 
300 // template<typename OrdinalArrayType>
301 // KOKKOS_INLINE_FUNCTION
302 // void getDkMultiplicities(OrdinalArrayType partialMult,
303 // const ordinal_type derivativeEnum,
304 // const EOperator operatorType,
305 // const ordinal_type spaceDim) {
306 
307 // /* Hash table to convert enumeration of partial derivative to multiplicities of dx,dy,dz in 3D.
308 // Multiplicities {mx,my,mz} are arranged lexicographically in bins numbered from 0 to 10.
309 // The size of bins is an arithmetic progression, i.e., 1,2,3,4,5,...,11. Conversion formula is:
310 // \verbatim
311 // mx = derivativeOrder - binNumber
312 // mz = derivativeEnum - binBegin
313 // my = derivativeOrder - mx - mz = binNumber + binBegin - derivativeEnum
314 // \endverbatim
315 // where binBegin is the enumeration of the first element in the bin. Bin numbers and binBegin
316 // values are stored in hash tables for quick access by derivative enumeration value.
317 // */
318 
319 // // Returns the bin number for the specified derivative enumeration
320 // const ordinal_type binNumber[66] = {
321 // 0,
322 // 1, 1,
323 // 2, 2, 2,
324 // 3, 3, 3, 3,
325 // 4, 4, 4, 4, 4,
326 // 5, 5, 5, 5, 5, 5,
327 // 6, 6, 6, 6, 6, 6, 6,
328 // 7, 7, 7, 7, 7, 7, 7, 7,
329 // 8, 8, 8, 8, 8, 8, 8, 8, 8,
330 // 9, 9, 9, 9, 9, 9, 9, 9, 9, 9,
331 // 10,10,10,10,10,10,10,10,10,10,10
332 // };
333 
334 // // Returns the binBegin value for the specified derivative enumeration
335 // const ordinal_type binBegin[66] = {
336 // 0,
337 // 1, 1,
338 // 3, 3 ,3,
339 // 6, 6, 6, 6,
340 // 10,10,10,10,10,
341 // 15,15,15,15,15,15,
342 // 21,21,21,21,21,21,21,
343 // 28,28,28,28,28,28,28,28,
344 // 36,36,36,36,36,36,36,36,36,
345 // 45,45,45,45,45,45,45,45,45,45,
346 // 55,55,55,55,55,55,55,55,55,55,55
347 // };
348 
349 // #ifdef HAVE_INTREPID2_DEBUG
350 // // Enumeration value must be between 0 and the cardinality of the derivative set
351 // INTREPID2_TEST_FOR_ABORT( !( (0 <= derivativeEnum) && (derivativeEnum < getDkCardinality(operatorType,spaceDim) ) ),
352 // ">>> ERROR (Intrepid2::getDkMultiplicities): Invalid derivative enumeration value for this order and space dimension");
353 // #endif
354 
355 // // This method should only be called for Dk operators
356 // ordinal_type derivativeOrder;
357 // switch(operatorType) {
358 
359 // case OPERATOR_D1:
360 // case OPERATOR_D2:
361 // case OPERATOR_D3:
362 // case OPERATOR_D4:
363 // case OPERATOR_D5:
364 // case OPERATOR_D6:
365 // case OPERATOR_D7:
366 // case OPERATOR_D8:
367 // case OPERATOR_D9:
368 // case OPERATOR_D10:
369 // derivativeOrder = Intrepid2::getOperatorOrder(operatorType);
370 // break;
371 
372 // default:
373 // INTREPID2_TEST_FOR_ABORT(true,
374 // ">>> ERROR (Intrepid2::getDkMultiplicities): operator type Dk required for this method");
375 // }// switch
376 
377 // #ifdef HAVE_INTREPID2_DEBUG
378 // INTREPID2_TEST_FOR_ABORT( partialMult.extent(0) != spaceDim,
379 // ">>> ERROR (Intrepid2::getDkMultiplicities): partialMult must have the same dimension as spaceDim" );
380 // #endif
381 
382 // switch(spaceDim) {
383 
384 // case 1:
385 // // Multiplicity of dx equals derivativeOrder
386 // partialMult(0) = derivativeOrder;
387 // break;
388 
389 // case 2:
390 // // Multiplicity of dy equals the enumeration of the derivative; of dx - the complement
391 // partialMult(1) = derivativeEnum;
392 // partialMult(0) = derivativeOrder - derivativeEnum;
393 // break;
394 
395 // case 3:
396 // // Recover multiplicities
397 // partialMult(0) = derivativeOrder - binNumber[derivativeEnum];
398 // partialMult(1) = binNumber[derivativeEnum] + binBegin[derivativeEnum] - derivativeEnum;
399 // partialMult(2) = derivativeEnum - binBegin[derivativeEnum];
400 // break;
401 
402 // default:
403 // INTREPID2_TEST_FOR_ABORT( !( (0 < spaceDim ) && (spaceDim < 4) ),
404 // ">>> ERROR (Intrepid2::getDkMultiplicities): Invalid space dimension");
405 // }
406 // }
407 
408 
409  KOKKOS_INLINE_FUNCTION
410  ordinal_type getDkCardinality(const EOperator operatorType,
411  const ordinal_type spaceDim) {
412 
413 #ifdef HAVE_INTREPID2_DEBUG
414  INTREPID2_TEST_FOR_ABORT( !( (0 < spaceDim ) && (spaceDim < 4) ),
415  ">>> ERROR (Intrepid2::getDkcardinality): Invalid space dimension");
416  switch (operatorType) {
417  case OPERATOR_GRAD:
418  case OPERATOR_D1:
419  case OPERATOR_D2:
420  case OPERATOR_D3:
421  case OPERATOR_D4:
422  case OPERATOR_D5:
423  case OPERATOR_D6:
424  case OPERATOR_D7:
425  case OPERATOR_D8:
426  case OPERATOR_D9:
427  case OPERATOR_D10:
428  break;
429  default:
430  INTREPID2_TEST_FOR_ABORT( true, ">>> ERROR (Intrepid2::getDkCardinality): Cannot be used for this operator ");
431  break;
432  }
433 #endif
434 
435  ordinal_type n = Intrepid2::getOperatorOrder(operatorType);
436  return (spaceDim==1) ? 1 :
437  (spaceDim==2) ? n+1 :
438  (n + 1) * (n + 2) / 2;
439  }
440 
441  template<EOperator operatorType, ordinal_type spaceDim>
442  KOKKOS_INLINE_FUNCTION
443  constexpr ordinal_type getDkCardinality() {
444  return getPnCardinality<spaceDim-1,Intrepid2::getOperatorOrder<operatorType>()>();
445  }
446 
447  template<ordinal_type spaceDim>
448  KOKKOS_INLINE_FUNCTION
449  ordinal_type getPnCardinality (ordinal_type n) {
450 
451 #ifdef HAVE_INTREPID2_DEBUG
452  INTREPID2_TEST_FOR_ABORT( !( (0 <= spaceDim ) && (spaceDim < 4) ),
453  ">>> ERROR (Intrepid2::getPnCardinality): Invalid space dimension");
454 #endif
455 
456  return (spaceDim==0) ? 1 :
457  (spaceDim==1) ? n+1 :
458  (spaceDim==2) ? (n + 1) * (n + 2) / 2 :
459  (n + 1) * (n + 2) * (n + 3) / 6;
460  }
461 
462 
463  template<ordinal_type spaceDim, ordinal_type n>
464  KOKKOS_INLINE_FUNCTION
465  constexpr ordinal_type getPnCardinality () {
466 
467  return (spaceDim==0) ? 1 :
468  (spaceDim==1) ? n+1 :
469  (spaceDim==2) ? (n + 1) * (n + 2) / 2 :
470  (n + 1) * (n + 2) * (n + 3) / 6;
471 
472  }
473 
474 
475 
476 
477  //
478  // Argument check
479  //
480 
481 
482  template<typename outputValueViewType,
483  typename inputPointViewType>
484  void getValues_HGRAD_Args( const outputValueViewType outputValues,
485  const inputPointViewType inputPoints,
486  const EOperator operatorType,
487  const shards::CellTopology cellTopo,
488  const ordinal_type basisCard ) {
489  const auto spaceDim = cellTopo.getDimension();
490 
491  // Verify inputPoints array
492  INTREPID2_TEST_FOR_EXCEPTION( !(inputPoints.rank() == 2), std::invalid_argument,
493  ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) rank = 2 required for inputPoints array");
494 
495  INTREPID2_TEST_FOR_EXCEPTION( (inputPoints.extent(0) <= 0), std::invalid_argument,
496  ">>> ERROR (Intrepid2::getValues_HGRAD_Args): dim 0 (number of points) > 0 required for inputPoints array");
497 
498  INTREPID2_TEST_FOR_EXCEPTION( !(inputPoints.extent(1) == spaceDim), std::invalid_argument,
499  ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) dim 1 (spatial dimension) of inputPoints array does not match cell dimension");
500 
501 
502  // Verify that all inputPoints are in the reference cell
503  /*
504  INTREPID2_TEST_FOR_EXCEPTION( !CellTools<Scalar>::checkPointSetInclusion(inputPoints, cellTopo), std::invalid_argument,
505  ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) One or more points are outside the "
506  << cellTopo <<" reference cell");
507  */
508 
509 
510  // Verify that operatorType is admissible for HGRAD fields
511  INTREPID2_TEST_FOR_EXCEPTION( ( (spaceDim == 2) && (operatorType == OPERATOR_DIV) ), std::invalid_argument,
512  ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) DIV is invalid operator for rank-0 (scalar) fields in 2D.");
513 
514  INTREPID2_TEST_FOR_EXCEPTION( ( (spaceDim == 3) && ( (operatorType == OPERATOR_DIV) ||
515  (operatorType == OPERATOR_CURL) ) ), std::invalid_argument,
516  ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) DIV and CURL are invalid operators for rank-0 (scalar) fields in 3D.");
517 
518 
519  // Check rank of outputValues (all operators are admissible in 1D) and its dim 2 when operator is
520  // GRAD, CURL (only in 2D), or Dk.
521 
522  if(spaceDim == 1) {
523  switch(operatorType){
524  case OPERATOR_VALUE:
525  INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.rank() == 2), std::invalid_argument,
526  ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) rank = 2 required for outputValues when operator = VALUE.");
527  break;
528  case OPERATOR_GRAD:
529  case OPERATOR_CURL:
530  case OPERATOR_DIV:
531  case OPERATOR_D1:
532  case OPERATOR_D2:
533  case OPERATOR_D3:
534  case OPERATOR_D4:
535  case OPERATOR_D5:
536  case OPERATOR_D6:
537  case OPERATOR_D7:
538  case OPERATOR_D8:
539  case OPERATOR_D9:
540  case OPERATOR_D10:
541  INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.rank() == 3), std::invalid_argument,
542  ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) rank = 3 required for outputValues in 1D when operator = GRAD, CURL, DIV, or Dk.");
543 
544  INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.extent(2) == 1 ),
545  std::invalid_argument,
546  ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) dim 2 of outputValues must equal 1 when operator = GRAD, CURL, DIV, or Dk.");
547  break;
548  default:
549  INTREPID2_TEST_FOR_EXCEPTION( (true), std::invalid_argument, ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) Invalid operator");
550  }
551  }
552  else if(spaceDim > 1) {
553  switch(operatorType){
554  case OPERATOR_VALUE:
555  INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.rank() == 2), std::invalid_argument,
556  ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) rank = 2 required for outputValues when operator = VALUE.");
557  break;
558  case OPERATOR_GRAD:
559  case OPERATOR_CURL:
560  case OPERATOR_D1:
561  INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.rank() == 3), std::invalid_argument,
562  ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) rank = 3 required for outputValues in 2D and 3D when operator = GRAD, CURL (in 2D), or Dk.");
563 
564  INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.extent(2) == spaceDim ),
565  std::invalid_argument,
566  ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) dim 2 of outputValues must equal cell dimension when operator = GRAD, CURL (in 2D), or D1.");
567  break;
568  case OPERATOR_D2:
569  case OPERATOR_D3:
570  case OPERATOR_D4:
571  case OPERATOR_D5:
572  case OPERATOR_D6:
573  case OPERATOR_D7:
574  case OPERATOR_D8:
575  case OPERATOR_D9:
576  case OPERATOR_D10:
577  INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.rank() == 3), std::invalid_argument,
578  ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) rank = 3 required for outputValues in 2D and 3D when operator = GRAD, CURL (in 2D), or Dk.");
579 
580  INTREPID2_TEST_FOR_EXCEPTION( !(static_cast<ordinal_type>(outputValues.extent(2)) == getDkCardinality(operatorType, spaceDim)),
581  std::invalid_argument,
582  ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) dim 2 of outputValues must equal cardinality of the Dk multiset.");
583  break;
584  default:
585  INTREPID2_TEST_FOR_EXCEPTION( (true), std::invalid_argument, ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) Invalid operator");
586  }
587  }
588 
589 
590  // Verify dim 0 and dim 1 of outputValues
591  INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.extent(1) == inputPoints.extent(0) ),
592  std::invalid_argument,
593  ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) dim 1 (number of points) of outputValues must equal dim 0 of inputPoints.");
594 
595  INTREPID2_TEST_FOR_EXCEPTION( !(static_cast<ordinal_type>(outputValues.extent(0)) == basisCard ),
596  std::invalid_argument,
597  ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) dim 0 (number of basis functions) of outputValues must equal basis cardinality.");
598  }
599 
600 
601  template<typename outputValueViewType,
602  typename inputPointViewType>
603  void getValues_HCURL_Args( const outputValueViewType outputValues,
604  const inputPointViewType inputPoints,
605  const EOperator operatorType,
606  const shards::CellTopology cellTopo,
607  const ordinal_type basisCard ) {
608 
609  const auto spaceDim = cellTopo.getDimension();
610 
611  // Verify that cell is 2D or 3D (this is redundant for default bases where we use correct cells,
612  // but will force any user-defined basis for HCURL spaces to use 2D or 3D cells
613  INTREPID2_TEST_FOR_EXCEPTION( !( (spaceDim == 2) || (spaceDim == 3) ), std::invalid_argument,
614  ">>> ERROR: (Intrepid2::getValues_HCURL_Args) cell dimension = 2 or 3 required for HCURL basis");
615 
616 
617  // Verify inputPoints array
618  INTREPID2_TEST_FOR_EXCEPTION( !(inputPoints.rank() == 2), std::invalid_argument,
619  ">>> ERROR: (Intrepid2::getValues_HCURL_Args) rank = 2 required for inputPoints array");
620  INTREPID2_TEST_FOR_EXCEPTION( (inputPoints.extent(0) <= 0), std::invalid_argument,
621  ">>> ERROR (Intrepid2::getValues_HCURL_Args): dim 0 (number of points) > 0 required for inputPoints array");
622 
623  INTREPID2_TEST_FOR_EXCEPTION( !(inputPoints.extent(1) == spaceDim), std::invalid_argument,
624  ">>> ERROR: (Intrepid2::getValues_HCURL_Args) dim 1 (spatial dimension) of inputPoints array does not match cell dimension");
625 
626  // Verify that all inputPoints are in the reference cell
627  /*
628  INTREPID2_TEST_FOR_EXCEPTION( !CellTools<Scalar>::checkPointSetInclusion(inputPoints, cellTopo), std::invalid_argument,
629  ">>> ERROR: (Intrepid2::getValues_HCURL_Args) One or more points are outside the "
630  << cellTopo <<" reference cell");
631  */
632 
633  // Verify that operatorType is admissible for HCURL fields
634  INTREPID2_TEST_FOR_EXCEPTION( !( (operatorType == OPERATOR_VALUE) || (operatorType == OPERATOR_CURL) ), std::invalid_argument,
635  ">>> ERROR: (Intrepid2::getValues_HCURL_Args) operator = VALUE or CURL required for HCURL fields.");
636 
637 
638  // Check rank of outputValues: for VALUE should always be rank-3 array with (F,P,D) layout
639  switch(operatorType) {
640 
641  case OPERATOR_VALUE:
642  INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.rank() == 3), std::invalid_argument,
643  ">>> ERROR: (Intrepid2::getValues_HCURL_Args) rank = 3 required for outputValues when operator is VALUE");
644  INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.extent(2) == spaceDim ),
645  std::invalid_argument,
646  ">>> ERROR: (Intrepid2::getValues_HCURL_Args) dim 2 of outputValues must equal cell dimension when operator is VALUE.");
647  break;
648 
649  case OPERATOR_CURL:
650 
651  // in 3D we need an (F,P,D) container because CURL gives a vector field:
652  if(spaceDim == 3) {
653  INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.rank() == 3 ) ,
654  std::invalid_argument,
655  ">>> ERROR: (Intrepid2::getValues_HCURL_Args) rank = 3 required for outputValues in 3D when operator is CURL");
656  INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.extent(2) == spaceDim),
657  std::invalid_argument,
658  ">>> ERROR: (Intrepid2::getValues_HCURL_Args) dim 2 of outputValues must equal cell dimension in 3D when operator is CURL.");
659  }
660  // In 2D we need an (F,P) container because CURL gives a scalar field
661  else if(spaceDim == 2) {
662  INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.rank() == 2 ) ,
663  std::invalid_argument,
664  ">>> ERROR: (Intrepid2::getValues_HCURL_Args) rank = 2 required for outputValues in 2D when operator is CURL");
665  }
666  break;
667 
668  default:
669  INTREPID2_TEST_FOR_EXCEPTION( (true), std::invalid_argument, ">>> ERROR: (Intrepid2::getValues_HCURL_Args) Invalid operator");
670  }
671 
672 
673  // Verify dim 0 and dim 1 of outputValues
674  INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.extent(1) == inputPoints.extent(0) ),
675  std::invalid_argument,
676  ">>> ERROR: (Intrepid2::getValues_HCURL_Args) dim 1 (number of points) of outputValues must equal dim 0 of inputPoints.");
677 
678  INTREPID2_TEST_FOR_EXCEPTION( !(static_cast<ordinal_type>(outputValues.extent(0)) == basisCard ),
679  std::invalid_argument,
680  ">>> ERROR: (Intrepid2::getValues_HCURL_Args) dim 0 (number of basis functions) of outputValues must equal basis cardinality.");
681 
682  }
683 
684 
685 
686  template<typename outputValueViewType,
687  typename inputPointViewType>
688  void getValues_HDIV_Args( const outputValueViewType outputValues,
689  const inputPointViewType inputPoints,
690  const EOperator operatorType,
691  const shards::CellTopology cellTopo,
692  const ordinal_type basisCard ) {
693 
694  const auto spaceDim = cellTopo.getDimension();
695 
696  // Verify inputPoints array
697  INTREPID2_TEST_FOR_EXCEPTION( !(inputPoints.rank() == 2), std::invalid_argument,
698  ">>> ERROR: (Intrepid2::getValues_HDIV_Args) rank = 2 required for inputPoints array");
699  INTREPID2_TEST_FOR_EXCEPTION( (inputPoints.extent(0) <= 0), std::invalid_argument,
700  ">>> ERROR (Intrepid2::getValues_HDIV_Args): dim 0 (number of points) > 0 required for inputPoints array");
701 
702  INTREPID2_TEST_FOR_EXCEPTION( !(inputPoints.extent(1) == spaceDim), std::invalid_argument,
703  ">>> ERROR: (Intrepid2::getValues_HDIV_Args) dim 1 (spatial dimension) of inputPoints array does not match cell dimension");
704 
705  // Verify that all inputPoints are in the reference cell
706  /*
707  INTREPID2_TEST_FOR_EXCEPTION( !CellTools<Scalar>::checkPointSetInclusion(inputPoints, cellTopo), std::invalid_argument,
708  ">>> ERROR: (Intrepid2::getValues_HDIV_Args) One or more points are outside the "
709  << cellTopo <<" reference cell");
710  */
711 
712  // Verify that operatorType is admissible for HDIV fields
713  INTREPID2_TEST_FOR_EXCEPTION( !( (operatorType == OPERATOR_VALUE) || (operatorType == OPERATOR_DIV) ), std::invalid_argument,
714  ">>> ERROR: (Intrepid2::getValues_HDIV_Args) operator = VALUE or DIV required for HDIV fields.");
715 
716 
717  // Check rank of outputValues
718  switch(operatorType) {
719  case OPERATOR_VALUE:
720  INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.rank() == 3), std::invalid_argument,
721  ">>> ERROR: (Intrepid2::getValues_HDIV_Args) rank = 3 required for outputValues when operator is VALUE.");
722 
723  INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.extent(2) == spaceDim ),
724  std::invalid_argument,
725  ">>> ERROR: (Intrepid2::getValues_HDIV_Args) dim 2 of outputValues must equal cell dimension for operator VALUE.");
726  break;
727  case OPERATOR_DIV:
728  INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.rank() == 2), std::invalid_argument,
729  ">>> ERROR: (Intrepid2::getValues_HDIV_Args) rank = 2 required for outputValues when operator is DIV.");
730  break;
731 
732  default:
733  INTREPID2_TEST_FOR_EXCEPTION( (true), std::invalid_argument, ">>> ERROR: (Intrepid2::getValues_HDIV_Args) Invalid operator");
734  }
735 
736 
737  // Verify dim 0 and dim 1 of outputValues
738  INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.extent(1) == inputPoints.extent(0) ),
739  std::invalid_argument,
740  ">>> ERROR: (Intrepid2::getValues_HDIV_Args) dim 1 (number of points) of outputValues must equal dim 0 of inputPoints.");
741 
742  INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.extent(0) == static_cast<size_type>(basisCard) ),
743  std::invalid_argument,
744  ">>> ERROR: (Intrepid2::getValues_HDIV_Args) dim 0 (number of basis functions) of outputValues must equal basis cardinality.");
745  }
746 
747  template<typename outputValueViewType,
748  typename inputPointViewType>
749  void getValues_HVOL_Args( const outputValueViewType outputValues,
750  const inputPointViewType inputPoints,
751  const EOperator operatorType,
752  const shards::CellTopology cellTopo,
753  const ordinal_type basisCard ) {
754  const auto spaceDim = cellTopo.getDimension();
755 
756  // Verify inputPoints array
757  INTREPID2_TEST_FOR_EXCEPTION( !(inputPoints.rank() == 2), std::invalid_argument,
758  ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) rank = 2 required for inputPoints array");
759 
760  INTREPID2_TEST_FOR_EXCEPTION( (inputPoints.extent(0) <= 0), std::invalid_argument,
761  ">>> ERROR (Intrepid2::getValues_HGRAD_Args): dim 0 (number of points) > 0 required for inputPoints array");
762 
763  INTREPID2_TEST_FOR_EXCEPTION( !(inputPoints.extent(1) == spaceDim), std::invalid_argument,
764  ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) dim 1 (spatial dimension) of inputPoints array does not match cell dimension");
765 
766 
767  // Verify that all inputPoints are in the reference cell
768  /*
769  INTREPID2_TEST_FOR_EXCEPTION( !CellTools<Scalar>::checkPointSetInclusion(inputPoints, cellTopo), std::invalid_argument,
770  ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) One or more points are outside the "
771  << cellTopo <<" reference cell");
772  */
773 
774 
775  // Verify that operatorType is admissible for HGRAD fields
776  INTREPID2_TEST_FOR_EXCEPTION( ( (spaceDim == 2) && (operatorType == OPERATOR_DIV) ), std::invalid_argument,
777  ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) DIV is invalid operator for rank-0 (scalar) fields in 2D.");
778 
779  INTREPID2_TEST_FOR_EXCEPTION( ( (spaceDim == 3) && ( (operatorType == OPERATOR_DIV) ||
780  (operatorType == OPERATOR_CURL) ) ), std::invalid_argument,
781  ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) DIV and CURL are invalid operators for rank-0 (scalar) fields in 3D.");
782 
783 
784  // Check rank of outputValues (all operators are admissible in 1D) and its dim 2 when operator is
785  // GRAD, CURL (only in 2D), or Dk.
786 
787  if(spaceDim == 1) {
788  switch(operatorType){
789  case OPERATOR_VALUE:
790  INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.rank() == 2), std::invalid_argument,
791  ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) rank = 2 required for outputValues when operator = VALUE.");
792  break;
793  case OPERATOR_GRAD:
794  case OPERATOR_CURL:
795  case OPERATOR_DIV:
796  case OPERATOR_D1:
797  case OPERATOR_D2:
798  case OPERATOR_D3:
799  case OPERATOR_D4:
800  case OPERATOR_D5:
801  case OPERATOR_D6:
802  case OPERATOR_D7:
803  case OPERATOR_D8:
804  case OPERATOR_D9:
805  case OPERATOR_D10:
806  INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.rank() == 3), std::invalid_argument,
807  ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) rank = 3 required for outputValues in 1D when operator = GRAD, CURL, DIV, or Dk.");
808 
809  INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.extent(2) == 1 ),
810  std::invalid_argument,
811  ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) dim 2 of outputValues must equal 1 when operator = GRAD, CURL, DIV, or Dk.");
812  break;
813  default:
814  INTREPID2_TEST_FOR_EXCEPTION( (true), std::invalid_argument, ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) Invalid operator");
815  }
816  }
817  else if(spaceDim > 1) {
818  switch(operatorType){
819  case OPERATOR_VALUE:
820  INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.rank() == 2), std::invalid_argument,
821  ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) rank = 2 required for outputValues when operator = VALUE.");
822  break;
823  case OPERATOR_GRAD:
824  case OPERATOR_CURL:
825  case OPERATOR_D1:
826  INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.rank() == 3), std::invalid_argument,
827  ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) rank = 3 required for outputValues in 2D and 3D when operator = GRAD, CURL (in 2D), or Dk.");
828 
829  INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.extent(2) == spaceDim ),
830  std::invalid_argument,
831  ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) dim 2 of outputValues must equal cell dimension when operator = GRAD, CURL (in 2D), or D1.");
832  break;
833  case OPERATOR_D2:
834  case OPERATOR_D3:
835  case OPERATOR_D4:
836  case OPERATOR_D5:
837  case OPERATOR_D6:
838  case OPERATOR_D7:
839  case OPERATOR_D8:
840  case OPERATOR_D9:
841  case OPERATOR_D10:
842  INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.rank() == 3), std::invalid_argument,
843  ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) rank = 3 required for outputValues in 2D and 3D when operator = GRAD, CURL (in 2D), or Dk.");
844 
845  INTREPID2_TEST_FOR_EXCEPTION( !(static_cast<ordinal_type>(outputValues.extent(2)) == getDkCardinality(operatorType, spaceDim)),
846  std::invalid_argument,
847  ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) dim 2 of outputValues must equal cardinality of the Dk multiset.");
848  break;
849  default:
850  INTREPID2_TEST_FOR_EXCEPTION( (true), std::invalid_argument, ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) Invalid operator");
851  }
852  }
853 
854 
855  // Verify dim 0 and dim 1 of outputValues
856  INTREPID2_TEST_FOR_EXCEPTION( !(outputValues.extent(1) == inputPoints.extent(0) ),
857  std::invalid_argument,
858  ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) dim 1 (number of points) of outputValues must equal dim 0 of inputPoints.");
859 
860  INTREPID2_TEST_FOR_EXCEPTION( !(static_cast<ordinal_type>(outputValues.extent(0)) == basisCard ),
861  std::invalid_argument,
862  ">>> ERROR: (Intrepid2::getValues_HGRAD_Args) dim 0 (number of basis functions) of outputValues must equal basis cardinality.");
863  }
864 
865 }
866 
867 #endif
Generated by   doxygen 1.8.5