|
Intrepid
|
00001 // @HEADER 00002 // ************************************************************************ 00003 // 00004 // Intrepid Package 00005 // Copyright (2007) Sandia Corporation 00006 // 00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive 00008 // license for use of this work by or on behalf of the U.S. Government. 00009 // 00010 // This library is free software; you can redistribute it and/or modify 00011 // it under the terms of the GNU Lesser General Public License as 00012 // published by the Free Software Foundation; either version 2.1 of the 00013 // License, or (at your option) any later version. 00014 // 00015 // This library is distributed in the hope that it will be useful, but 00016 // WITHOUT ANY WARRANTY; without even the implied warranty of 00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00018 // Lesser General Public License for more details. 00019 // 00020 // You should have received a copy of the GNU Lesser General Public 00021 // License along with this library; if not, write to the Free Software 00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 00023 // USA 00024 // Questions? Contact Pavel Bochev (pbboche@sandia.gov) or 00025 // Denis Ridzal (dridzal@sandia.gov). 00026 // 00027 // ************************************************************************ 00028 // @HEADER 00029 00035 #ifndef INTREPID_UTILS_CPP 00036 #define INTREPID_UTILS_CPP 00037 00038 #include "Intrepid_Utils.hpp" 00039 00040 namespace Intrepid { 00041 00042 //--------------------------------------------------------------------------------------------// 00043 // // 00044 // Definitions: functions for orders, cardinality and etc. of differential operators // 00045 // // 00046 //--------------------------------------------------------------------------------------------// 00047 00048 int getFieldRank(const EFunctionSpace spaceType) { 00049 int fieldRank = -1; 00050 00051 switch(spaceType){ 00052 00053 case FUNCTION_SPACE_HGRAD: 00054 case FUNCTION_SPACE_HVOL: 00055 fieldRank = 0; 00056 break; 00057 00058 case FUNCTION_SPACE_HCURL: 00059 case FUNCTION_SPACE_HDIV: 00060 case FUNCTION_SPACE_VECTOR_HGRAD: 00061 fieldRank = 1; 00062 break; 00063 00064 case FUNCTION_SPACE_TENSOR_HGRAD: 00065 fieldRank = 2; 00066 break; 00067 00068 default: 00069 TEST_FOR_EXCEPTION( !( isValidFunctionSpace(spaceType) ), std::invalid_argument, 00070 ">>> ERROR (Intrepid::getFieldRank): Invalid function space type"); 00071 } 00072 return fieldRank; 00073 } 00074 00075 00076 00077 int getOperatorRank(const EFunctionSpace spaceType, 00078 const EOperator operatorType, 00079 const int spaceDim) { 00080 00081 int fieldRank = Intrepid::getFieldRank(spaceType); 00082 00083 // Verify arguments: field rank can be 0,1, or 2, spaceDim can be 1,2, or 3. 00084 #ifdef HAVE_INTREPID_DEBUG 00085 TEST_FOR_EXCEPTION( !( (0 <= fieldRank) && (fieldRank <= 2) ), 00086 std::invalid_argument, 00087 ">>> ERROR (Intrepid::getOperatorRank): Invalid field rank"); 00088 TEST_FOR_EXCEPTION( !( (1 <= spaceDim ) && (spaceDim <= 3) ), 00089 std::invalid_argument, 00090 ">>> ERROR (Intrepid::getOperatorRank): Invalid space dimension"); 00091 #endif 00092 int operatorRank = -999; 00093 00094 // In 1D GRAD, CURL, and DIV default to d/dx; Dk defaults to d^k/dx^k, no casing needed. 00095 if(spaceDim == 1) { 00096 if(fieldRank == 0) { 00097 00098 // By default, in 1D any operator other than VALUE has rank 1 00099 if(operatorType == OPERATOR_VALUE) { 00100 operatorRank = 0; 00101 } 00102 else { 00103 operatorRank = 1; 00104 } 00105 } 00106 00107 // Only scalar fields are allowed in 1D 00108 else { 00109 TEST_FOR_EXCEPTION( ( fieldRank > 0 ), 00110 std::invalid_argument, 00111 ">>> ERROR (getOperatorRank): Only scalar fields are allowed in 1D"); 00112 } // fieldRank == 0 00113 } // spaceDim == 1 00114 00115 // We are either in 2D or 3D 00116 else { 00117 switch(operatorType) { 00118 case OPERATOR_VALUE: 00119 operatorRank = 0; 00120 break; 00121 00122 case OPERATOR_GRAD: 00123 case OPERATOR_D1: 00124 case OPERATOR_D2: 00125 case OPERATOR_D3: 00126 case OPERATOR_D4: 00127 case OPERATOR_D5: 00128 case OPERATOR_D6: 00129 case OPERATOR_D7: 00130 case OPERATOR_D8: 00131 case OPERATOR_D9: 00132 case OPERATOR_D10: 00133 operatorRank = 1; 00134 break; 00135 00136 case OPERATOR_CURL: 00137 00138 // operator rank for vector and tensor fields equals spaceDim - 3 (-1 in 2D and 0 in 3D) 00139 if(fieldRank > 0) { 00140 operatorRank = spaceDim - 3; 00141 } 00142 else { 00143 00144 // CURL can be applied to scalar functions (rank = 0) in 2D and gives a vector (rank = 1) 00145 if(spaceDim == 2) { 00146 operatorRank = 3 - spaceDim; 00147 } 00148 00149 // If we are here, fieldRank=0, spaceDim=3: CURL is undefined for 3D scalar functions 00150 else { 00151 TEST_FOR_EXCEPTION( ( (spaceDim == 3) && (fieldRank == 0) ), std::invalid_argument, 00152 ">>> ERROR (Intrepid::getOperatorRank): CURL cannot be applied to scalar fields in 3D"); 00153 } 00154 } 00155 break; 00156 00157 case OPERATOR_DIV: 00158 00159 // DIV can be applied to vectors and tensors and has rank -1 in 2D and 3D 00160 if(fieldRank > 0) { 00161 operatorRank = -1; 00162 } 00163 00164 // DIV cannot be applied to scalar fields except in 1D where it defaults to d/dx 00165 else { 00166 TEST_FOR_EXCEPTION( ( (spaceDim > 1) && (fieldRank == 0) ), std::invalid_argument, 00167 ">>> ERROR (Intrepid::getOperatorRank): DIV cannot be applied to scalar fields in 2D and 3D"); 00168 } 00169 break; 00170 00171 default: 00172 TEST_FOR_EXCEPTION( !( isValidOperator(operatorType) ), std::invalid_argument, 00173 ">>> ERROR (Intrepid::getOperatorRank): Invalid operator type"); 00174 } // switch 00175 }// 2D and 3D 00176 00177 return operatorRank; 00178 } 00179 00180 00181 00182 int getOperatorOrder(const EOperator operatorType) { 00183 int opOrder = -1; 00184 00185 switch(operatorType){ 00186 00187 case OPERATOR_VALUE: 00188 opOrder = 0; 00189 break; 00190 00191 case OPERATOR_GRAD: 00192 case OPERATOR_CURL: 00193 case OPERATOR_DIV: 00194 case OPERATOR_D1: 00195 opOrder = 1; 00196 break; 00197 00198 case OPERATOR_D2: 00199 case OPERATOR_D3: 00200 case OPERATOR_D4: 00201 case OPERATOR_D5: 00202 case OPERATOR_D6: 00203 case OPERATOR_D7: 00204 case OPERATOR_D8: 00205 case OPERATOR_D9: 00206 case OPERATOR_D10: 00207 opOrder = (int)operatorType - (int)OPERATOR_D1 + 1; 00208 break; 00209 00210 default: 00211 TEST_FOR_EXCEPTION( !( Intrepid::isValidOperator(operatorType) ), 00212 std::invalid_argument, 00213 ">>> ERROR (Intrepid::getOperatorOrder): Invalid operator type"); 00214 } 00215 return opOrder; 00216 } 00217 00218 00219 00220 int getDkEnumeration(const int xMult, 00221 const int yMult, 00222 const int zMult) { 00223 00224 if( (yMult < 0) && (zMult < 0)) { 00225 00226 #ifdef HAVE_INTREPID_DEBUG 00227 // We are in 1D: verify input - xMult is non-negative and total order <= 10: 00228 TEST_FOR_EXCEPTION( !( (0 <= xMult) && (xMult <= INTREPID_MAX_DERIVATIVE) ), std::out_of_range, 00229 ">>> ERROR (Intrepid::getDkEnumeration): Derivative order out of range"); 00230 #endif 00231 00232 // there's only one derivative of order xMult 00233 return 0; 00234 } 00235 else { 00236 if( zMult < 0 ) { 00237 00238 #ifdef HAVE_INTREPID_DEBUG 00239 // We are in 2D: verify input - xMult and yMult are non-negative and total order <= 10: 00240 TEST_FOR_EXCEPTION( !( (0 <= xMult) && (0 <= yMult) && 00241 ( (xMult + yMult) <= INTREPID_MAX_DERIVATIVE) ), std::out_of_range, 00242 ">>> ERROR (Intrepid::getDkEnumeration): Derivative order out of range"); 00243 #endif 00244 00245 // enumeration is the value of yMult 00246 return yMult; 00247 } 00248 00249 // We are in 3D: verify input - xMult, yMult and zMult are non-negative and total order <= 10: 00250 else { 00251 int order = xMult + yMult + zMult; 00252 00253 #ifdef HAVE_INTREPID_DEBUG 00254 // Verify input: total order cannot exceed 10: 00255 TEST_FOR_EXCEPTION( !( (0 <= xMult) && (0 <= yMult) && (0 <= zMult) && 00256 (order <= INTREPID_MAX_DERIVATIVE) ), std::out_of_range, 00257 ">>> ERROR (Intrepid::getDkEnumeration): Derivative order out of range"); 00258 #endif 00259 int enumeration = zMult; 00260 for(int i = 0; i < order - xMult + 1; i++){ 00261 enumeration += i; 00262 } 00263 return enumeration; 00264 } 00265 } 00266 } 00267 00268 00269 00270 void getDkMultiplicities(Teuchos::Array<int>& partialMult, 00271 const int derivativeEnum, 00272 const EOperator operatorType, 00273 const int spaceDim) { 00274 00275 /* Hash table to convert enumeration of partial derivative to multiplicities of dx,dy,dz in 3D. 00276 Multiplicities {mx,my,mz} are arranged lexicographically in bins numbered from 0 to 10. 00277 The size of bins is an arithmetic progression, i.e., 1,2,3,4,5,...,11. Conversion formula is: 00278 \verbatim 00279 mx = derivativeOrder - binNumber 00280 mz = derivativeEnum - binBegin 00281 my = derivativeOrder - mx - mz = binNumber + binBegin - derivativeEnum 00282 \endverbatim 00283 where binBegin is the enumeration of the first element in the bin. Bin numbers and binBegin 00284 values are stored in hash tables for quick access by derivative enumeration value. 00285 */ 00286 00287 // Returns the bin number for the specified derivative enumeration 00288 static const int binNumber[66] = { 00289 0, 00290 1, 1, 00291 2, 2, 2, 00292 3, 3, 3, 3, 00293 4, 4, 4, 4, 4, 00294 5, 5, 5, 5, 5, 5, 00295 6, 6, 6, 6, 6, 6, 6, 00296 7, 7, 7, 7, 7, 7, 7, 7, 00297 8, 8, 8, 8, 8, 8, 8, 8, 8, 00298 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 00299 10,10,10,10,10,10,10,10,10,10,10 00300 }; 00301 00302 // Returns the binBegin value for the specified derivative enumeration 00303 static const int binBegin[66] ={ 00304 0, 00305 1, 1, 00306 3, 3 ,3, 00307 6, 6, 6, 6, 00308 10,10,10,10,10, 00309 15,15,15,15,15,15, 00310 21,21,21,21,21,21,21, 00311 28,28,28,28,28,28,28,28, 00312 36,36,36,36,36,36,36,36,36, 00313 45,45,45,45,45,45,45,45,45,45, 00314 55,55,55,55,55,55,55,55,55,55,55 00315 }; 00316 00317 #ifdef HAVE_INTREPID_DEBUG 00318 // Enumeration value must be between 0 and the cardinality of the derivative set 00319 TEST_FOR_EXCEPTION( !( (0 <= derivativeEnum) && (derivativeEnum < getDkCardinality(operatorType,spaceDim) ) ), 00320 std::invalid_argument, 00321 ">>> ERROR (Intrepid::getDkMultiplicities): Invalid derivative enumeration value for this order and space dimension"); 00322 #endif 00323 00324 // This method should only be called for Dk operators 00325 int derivativeOrder; 00326 switch(operatorType){ 00327 00328 case OPERATOR_D1: 00329 case OPERATOR_D2: 00330 case OPERATOR_D3: 00331 case OPERATOR_D4: 00332 case OPERATOR_D5: 00333 case OPERATOR_D6: 00334 case OPERATOR_D7: 00335 case OPERATOR_D8: 00336 case OPERATOR_D9: 00337 case OPERATOR_D10: 00338 derivativeOrder = Intrepid::getOperatorOrder(operatorType); 00339 break; 00340 00341 default: 00342 TEST_FOR_EXCEPTION(true, std::invalid_argument, 00343 ">>> ERROR (Intrepid::getDkMultiplicities): operator type Dk required for this method"); 00344 }// switch 00345 00346 switch(spaceDim) { 00347 00348 case 1: 00349 00350 // Resize return array for multiplicity of {dx} 00351 partialMult.resize(1); 00352 00353 // Multiplicity of dx equals derivativeOrder 00354 partialMult[0] = derivativeOrder; 00355 break; 00356 00357 case 2: 00358 00359 // Resize array for multiplicities of {dx,dy} 00360 partialMult.resize(2); 00361 00362 // Multiplicity of dy equals the enumeration of the derivative; of dx - the complement 00363 partialMult[1] = derivativeEnum; 00364 partialMult[0] = derivativeOrder - derivativeEnum; 00365 break; 00366 00367 case 3: 00368 00369 // Resize array for multiplicities of {dx,dy,dz} 00370 partialMult.resize(3); 00371 00372 // Recover multiplicities 00373 partialMult[0] = derivativeOrder - binNumber[derivativeEnum]; 00374 partialMult[1] = binNumber[derivativeEnum] + binBegin[derivativeEnum] - derivativeEnum; 00375 partialMult[2] = derivativeEnum - binBegin[derivativeEnum]; 00376 break; 00377 00378 default: 00379 TEST_FOR_EXCEPTION( !( (0 < spaceDim ) && (spaceDim < 4) ), std::invalid_argument, 00380 ">>> ERROR (Intrepid::getDkMultiplicities): Invalid space dimension"); 00381 } 00382 } 00383 00384 00385 00386 int getDkCardinality(const EOperator operatorType, 00387 const int spaceDim) { 00388 00389 // This should only be called for Dk operators 00390 int derivativeOrder; 00391 switch(operatorType){ 00392 00393 case OPERATOR_D1: 00394 case OPERATOR_D2: 00395 case OPERATOR_D3: 00396 case OPERATOR_D4: 00397 case OPERATOR_D5: 00398 case OPERATOR_D6: 00399 case OPERATOR_D7: 00400 case OPERATOR_D8: 00401 case OPERATOR_D9: 00402 case OPERATOR_D10: 00403 derivativeOrder = Intrepid::getOperatorOrder(operatorType); 00404 break; 00405 00406 default: 00407 TEST_FOR_EXCEPTION(true, std::invalid_argument, 00408 ">>> ERROR (Intrepid::getDkCardinality): operator type Dk required for this method"); 00409 }// switch 00410 00411 int cardinality = -999; 00412 switch(spaceDim) { 00413 00414 case 1: 00415 cardinality = 1; 00416 break; 00417 00418 case 2: 00419 cardinality = derivativeOrder + 1; 00420 break; 00421 00422 case 3: 00423 cardinality = (derivativeOrder + 1)*(derivativeOrder + 2)/2; 00424 break; 00425 00426 default: 00427 TEST_FOR_EXCEPTION( !( (0 < spaceDim ) && (spaceDim < 4) ), std::invalid_argument, 00428 ">>> ERROR (Intrepid::getDkcardinality): Invalid space dimension"); 00429 } 00430 return cardinality; 00431 } 00432 00433 00434 00435 //--------------------------------------------------------------------------------------------// 00436 // // 00437 // Definitions: Helper functions of the Basis class // 00438 // // 00439 //--------------------------------------------------------------------------------------------// 00440 00441 void setOrdinalTagData(std::vector<std::vector<std::vector<int> > > &tagToOrdinal, 00442 std::vector<std::vector<int> > &ordinalToTag, 00443 const int *tags, 00444 const int basisCard, 00445 const int tagSize, 00446 const int posScDim, 00447 const int posScOrd, 00448 const int posDfOrd) { 00449 00450 00451 // Resize ordinalToTag to a rank-2 array with dimensions (basisCardinality_, 4) and copy tag data 00452 ordinalToTag.resize(basisCard); 00453 for (int i = 0; i < basisCard; i++) { 00454 ordinalToTag[i].resize(4); 00455 for (int j = 0; j < tagSize; j++) { 00456 ordinalToTag[i][j] = tags[i*tagSize + j]; 00457 } 00458 } 00459 00460 // Resize tagToOrdinal to a rank-3 array with dimensions (maxScDim + 1, maxScOrd + 1 , maxDfOrd +1) 00461 // The 1st dimension of tagToOrdinal is the max value of the 1st column (max subcell dim) in the tag + 1 00462 int maxScDim = 0; 00463 for (int i = 0; i < basisCard; i++) { 00464 if (maxScDim < tags[i*tagSize + posScDim]) { 00465 maxScDim = tags[i*tagSize + posScDim]; 00466 } 00467 } 00468 maxScDim += 1; 00469 00470 // The 2nd dimension of tagToOrdinal is the max value of the 2nd column (max subcell id) in the tag + 1 00471 int maxScOrd = 0; 00472 for (int i = 0; i < basisCard; i++) { 00473 if (maxScOrd < tags[i*tagSize + posScOrd]) { 00474 maxScOrd = tags[i*tagSize + posScOrd]; 00475 } 00476 } 00477 maxScOrd += 1; 00478 00479 // The 3rd dimension of tagToOrdinal is the max value of the 3rd column (max subcell DofId in the tag) + 1 00480 int maxDfOrd = 0; 00481 for (int i = 0; i < basisCard; i++) { 00482 if (maxDfOrd < tags[i*tagSize + posDfOrd]) { 00483 maxDfOrd = tags[i*tagSize + posDfOrd]; 00484 } 00485 } 00486 maxDfOrd += 1; 00487 00488 // Create rank-1 array with dimension maxDfOrd (the 3rd dimension of tagToOrdinal) filled with -1 00489 std::vector<int> rank1Array(maxDfOrd, -1); 00490 00491 // Create rank-2 array with dimensions (maxScOrd, maxDfOrd) (2nd and 3rd dimensions of tagToOrdinal) 00492 std::vector<std::vector<int> > rank2Array(maxScOrd, rank1Array); 00493 00494 // Resize tagToOrdinal to a rank-3 array with dimensions (maxScDim, maxScOrd, maxDfOrd) 00495 tagToOrdinal.assign(maxScDim, rank2Array); 00496 00497 // Overwrite elements of the array corresponding to tags with local DoF Id's, leave all other = -1 00498 for (int i = 0; i < basisCard; i++) { 00499 tagToOrdinal[tags[i*tagSize]][tags[i*tagSize+1]][tags[i*tagSize+2]] = i; 00500 } 00501 } 00502 00503 00504 00505 } // end namespace Intrepid 00506 00507 #endif
1.7.4