Intrepid2
Intrepid2_UtilsDef.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
48//
49// functions here are moved to basis class
50//
51
52#ifndef __INTREPID2_UTILS_DEF_HPP__
53#define __INTREPID2_UTILS_DEF_HPP__
54
55#include "Intrepid2_Utils.hpp"
56
57namespace Intrepid2 {
58
59 //--------------------------------------------------------------------------------------------//
60 // //
61 // Definitions: functions for orders, cardinality and etc. of differential operators //
62 // //
63 //--------------------------------------------------------------------------------------------//
64
65 KOKKOS_INLINE_FUNCTION
66 ordinal_type getFieldRank(const EFunctionSpace spaceType) {
67 ordinal_type fieldRank = -1;
68
69 switch (spaceType) {
70
71 case FUNCTION_SPACE_HGRAD:
72 case FUNCTION_SPACE_HVOL:
73 fieldRank = 0;
74 break;
75
76 case FUNCTION_SPACE_HCURL:
77 case FUNCTION_SPACE_HDIV:
78 case FUNCTION_SPACE_VECTOR_HGRAD:
79 fieldRank = 1;
80 break;
81
82 case FUNCTION_SPACE_TENSOR_HGRAD:
83 fieldRank = 2;
84 break;
85
86 default:
87 INTREPID2_TEST_FOR_ABORT( !isValidFunctionSpace(spaceType),
88 ">>> ERROR (Intrepid2::getFieldRank): Invalid function space type");
89 }
90 return fieldRank;
91 }
92
93
94 KOKKOS_INLINE_FUNCTION
95 ordinal_type getOperatorRank(const EFunctionSpace spaceType,
96 const EOperator operatorType,
97 const ordinal_type spaceDim) {
98
99 const auto fieldRank = Intrepid2::getFieldRank(spaceType);
100#ifdef HAVE_INTREPID2_DEBUG
101 // Verify arguments: field rank can be 0,1, or 2, spaceDim can be 1,2, or 3.
102 INTREPID2_TEST_FOR_ABORT( !(0 <= fieldRank && fieldRank <= 2),
103 ">>> ERROR (Intrepid2::getOperatorRank): Invalid field rank");
104 INTREPID2_TEST_FOR_ABORT( !(1 <= spaceDim && spaceDim <= 3),
105 ">>> ERROR (Intrepid2::getOperatorRank): Invalid space dimension");
106#endif
107 ordinal_type operatorRank = -999;
108
109 // In 1D GRAD, CURL, and DIV default to d/dx; Dk defaults to d^k/dx^k, no casing needed.
110 if (spaceDim == 1) {
111 if (fieldRank == 0) {
112 // By default, in 1D any operator other than VALUE has rank 1
113 if (operatorType == OPERATOR_VALUE) {
114 operatorRank = 0;
115 } else {
116 operatorRank = 1;
117 }
118 }
119
120 // Only scalar fields are allowed in 1D
121 else {
122 INTREPID2_TEST_FOR_ABORT( fieldRank > 0,
123 ">>> ERROR (getOperatorRank): Only scalar fields are allowed in 1D");
124 } // fieldRank == 0
125 } // spaceDim == 1
126
127 // We are either in 2D or 3D
128 else {
129 switch (operatorType) {
130 case OPERATOR_VALUE:
131 operatorRank = 0;
132 break;
133
134 case OPERATOR_GRAD:
135 case OPERATOR_D1:
136 case OPERATOR_D2:
137 case OPERATOR_D3:
138 case OPERATOR_D4:
139 case OPERATOR_D5:
140 case OPERATOR_D6:
141 case OPERATOR_D7:
142 case OPERATOR_D8:
143 case OPERATOR_D9:
144 case OPERATOR_D10:
145 operatorRank = 1;
146 break;
147
148 case OPERATOR_CURL:
149
150 // operator rank for vector and tensor fields equals spaceDim - 3 (-1 in 2D and 0 in 3D)
151 if (fieldRank > 0) {
152 operatorRank = spaceDim - 3;
153 } else {
154
155 // CURL can be applied to scalar functions (rank = 0) in 2D and gives a vector (rank = 1)
156 if (spaceDim == 2) {
157 operatorRank = 3 - spaceDim;
158 }
159
160 // If we are here, fieldRank=0, spaceDim=3: CURL is undefined for 3D scalar functions
161 else {
162 INTREPID2_TEST_FOR_ABORT( ( (spaceDim == 3) && (fieldRank == 0) ),
163 ">>> ERROR (Intrepid2::getOperatorRank): CURL cannot be applied to scalar fields in 3D");
164 }
165 }
166 break;
167
168 case OPERATOR_DIV:
169
170 // DIV can be applied to vectors and tensors and has rank -1 in 2D and 3D
171 if (fieldRank > 0) {
172 operatorRank = -1;
173 }
174
175 // DIV cannot be applied to scalar fields except in 1D where it defaults to d/dx
176 else {
177 INTREPID2_TEST_FOR_ABORT( ( (spaceDim > 1) && (fieldRank == 0) ),
178 ">>> ERROR (Intrepid2::getOperatorRank): DIV cannot be applied to scalar fields in 2D and 3D");
179 }
180 break;
181
182 default:
183 INTREPID2_TEST_FOR_ABORT( !( isValidOperator(operatorType) ),
184 ">>> ERROR (Intrepid2::getOperatorRank): Invalid operator type");
185 } // switch
186 }// 2D and 3D
187
188 return operatorRank;
189 }
190
191
192 KOKKOS_INLINE_FUNCTION
193 ordinal_type getOperatorOrder(const EOperator operatorType) {
194 ordinal_type opOrder = -1;
195
196 switch (operatorType) {
197
198 case OPERATOR_VALUE:
199 opOrder = 0;
200 break;
201
202 case OPERATOR_GRAD:
203 case OPERATOR_CURL:
204 case OPERATOR_DIV:
205 case OPERATOR_D1:
206 opOrder = 1;
207 break;
208
209 case OPERATOR_D2:
210 case OPERATOR_D3:
211 case OPERATOR_D4:
212 case OPERATOR_D5:
213 case OPERATOR_D6:
214 case OPERATOR_D7:
215 case OPERATOR_D8:
216 case OPERATOR_D9:
217 case OPERATOR_D10:
218 opOrder = (ordinal_type)operatorType - (ordinal_type)OPERATOR_D1 + 1;
219 break;
220
221 default:
222 INTREPID2_TEST_FOR_ABORT( !( Intrepid2::isValidOperator(operatorType) ),
223 ">>> ERROR (Intrepid2::getOperatorOrder): Invalid operator type");
224 }
225 return opOrder;
226 }
227
228
229 KOKKOS_INLINE_FUNCTION
230 ordinal_type getDkEnumeration(const ordinal_type xMult,
231 const ordinal_type yMult,
232 const ordinal_type zMult) {
233
234 if (yMult < 0 && zMult < 0) {
235
236#ifdef HAVE_INTREPID2_DEBUG
237 // We are in 1D: verify input - xMult is non-negative and total order <= 10:
238 INTREPID2_TEST_FOR_ABORT( !( (0 <= xMult) && (xMult <= Parameters::MaxDerivative) ),
239 ">>> ERROR (Intrepid2::getDkEnumeration): Derivative order out of range");
240#endif
241
242 // there's only one derivative of order xMult
243 return 0;
244 } else {
245 if (zMult < 0) {
246
247#ifdef HAVE_INTREPID2_DEBUG
248 // We are in 2D: verify input - xMult and yMult are non-negative and total order <= 10:
249 INTREPID2_TEST_FOR_ABORT( !(0 <= xMult && 0 <= yMult && (xMult + yMult) <= Parameters::MaxDerivative),
250 ">>> ERROR (Intrepid2::getDkEnumeration): Derivative order out of range");
251#endif
252
253 // enumeration is the value of yMult
254 return yMult;
255 }
256
257 // We are in 3D: verify input - xMult, yMult and zMult are non-negative and total order <= 10:
258 else {
259 const auto order = xMult + yMult + zMult;
260#ifdef HAVE_INTREPID2_DEBUG
261 // Verify input: total order cannot exceed 10:
262 INTREPID2_TEST_FOR_ABORT( !( (0 <= xMult) && (0 <= yMult) && (0 <= zMult) &&
263 (order <= Parameters::MaxDerivative) ),
264 ">>> ERROR (Intrepid2::getDkEnumeration): Derivative order out of range");
265#endif
266 ordinal_type enumeration = zMult;
267 const ordinal_type iend = order-xMult+1;
268 for(ordinal_type i=0;i<iend;++i) {
269 enumeration += i;
270 }
271 return enumeration;
272 }
273 }
274 }
275
276// template<typename OrdinalArrayType>
277// KOKKOS_INLINE_FUNCTION
278// void getDkMultiplicities(OrdinalArrayType partialMult,
279// const ordinal_type derivativeEnum,
280// const EOperator operatorType,
281// const ordinal_type spaceDim) {
282
283// /* Hash table to convert enumeration of partial derivative to multiplicities of dx,dy,dz in 3D.
284// Multiplicities {mx,my,mz} are arranged lexicographically in bins numbered from 0 to 10.
285// The size of bins is an arithmetic progression, i.e., 1,2,3,4,5,...,11. Conversion formula is:
286// \verbatim
287// mx = derivativeOrder - binNumber
288// mz = derivativeEnum - binBegin
289// my = derivativeOrder - mx - mz = binNumber + binBegin - derivativeEnum
290// \endverbatim
291// where binBegin is the enumeration of the first element in the bin. Bin numbers and binBegin
292// values are stored in hash tables for quick access by derivative enumeration value.
293// */
294
295// // Returns the bin number for the specified derivative enumeration
296// const ordinal_type binNumber[66] = {
297// 0,
298// 1, 1,
299// 2, 2, 2,
300// 3, 3, 3, 3,
301// 4, 4, 4, 4, 4,
302// 5, 5, 5, 5, 5, 5,
303// 6, 6, 6, 6, 6, 6, 6,
304// 7, 7, 7, 7, 7, 7, 7, 7,
305// 8, 8, 8, 8, 8, 8, 8, 8, 8,
306// 9, 9, 9, 9, 9, 9, 9, 9, 9, 9,
307// 10,10,10,10,10,10,10,10,10,10,10
308// };
309
310// // Returns the binBegin value for the specified derivative enumeration
311// const ordinal_type binBegin[66] = {
312// 0,
313// 1, 1,
314// 3, 3 ,3,
315// 6, 6, 6, 6,
316// 10,10,10,10,10,
317// 15,15,15,15,15,15,
318// 21,21,21,21,21,21,21,
319// 28,28,28,28,28,28,28,28,
320// 36,36,36,36,36,36,36,36,36,
321// 45,45,45,45,45,45,45,45,45,45,
322// 55,55,55,55,55,55,55,55,55,55,55
323// };
324
325// #ifdef HAVE_INTREPID2_DEBUG
326// // Enumeration value must be between 0 and the cardinality of the derivative set
327// INTREPID2_TEST_FOR_ABORT( !( (0 <= derivativeEnum) && (derivativeEnum < getDkCardinality(operatorType,spaceDim) ) ),
328// ">>> ERROR (Intrepid2::getDkMultiplicities): Invalid derivative enumeration value for this order and space dimension");
329// #endif
330
331// // This method should only be called for Dk operators
332// ordinal_type derivativeOrder;
333// switch(operatorType) {
334
335// case OPERATOR_D1:
336// case OPERATOR_D2:
337// case OPERATOR_D3:
338// case OPERATOR_D4:
339// case OPERATOR_D5:
340// case OPERATOR_D6:
341// case OPERATOR_D7:
342// case OPERATOR_D8:
343// case OPERATOR_D9:
344// case OPERATOR_D10:
345// derivativeOrder = Intrepid2::getOperatorOrder(operatorType);
346// break;
347
348// default:
349// INTREPID2_TEST_FOR_ABORT(true,
350// ">>> ERROR (Intrepid2::getDkMultiplicities): operator type Dk required for this method");
351// }// switch
352
353// #ifdef HAVE_INTREPID2_DEBUG
354// INTREPID2_TEST_FOR_ABORT( partialMult.extent(0) != spaceDim,
355// ">>> ERROR (Intrepid2::getDkMultiplicities): partialMult must have the same dimension as spaceDim" );
356// #endif
357
358// switch(spaceDim) {
359
360// case 1:
361// // Multiplicity of dx equals derivativeOrder
362// partialMult(0) = derivativeOrder;
363// break;
364
365// case 2:
366// // Multiplicity of dy equals the enumeration of the derivative; of dx - the complement
367// partialMult(1) = derivativeEnum;
368// partialMult(0) = derivativeOrder - derivativeEnum;
369// break;
370
371// case 3:
372// // Recover multiplicities
373// partialMult(0) = derivativeOrder - binNumber[derivativeEnum];
374// partialMult(1) = binNumber[derivativeEnum] + binBegin[derivativeEnum] - derivativeEnum;
375// partialMult(2) = derivativeEnum - binBegin[derivativeEnum];
376// break;
377
378// default:
379// INTREPID2_TEST_FOR_ABORT( !( (0 < spaceDim ) && (spaceDim < 4) ),
380// ">>> ERROR (Intrepid2::getDkMultiplicities): Invalid space dimension");
381// }
382// }
383
384
385 KOKKOS_INLINE_FUNCTION
386 ordinal_type getDkCardinality(const EOperator operatorType,
387 const ordinal_type spaceDim) {
388
389 // This should only be called for Dk operators
390 ordinal_type derivativeOrder;
391 switch(operatorType) {
392
393 case OPERATOR_D1:
394 case OPERATOR_D2:
395 case OPERATOR_D3:
396 case OPERATOR_D4:
397 case OPERATOR_D5:
398 case OPERATOR_D6:
399 case OPERATOR_D7:
400 case OPERATOR_D8:
401 case OPERATOR_D9:
402 case OPERATOR_D10:
403 derivativeOrder = Intrepid2::getOperatorOrder(operatorType);
404 break;
405
406 default:
407 INTREPID2_TEST_FOR_ABORT(true,
408 ">>> ERROR (Intrepid2::getDkCardinality): operator type Dk required for this method");
409 }// switch
410
411 ordinal_type cardinality = -999;
412 switch(spaceDim) {
413
414 case 1:
415 cardinality = 1;
416 break;
417
418 case 2:
419 cardinality = derivativeOrder + 1;
420 break;
421
422 case 3:
423 cardinality = (derivativeOrder + 1)*(derivativeOrder + 2)/2;
424 break;
425
426 default:
427 INTREPID2_TEST_FOR_ABORT( !( (0 < spaceDim ) && (spaceDim < 4) ),
428 ">>> ERROR (Intrepid2::getDkcardinality): Invalid space dimension");
429 }
430
431 return cardinality;
432 }
433
434} // end namespace Intrepid2
435
436#endif
KOKKOS_INLINE_FUNCTION ordinal_type getOperatorRank(const EFunctionSpace spaceType, const EOperator operatorType, const ordinal_type spaceDim)
Returns rank of an operator.
KOKKOS_INLINE_FUNCTION ordinal_type getFieldRank(const EFunctionSpace spaceType)
Returns the rank of fields in a function space of the specified type.
KOKKOS_INLINE_FUNCTION ordinal_type getOperatorOrder(const EOperator operatorType)
Returns order of an operator.
KOKKOS_INLINE_FUNCTION ordinal_type getDkCardinality(const EOperator operatorType, const ordinal_type spaceDim)
Returns cardinality of Dk, i.e., the number of all derivatives of order k.
KOKKOS_INLINE_FUNCTION ordinal_type getDkEnumeration(const ordinal_type xMult, const ordinal_type yMult=-1, const ordinal_type zMult=-1)
Returns the ordinal of a partial derivative of order k based on the multiplicities of the partials dx...
KOKKOS_FORCEINLINE_FUNCTION bool isValidFunctionSpace(const EFunctionSpace spaceType)
Verifies validity of a function space enum.
KOKKOS_FORCEINLINE_FUNCTION bool isValidOperator(const EOperator operatorType)
Verifies validity of an operator enum.
Header function for Intrepid2::Util class and other utility functions.
static constexpr ordinal_type MaxDerivative
Maximum order of derivatives allowed in intrepid.