Intrepid
Intrepid_Utils.cpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Intrepid 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 Pavel Bochev (pbboche@sandia.gov)
38// Denis Ridzal (dridzal@sandia.gov), or
39// Kara Peterson (kjpeter@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
49#ifndef INTREPID_UTILS_CPP
50#define INTREPID_UTILS_CPP
51
52#include "Intrepid_Utils.hpp"
53
54namespace Intrepid {
55
56//--------------------------------------------------------------------------------------------//
57// //
58// Definitions: functions for orders, cardinality and etc. of differential operators //
59// //
60//--------------------------------------------------------------------------------------------//
61
62int getFieldRank(const EFunctionSpace spaceType) {
63 int fieldRank = -1;
64
65 switch(spaceType){
66
67 case FUNCTION_SPACE_HGRAD:
68 case FUNCTION_SPACE_HVOL:
69 fieldRank = 0;
70 break;
71
72 case FUNCTION_SPACE_HCURL:
73 case FUNCTION_SPACE_HDIV:
74 case FUNCTION_SPACE_VECTOR_HGRAD:
75 fieldRank = 1;
76 break;
77
78 case FUNCTION_SPACE_TENSOR_HGRAD:
79 fieldRank = 2;
80 break;
81
82 default:
83 TEUCHOS_TEST_FOR_EXCEPTION( !( isValidFunctionSpace(spaceType) ), std::invalid_argument,
84 ">>> ERROR (Intrepid::getFieldRank): Invalid function space type");
85 }
86 return fieldRank;
87}
88
89
90
91int getOperatorRank(const EFunctionSpace spaceType,
92 const EOperator operatorType,
93 const int spaceDim) {
94
95 int fieldRank = Intrepid::getFieldRank(spaceType);
96
97 // Verify arguments: field rank can be 0,1, or 2, spaceDim can be 1,2, or 3.
98#ifdef HAVE_INTREPID_DEBUG
99 TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= fieldRank) && (fieldRank <= 2) ),
100 std::invalid_argument,
101 ">>> ERROR (Intrepid::getOperatorRank): Invalid field rank");
102 TEUCHOS_TEST_FOR_EXCEPTION( !( (1 <= spaceDim ) && (spaceDim <= 3) ),
103 std::invalid_argument,
104 ">>> ERROR (Intrepid::getOperatorRank): Invalid space dimension");
105#endif
106 int operatorRank = -999;
107
108 // In 1D GRAD, CURL, and DIV default to d/dx; Dk defaults to d^k/dx^k, no casing needed.
109 if(spaceDim == 1) {
110 if(fieldRank == 0) {
111
112 // By default, in 1D any operator other than VALUE has rank 1
113 if(operatorType == OPERATOR_VALUE) {
114 operatorRank = 0;
115 }
116 else {
117 operatorRank = 1;
118 }
119 }
120
121 // Only scalar fields are allowed in 1D
122 else {
123 TEUCHOS_TEST_FOR_EXCEPTION( ( fieldRank > 0 ),
124 std::invalid_argument,
125 ">>> ERROR (getOperatorRank): Only scalar fields are allowed in 1D");
126 } // fieldRank == 0
127 } // spaceDim == 1
128
129 // We are either in 2D or 3D
130 else {
131 switch(operatorType) {
132 case OPERATOR_VALUE:
133 operatorRank = 0;
134 break;
135
136 case OPERATOR_GRAD:
137 case OPERATOR_D1:
138 case OPERATOR_D2:
139 case OPERATOR_D3:
140 case OPERATOR_D4:
141 case OPERATOR_D5:
142 case OPERATOR_D6:
143 case OPERATOR_D7:
144 case OPERATOR_D8:
145 case OPERATOR_D9:
146 case OPERATOR_D10:
147 operatorRank = 1;
148 break;
149
150 case OPERATOR_CURL:
151
152 // operator rank for vector and tensor fields equals spaceDim - 3 (-1 in 2D and 0 in 3D)
153 if(fieldRank > 0) {
154 operatorRank = spaceDim - 3;
155 }
156 else {
157
158 // CURL can be applied to scalar functions (rank = 0) in 2D and gives a vector (rank = 1)
159 if(spaceDim == 2) {
160 operatorRank = 3 - spaceDim;
161 }
162
163 // If we are here, fieldRank=0, spaceDim=3: CURL is undefined for 3D scalar functions
164 else {
165 TEUCHOS_TEST_FOR_EXCEPTION( ( (spaceDim == 3) && (fieldRank == 0) ), std::invalid_argument,
166 ">>> ERROR (Intrepid::getOperatorRank): CURL cannot be applied to scalar fields in 3D");
167 }
168 }
169 break;
170
171 case OPERATOR_DIV:
172
173 // DIV can be applied to vectors and tensors and has rank -1 in 2D and 3D
174 if(fieldRank > 0) {
175 operatorRank = -1;
176 }
177
178 // DIV cannot be applied to scalar fields except in 1D where it defaults to d/dx
179 else {
180 TEUCHOS_TEST_FOR_EXCEPTION( ( (spaceDim > 1) && (fieldRank == 0) ), std::invalid_argument,
181 ">>> ERROR (Intrepid::getOperatorRank): DIV cannot be applied to scalar fields in 2D and 3D");
182 }
183 break;
184
185 default:
186 TEUCHOS_TEST_FOR_EXCEPTION( !( isValidOperator(operatorType) ), std::invalid_argument,
187 ">>> ERROR (Intrepid::getOperatorRank): Invalid operator type");
188 } // switch
189 }// 2D and 3D
190
191 return operatorRank;
192}
193
194
195
196int getOperatorOrder(const EOperator operatorType) {
197 int opOrder = -1;
198
199 switch(operatorType){
200
201 case OPERATOR_VALUE:
202 opOrder = 0;
203 break;
204
205 case OPERATOR_GRAD:
206 case OPERATOR_CURL:
207 case OPERATOR_DIV:
208 case OPERATOR_D1:
209 opOrder = 1;
210 break;
211
212 case OPERATOR_D2:
213 case OPERATOR_D3:
214 case OPERATOR_D4:
215 case OPERATOR_D5:
216 case OPERATOR_D6:
217 case OPERATOR_D7:
218 case OPERATOR_D8:
219 case OPERATOR_D9:
220 case OPERATOR_D10:
221 opOrder = (int)operatorType - (int)OPERATOR_D1 + 1;
222 break;
223
224 default:
225 TEUCHOS_TEST_FOR_EXCEPTION( !( Intrepid::isValidOperator(operatorType) ),
226 std::invalid_argument,
227 ">>> ERROR (Intrepid::getOperatorOrder): Invalid operator type");
228 }
229 return opOrder;
230}
231
232
233
234int getDkEnumeration(const int xMult,
235 const int yMult,
236 const int zMult) {
237
238 if( (yMult < 0) && (zMult < 0)) {
239
240#ifdef HAVE_INTREPID_DEBUG
241 // We are in 1D: verify input - xMult is non-negative and total order <= 10:
242 TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= xMult) && (xMult <= INTREPID_MAX_DERIVATIVE) ), std::out_of_range,
243 ">>> ERROR (Intrepid::getDkEnumeration): Derivative order out of range");
244#endif
245
246 // there's only one derivative of order xMult
247 return 0;
248 }
249 else {
250 if( zMult < 0 ) {
251
252#ifdef HAVE_INTREPID_DEBUG
253 // We are in 2D: verify input - xMult and yMult are non-negative and total order <= 10:
254 TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= xMult) && (0 <= yMult) &&
255 ( (xMult + yMult) <= INTREPID_MAX_DERIVATIVE) ), std::out_of_range,
256 ">>> ERROR (Intrepid::getDkEnumeration): Derivative order out of range");
257#endif
258
259 // enumeration is the value of yMult
260 return yMult;
261 }
262
263 // We are in 3D: verify input - xMult, yMult and zMult are non-negative and total order <= 10:
264 else {
265 int order = xMult + yMult + zMult;
266
267#ifdef HAVE_INTREPID_DEBUG
268 // Verify input: total order cannot exceed 10:
269 TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= xMult) && (0 <= yMult) && (0 <= zMult) &&
270 (order <= INTREPID_MAX_DERIVATIVE) ), std::out_of_range,
271 ">>> ERROR (Intrepid::getDkEnumeration): Derivative order out of range");
272#endif
273 int enumeration = zMult;
274 for(int i = 0; i < order - xMult + 1; i++){
275 enumeration += i;
276 }
277 return enumeration;
278 }
279 }
280}
281
282
283
284void getDkMultiplicities(Teuchos::Array<int>& partialMult,
285 const int derivativeEnum,
286 const EOperator operatorType,
287 const int spaceDim) {
288
289 /* Hash table to convert enumeration of partial derivative to multiplicities of dx,dy,dz in 3D.
290 Multiplicities {mx,my,mz} are arranged lexicographically in bins numbered from 0 to 10.
291 The size of bins is an arithmetic progression, i.e., 1,2,3,4,5,...,11. Conversion formula is:
292 \verbatim
293 mx = derivativeOrder - binNumber
294 mz = derivativeEnum - binBegin
295 my = derivativeOrder - mx - mz = binNumber + binBegin - derivativeEnum
296 \endverbatim
297 where binBegin is the enumeration of the first element in the bin. Bin numbers and binBegin
298 values are stored in hash tables for quick access by derivative enumeration value.
299 */
300
301 // Returns the bin number for the specified derivative enumeration
302 static const int binNumber[66] = {
303 0,
304 1, 1,
305 2, 2, 2,
306 3, 3, 3, 3,
307 4, 4, 4, 4, 4,
308 5, 5, 5, 5, 5, 5,
309 6, 6, 6, 6, 6, 6, 6,
310 7, 7, 7, 7, 7, 7, 7, 7,
311 8, 8, 8, 8, 8, 8, 8, 8, 8,
312 9, 9, 9, 9, 9, 9, 9, 9, 9, 9,
313 10,10,10,10,10,10,10,10,10,10,10
314 };
315
316 // Returns the binBegin value for the specified derivative enumeration
317 static const int binBegin[66] ={
318 0,
319 1, 1,
320 3, 3 ,3,
321 6, 6, 6, 6,
322 10,10,10,10,10,
323 15,15,15,15,15,15,
324 21,21,21,21,21,21,21,
325 28,28,28,28,28,28,28,28,
326 36,36,36,36,36,36,36,36,36,
327 45,45,45,45,45,45,45,45,45,45,
328 55,55,55,55,55,55,55,55,55,55,55
329 };
330
331#ifdef HAVE_INTREPID_DEBUG
332 // Enumeration value must be between 0 and the cardinality of the derivative set
333 TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= derivativeEnum) && (derivativeEnum < getDkCardinality(operatorType,spaceDim) ) ),
334 std::invalid_argument,
335 ">>> ERROR (Intrepid::getDkMultiplicities): Invalid derivative enumeration value for this order and space dimension");
336#endif
337
338 // This method should only be called for Dk operators
339 int derivativeOrder;
340 switch(operatorType){
341
342 case OPERATOR_D1:
343 case OPERATOR_D2:
344 case OPERATOR_D3:
345 case OPERATOR_D4:
346 case OPERATOR_D5:
347 case OPERATOR_D6:
348 case OPERATOR_D7:
349 case OPERATOR_D8:
350 case OPERATOR_D9:
351 case OPERATOR_D10:
352 derivativeOrder = Intrepid::getOperatorOrder(operatorType);
353 break;
354
355 default:
356 TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
357 ">>> ERROR (Intrepid::getDkMultiplicities): operator type Dk required for this method");
358 }// switch
359
360 switch(spaceDim) {
361
362 case 1:
363
364 // Resize return array for multiplicity of {dx}
365 partialMult.resize(1);
366
367 // Multiplicity of dx equals derivativeOrder
368 partialMult[0] = derivativeOrder;
369 break;
370
371 case 2:
372
373 // Resize array for multiplicities of {dx,dy}
374 partialMult.resize(2);
375
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
383 // Resize array for multiplicities of {dx,dy,dz}
384 partialMult.resize(3);
385
386 // Recover multiplicities
387 partialMult[0] = derivativeOrder - binNumber[derivativeEnum];
388 partialMult[1] = binNumber[derivativeEnum] + binBegin[derivativeEnum] - derivativeEnum;
389 partialMult[2] = derivativeEnum - binBegin[derivativeEnum];
390 break;
391
392 default:
393 TEUCHOS_TEST_FOR_EXCEPTION( !( (0 < spaceDim ) && (spaceDim < 4) ), std::invalid_argument,
394 ">>> ERROR (Intrepid::getDkMultiplicities): Invalid space dimension");
395 }
396}
397
398
399
400int getDkCardinality(const EOperator operatorType,
401 const int spaceDim) {
402
403 // This should only be called for Dk operators
404 int derivativeOrder;
405 switch(operatorType){
406
407 case OPERATOR_D1:
408 case OPERATOR_D2:
409 case OPERATOR_D3:
410 case OPERATOR_D4:
411 case OPERATOR_D5:
412 case OPERATOR_D6:
413 case OPERATOR_D7:
414 case OPERATOR_D8:
415 case OPERATOR_D9:
416 case OPERATOR_D10:
417 derivativeOrder = Intrepid::getOperatorOrder(operatorType);
418 break;
419
420 default:
421 TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
422 ">>> ERROR (Intrepid::getDkCardinality): operator type Dk required for this method");
423 }// switch
424
425 int cardinality = -999;
426 switch(spaceDim) {
427
428 case 1:
429 cardinality = 1;
430 break;
431
432 case 2:
433 cardinality = derivativeOrder + 1;
434 break;
435
436 case 3:
437 cardinality = (derivativeOrder + 1)*(derivativeOrder + 2)/2;
438 break;
439
440 default:
441 TEUCHOS_TEST_FOR_EXCEPTION( !( (0 < spaceDim ) && (spaceDim < 4) ), std::invalid_argument,
442 ">>> ERROR (Intrepid::getDkcardinality): Invalid space dimension");
443 }
444 return cardinality;
445}
446
447
448
449//--------------------------------------------------------------------------------------------//
450// //
451// Definitions: Helper functions of the Basis class //
452// //
453//--------------------------------------------------------------------------------------------//
454
455void setOrdinalTagData(std::vector<std::vector<std::vector<int> > > &tagToOrdinal,
456 std::vector<std::vector<int> > &ordinalToTag,
457 const int *tags,
458 const int basisCard,
459 const int tagSize,
460 const int posScDim,
461 const int posScOrd,
462 const int posDfOrd) {
463
464
465 // Resize ordinalToTag to a rank-2 array with dimensions (basisCardinality_, 4) and copy tag data
466 ordinalToTag.resize(basisCard);
467 for (int i = 0; i < basisCard; i++) {
468 ordinalToTag[i].resize(4);
469 for (int j = 0; j < tagSize; j++) {
470 ordinalToTag[i][j] = tags[i*tagSize + j];
471 }
472 }
473
474 // Resize tagToOrdinal to a rank-3 array with dimensions (maxScDim + 1, maxScOrd + 1 , maxDfOrd +1)
475 // The 1st dimension of tagToOrdinal is the max value of the 1st column (max subcell dim) in the tag + 1
476 int maxScDim = 0;
477 for (int i = 0; i < basisCard; i++) {
478 if (maxScDim < tags[i*tagSize + posScDim]) {
479 maxScDim = tags[i*tagSize + posScDim];
480 }
481 }
482 maxScDim += 1;
483
484 // The 2nd dimension of tagToOrdinal is the max value of the 2nd column (max subcell id) in the tag + 1
485 int maxScOrd = 0;
486 for (int i = 0; i < basisCard; i++) {
487 if (maxScOrd < tags[i*tagSize + posScOrd]) {
488 maxScOrd = tags[i*tagSize + posScOrd];
489 }
490 }
491 maxScOrd += 1;
492
493 // The 3rd dimension of tagToOrdinal is the max value of the 3rd column (max subcell DofId in the tag) + 1
494 int maxDfOrd = 0;
495 for (int i = 0; i < basisCard; i++) {
496 if (maxDfOrd < tags[i*tagSize + posDfOrd]) {
497 maxDfOrd = tags[i*tagSize + posDfOrd];
498 }
499 }
500 maxDfOrd += 1;
501
502 // Create rank-1 array with dimension maxDfOrd (the 3rd dimension of tagToOrdinal) filled with -1
503 std::vector<int> rank1Array(maxDfOrd, -1);
504
505 // Create rank-2 array with dimensions (maxScOrd, maxDfOrd) (2nd and 3rd dimensions of tagToOrdinal)
506 std::vector<std::vector<int> > rank2Array(maxScOrd, rank1Array);
507
508 // Resize tagToOrdinal to a rank-3 array with dimensions (maxScDim, maxScOrd, maxDfOrd)
509 tagToOrdinal.assign(maxScDim, rank2Array);
510
511 // Overwrite elements of the array corresponding to tags with local DoF Id's, leave all other = -1
512 for (int i = 0; i < basisCard; i++) {
513 tagToOrdinal[tags[i*tagSize]][tags[i*tagSize+1]][tags[i*tagSize+2]] = i;
514 }
515}
516
517
518
519} // end namespace Intrepid
520
521#endif
int isValidOperator(const EOperator operatorType)
Verifies validity of an operator enum.
#define INTREPID_MAX_DERIVATIVE
Maximum order of derivatives allowed in intrepid.
int isValidFunctionSpace(const EFunctionSpace spaceType)
Verifies validity of a function space enum.
void setOrdinalTagData(std::vector< std::vector< std::vector< int > > > &tagToOrdinal, std::vector< std::vector< int > > &ordinalToTag, const int *tags, const int basisCard, const int tagSize, const int posScDim, const int posScOrd, const int posDfOrd)
Fills ordinalToTag_ and tagToOrdinal_ by basis-specific tag data.
int getFieldRank(const EFunctionSpace spaceType)
Returns the rank of fields in a function space of the specified type.
int getOperatorOrder(const EOperator operatorType)
Returns order of an operator.
int getOperatorRank(const EFunctionSpace spaceType, const EOperator operatorType, const int spaceDim)
Returns rank of an operator.
int getDkEnumeration(const int xMult, const int yMult, const int zMult)
Returns the ordinal of a partial derivative of order k based on the multiplicities of the partials dx...
int getDkCardinality(const EOperator operatorType, const int spaceDim)
Returns cardinality of Dk, i.e., the number of all derivatives of order k.
void getDkMultiplicities(Teuchos::Array< int > &partialMult, const int derivativeEnum, const EOperator operatorType, const int spaceDim)
Returns multiplicities of dx, dy, and dz based on the enumeration of the partial derivative,...
Intrepid utilities.