Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/src/Shared/Intrepid_Utils.cpp
Go to the documentation of this file.
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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Pavel Bochev  (pbboche@sandia.gov)
00038 //                    Denis Ridzal  (dridzal@sandia.gov), or
00039 //                    Kara Peterson (kjpeter@sandia.gov)
00040 //
00041 // ************************************************************************
00042 // @HEADER
00043 
00049 #ifndef INTREPID_UTILS_CPP
00050 #define INTREPID_UTILS_CPP
00051 
00052 #include "Intrepid_Utils.hpp"
00053 
00054 namespace Intrepid {
00055   
00056 //--------------------------------------------------------------------------------------------//
00057 //                                                                                            //
00058 //    Definitions: functions for orders, cardinality and etc. of differential operators       //
00059 //                                                                                            //
00060 //--------------------------------------------------------------------------------------------//
00061 
00062 int getFieldRank(const EFunctionSpace spaceType) {
00063   int fieldRank = -1;
00064   
00065   switch(spaceType){
00066     
00067     case FUNCTION_SPACE_HGRAD:
00068     case FUNCTION_SPACE_HVOL:
00069       fieldRank = 0;
00070       break;
00071       
00072     case FUNCTION_SPACE_HCURL:
00073     case FUNCTION_SPACE_HDIV:
00074     case FUNCTION_SPACE_VECTOR_HGRAD:
00075       fieldRank = 1;
00076       break;
00077       
00078     case FUNCTION_SPACE_TENSOR_HGRAD:
00079       fieldRank = 2;
00080       break;
00081       
00082     default:
00083       TEUCHOS_TEST_FOR_EXCEPTION( !( isValidFunctionSpace(spaceType) ), std::invalid_argument,
00084                           ">>> ERROR (Intrepid::getFieldRank): Invalid function space type");
00085   }
00086   return fieldRank;
00087 }
00088 
00089 
00090 
00091 int getOperatorRank(const EFunctionSpace spaceType,
00092                     const EOperator      operatorType,
00093                     const int            spaceDim) {
00094   
00095   int fieldRank = Intrepid::getFieldRank(spaceType);
00096   
00097   // Verify arguments: field rank can be 0,1, or 2, spaceDim can be 1,2, or 3.
00098 #ifdef HAVE_INTREPID_DEBUG
00099   TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= fieldRank) && (fieldRank <= 2) ),
00100                       std::invalid_argument,
00101                       ">>> ERROR (Intrepid::getOperatorRank): Invalid field rank");
00102   TEUCHOS_TEST_FOR_EXCEPTION( !( (1 <= spaceDim ) && (spaceDim  <= 3) ),
00103                       std::invalid_argument,
00104                       ">>> ERROR (Intrepid::getOperatorRank): Invalid space dimension");
00105 #endif
00106   int operatorRank = -999;
00107   
00108   // In 1D GRAD, CURL, and DIV default to d/dx; Dk defaults to d^k/dx^k, no casing needed.
00109   if(spaceDim == 1) {
00110     if(fieldRank == 0) {
00111       
00112       // By default, in 1D any operator other than VALUE has rank 1
00113       if(operatorType == OPERATOR_VALUE) {
00114         operatorRank = 0;
00115       }
00116       else {
00117         operatorRank = 1;
00118       }
00119     }
00120     
00121     // Only scalar fields are allowed in 1D
00122     else {
00123       TEUCHOS_TEST_FOR_EXCEPTION( ( fieldRank > 0 ),
00124                           std::invalid_argument,
00125                           ">>> ERROR (getOperatorRank): Only scalar fields are allowed in 1D");  
00126     } // fieldRank == 0
00127   } // spaceDim == 1
00128   
00129   // We are either in 2D or 3D
00130   else {  
00131     switch(operatorType) {
00132       case OPERATOR_VALUE:
00133         operatorRank = 0;
00134         break;
00135         
00136       case OPERATOR_GRAD:
00137       case OPERATOR_D1:
00138       case OPERATOR_D2:
00139       case OPERATOR_D3:
00140       case OPERATOR_D4:
00141       case OPERATOR_D5:
00142       case OPERATOR_D6:
00143       case OPERATOR_D7:
00144       case OPERATOR_D8:
00145       case OPERATOR_D9:
00146       case OPERATOR_D10:
00147         operatorRank = 1;
00148         break;
00149         
00150       case OPERATOR_CURL:
00151         
00152         // operator rank for vector and tensor fields equals spaceDim - 3 (-1 in 2D and 0 in 3D)   
00153         if(fieldRank > 0) {
00154           operatorRank = spaceDim - 3;
00155         }
00156         else {
00157           
00158           // CURL can be applied to scalar functions (rank = 0) in 2D and gives a vector (rank = 1)
00159           if(spaceDim == 2) {
00160             operatorRank = 3 - spaceDim;
00161           }
00162           
00163           // If we are here, fieldRank=0, spaceDim=3: CURL is undefined for 3D scalar functions
00164           else {
00165             TEUCHOS_TEST_FOR_EXCEPTION( ( (spaceDim == 3) && (fieldRank == 0) ), std::invalid_argument,
00166                                 ">>> ERROR (Intrepid::getOperatorRank): CURL cannot be applied to scalar fields in 3D");  
00167           }
00168         }
00169         break;
00170         
00171       case OPERATOR_DIV:
00172         
00173         // DIV can be applied to vectors and tensors and has rank -1 in 2D and 3D
00174         if(fieldRank > 0) {
00175           operatorRank = -1; 
00176         }
00177         
00178         // DIV cannot be applied to scalar fields except in 1D where it defaults to d/dx
00179         else {
00180           TEUCHOS_TEST_FOR_EXCEPTION( ( (spaceDim > 1) && (fieldRank == 0) ), std::invalid_argument,
00181                               ">>> ERROR (Intrepid::getOperatorRank): DIV cannot be applied to scalar fields in 2D and 3D");  
00182         }
00183         break;
00184         
00185       default:
00186         TEUCHOS_TEST_FOR_EXCEPTION( !( isValidOperator(operatorType) ), std::invalid_argument,
00187                             ">>> ERROR (Intrepid::getOperatorRank): Invalid operator type");
00188     } // switch
00189   }// 2D and 3D
00190   
00191   return operatorRank; 
00192 }
00193 
00194 
00195 
00196 int getOperatorOrder(const EOperator operatorType) {
00197   int opOrder = -1;
00198   
00199   switch(operatorType){
00200     
00201     case OPERATOR_VALUE:
00202       opOrder = 0;
00203       break;
00204       
00205     case OPERATOR_GRAD:
00206     case OPERATOR_CURL:
00207     case OPERATOR_DIV:
00208     case OPERATOR_D1:
00209       opOrder = 1;
00210       break;
00211       
00212     case OPERATOR_D2:
00213     case OPERATOR_D3:
00214     case OPERATOR_D4:
00215     case OPERATOR_D5:
00216     case OPERATOR_D6:
00217     case OPERATOR_D7:
00218     case OPERATOR_D8:
00219     case OPERATOR_D9:
00220     case OPERATOR_D10:
00221       opOrder = (int)operatorType - (int)OPERATOR_D1 + 1;
00222       break;
00223       
00224     default:
00225       TEUCHOS_TEST_FOR_EXCEPTION( !( Intrepid::isValidOperator(operatorType) ),
00226                           std::invalid_argument,
00227                           ">>> ERROR (Intrepid::getOperatorOrder): Invalid operator type");
00228   }
00229   return opOrder;
00230 }    
00231 
00232 
00233 
00234 int getDkEnumeration(const int xMult,
00235                      const int yMult,
00236                      const int zMult) {
00237   
00238   if( (yMult < 0) && (zMult < 0)) {
00239     
00240 #ifdef HAVE_INTREPID_DEBUG
00241     // We are in 1D: verify input - xMult is non-negative  and total order <= 10:
00242     TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= xMult) && (xMult <= INTREPID_MAX_DERIVATIVE) ), std::out_of_range,
00243                         ">>> ERROR (Intrepid::getDkEnumeration): Derivative order out of range");
00244 #endif
00245     
00246     // there's only one derivative of order xMult
00247     return 0;
00248   }
00249   else {
00250     if( zMult < 0 ) {
00251       
00252 #ifdef HAVE_INTREPID_DEBUG
00253       // We are in 2D: verify input - xMult and yMult are non-negative and total order <= 10:
00254       TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= xMult) && (0 <= yMult) && 
00255                              ( (xMult + yMult) <= INTREPID_MAX_DERIVATIVE) ), std::out_of_range,
00256                           ">>> ERROR (Intrepid::getDkEnumeration): Derivative order out of range");
00257 #endif
00258       
00259       // enumeration is the value of yMult
00260       return yMult;
00261     }
00262     
00263     // We are in 3D: verify input - xMult, yMult and zMult are non-negative and total order <= 10:
00264     else {
00265       int order = xMult + yMult + zMult;
00266       
00267 #ifdef HAVE_INTREPID_DEBUG
00268       // Verify input:  total order cannot exceed 10:
00269       TEUCHOS_TEST_FOR_EXCEPTION(  !( (0 <= xMult) && (0 <= yMult) && (0 <= zMult) && 
00270                               (order <= INTREPID_MAX_DERIVATIVE) ), std::out_of_range,
00271                            ">>> ERROR (Intrepid::getDkEnumeration): Derivative order out of range");
00272 #endif
00273       int enumeration = zMult;
00274       for(int i = 0; i < order - xMult + 1; i++){
00275         enumeration += i; 
00276       }
00277       return enumeration;
00278     }
00279   }
00280 }
00281 
00282 
00283 
00284 void getDkMultiplicities(Teuchos::Array<int>&  partialMult,
00285                          const int             derivativeEnum,
00286                          const EOperator       operatorType,
00287                          const int             spaceDim) {
00288   
00289   /* Hash table to convert enumeration of partial derivative to multiplicities of dx,dy,dz in 3D.
00290   Multiplicities {mx,my,mz} are arranged lexicographically in bins numbered from 0 to 10. 
00291   The size of bins is an arithmetic progression, i.e., 1,2,3,4,5,...,11. Conversion formula is:
00292   \verbatim
00293     mx = derivativeOrder - binNumber
00294     mz = derivativeEnum  - binBegin
00295     my = derivativeOrder - mx - mz = binNumber + binBegin - derivativeEnum
00296   \endverbatim
00297   where binBegin is the enumeration of the first element in the bin. Bin numbers and binBegin 
00298   values are stored in hash tables for quick access by derivative enumeration value. 
00299   */ 
00300   
00301   // Returns the bin number for the specified derivative enumeration
00302   static const int binNumber[66] = { 
00303     0,
00304     1, 1,
00305     2, 2, 2,
00306     3, 3, 3, 3,
00307     4, 4, 4, 4, 4,
00308     5, 5, 5, 5, 5, 5,
00309     6, 6, 6, 6, 6, 6, 6,
00310     7, 7, 7, 7, 7, 7, 7, 7,
00311     8, 8, 8, 8, 8, 8, 8, 8, 8,
00312     9, 9, 9, 9, 9, 9, 9, 9, 9, 9,
00313     10,10,10,10,10,10,10,10,10,10,10
00314   };
00315   
00316   // Returns the binBegin value for the specified derivative enumeration 
00317   static const int binBegin[66] ={ 
00318     0,
00319     1, 1,
00320     3, 3 ,3,
00321     6, 6, 6, 6,
00322     10,10,10,10,10,
00323     15,15,15,15,15,15,
00324     21,21,21,21,21,21,21,
00325     28,28,28,28,28,28,28,28,
00326     36,36,36,36,36,36,36,36,36,
00327     45,45,45,45,45,45,45,45,45,45,
00328     55,55,55,55,55,55,55,55,55,55,55
00329   };
00330     
00331 #ifdef HAVE_INTREPID_DEBUG
00332   // Enumeration value must be between 0 and the cardinality of the derivative set
00333   TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= derivativeEnum) && (derivativeEnum < getDkCardinality(operatorType,spaceDim) ) ),
00334                       std::invalid_argument,
00335                       ">>> ERROR (Intrepid::getDkMultiplicities): Invalid derivative enumeration value for this order and space dimension");
00336 #endif
00337   
00338   // This method should only be called for Dk operators
00339   int derivativeOrder;
00340   switch(operatorType){
00341     
00342     case OPERATOR_D1:
00343     case OPERATOR_D2:
00344     case OPERATOR_D3:
00345     case OPERATOR_D4:
00346     case OPERATOR_D5:
00347     case OPERATOR_D6:
00348     case OPERATOR_D7:
00349     case OPERATOR_D8:
00350     case OPERATOR_D9:
00351     case OPERATOR_D10:
00352       derivativeOrder = Intrepid::getOperatorOrder(operatorType);
00353       break;
00354       
00355     default:
00356       TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
00357                          ">>> ERROR (Intrepid::getDkMultiplicities): operator type Dk required for this method");
00358   }// switch
00359   
00360   switch(spaceDim) {
00361     
00362     case 1:
00363       
00364       // Resize return array for multiplicity of {dx}
00365       partialMult.resize(1);
00366       
00367       // Multiplicity of dx equals derivativeOrder
00368       partialMult[0] = derivativeOrder;
00369       break;
00370       
00371     case 2:
00372       
00373       // Resize array for multiplicities of {dx,dy}
00374       partialMult.resize(2);
00375       
00376       // Multiplicity of dy equals the enumeration of the derivative; of dx - the complement
00377       partialMult[1] = derivativeEnum;
00378       partialMult[0] = derivativeOrder - derivativeEnum;
00379       break;
00380       
00381     case 3:
00382       
00383       // Resize array for multiplicities of {dx,dy,dz}
00384       partialMult.resize(3);
00385       
00386       // Recover multiplicities
00387       partialMult[0] = derivativeOrder - binNumber[derivativeEnum];
00388       partialMult[1] = binNumber[derivativeEnum] + binBegin[derivativeEnum] - derivativeEnum;
00389       partialMult[2] = derivativeEnum  -  binBegin[derivativeEnum];
00390       break;
00391       
00392     default:
00393       TEUCHOS_TEST_FOR_EXCEPTION( !( (0 < spaceDim ) && (spaceDim < 4) ), std::invalid_argument,
00394                           ">>> ERROR (Intrepid::getDkMultiplicities): Invalid space dimension");          
00395   }
00396 }
00397 
00398 
00399 
00400 int getDkCardinality(const EOperator operatorType,
00401                      const int       spaceDim) {
00402   
00403   // This should only be called for Dk operators
00404   int derivativeOrder;
00405   switch(operatorType){
00406     
00407     case OPERATOR_D1:
00408     case OPERATOR_D2:
00409     case OPERATOR_D3:
00410     case OPERATOR_D4:
00411     case OPERATOR_D5:
00412     case OPERATOR_D6:
00413     case OPERATOR_D7:
00414     case OPERATOR_D8:
00415     case OPERATOR_D9:
00416     case OPERATOR_D10:
00417       derivativeOrder = Intrepid::getOperatorOrder(operatorType);
00418       break;
00419       
00420     default:
00421       TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
00422                         ">>> ERROR (Intrepid::getDkCardinality): operator type Dk required for this method");
00423   }// switch
00424 
00425   int cardinality = -999;
00426   switch(spaceDim) {
00427     
00428     case 1:
00429       cardinality = 1;
00430       break;
00431       
00432     case 2:
00433       cardinality = derivativeOrder + 1;
00434       break;
00435       
00436     case 3:
00437       cardinality = (derivativeOrder + 1)*(derivativeOrder + 2)/2;
00438       break;
00439       
00440     default:
00441       TEUCHOS_TEST_FOR_EXCEPTION( !( (0 < spaceDim ) && (spaceDim < 4) ), std::invalid_argument,
00442                           ">>> ERROR (Intrepid::getDkcardinality): Invalid space dimension");          
00443   }
00444   return cardinality;
00445 }
00446 
00447 
00448 
00449 //--------------------------------------------------------------------------------------------//
00450 //                                                                                            //
00451 //            Definitions:    Helper functions of the Basis class                             //
00452 //                                                                                            //
00453 //--------------------------------------------------------------------------------------------//
00454 
00455 void setOrdinalTagData(std::vector<std::vector<std::vector<int> > >    &tagToOrdinal,
00456                        std::vector<std::vector<int> >                  &ordinalToTag,
00457                        const int                                       *tags,
00458                        const int                                       basisCard,
00459                        const int                                       tagSize,
00460                        const int                                       posScDim,
00461                        const int                                       posScOrd,
00462                        const int                                       posDfOrd) {
00463 
00464 
00465   // Resize ordinalToTag to a rank-2 array with dimensions (basisCardinality_, 4) and copy tag data
00466   ordinalToTag.resize(basisCard);
00467   for (int i = 0; i < basisCard; i++) {
00468     ordinalToTag[i].resize(4);
00469     for (int j = 0; j < tagSize; j++) {
00470       ordinalToTag[i][j] = tags[i*tagSize + j];
00471     }
00472   }
00473 
00474   // Resize tagToOrdinal to a rank-3 array with dimensions (maxScDim + 1, maxScOrd + 1 , maxDfOrd +1)
00475   // The 1st dimension of tagToOrdinal is the max value of the 1st column (max subcell dim) in the tag + 1
00476   int maxScDim = 0;
00477   for (int i = 0; i < basisCard; i++) {
00478     if (maxScDim < tags[i*tagSize + posScDim]) {
00479       maxScDim = tags[i*tagSize + posScDim];
00480     }
00481   }
00482   maxScDim += 1;
00483 
00484   // The 2nd dimension of tagToOrdinal is the max value of the 2nd column (max subcell id) in the tag  + 1
00485   int maxScOrd = 0;
00486   for (int i = 0; i < basisCard; i++) {
00487     if (maxScOrd < tags[i*tagSize + posScOrd]) {
00488       maxScOrd = tags[i*tagSize + posScOrd];
00489     }
00490   }
00491   maxScOrd += 1;
00492 
00493   // The 3rd dimension of tagToOrdinal is the max value of the 3rd column (max subcell DofId in the tag) + 1
00494   int maxDfOrd = 0;
00495   for (int i = 0; i < basisCard; i++) {
00496     if (maxDfOrd < tags[i*tagSize + posDfOrd]) {
00497       maxDfOrd = tags[i*tagSize + posDfOrd];
00498     }
00499   }
00500   maxDfOrd += 1;
00501 
00502   // Create rank-1 array with dimension maxDfOrd (the 3rd dimension of tagToOrdinal) filled with -1
00503   std::vector<int> rank1Array(maxDfOrd, -1);
00504 
00505   // Create rank-2 array with dimensions (maxScOrd, maxDfOrd) (2nd and 3rd dimensions of tagToOrdinal)
00506   std::vector<std::vector<int> > rank2Array(maxScOrd, rank1Array);
00507 
00508   // Resize tagToOrdinal to a rank-3 array with dimensions (maxScDim, maxScOrd, maxDfOrd)
00509   tagToOrdinal.assign(maxScDim, rank2Array);
00510 
00511   // Overwrite elements of the array corresponding to tags with local DoF Id's, leave all other = -1
00512   for (int i = 0; i < basisCard; i++) {
00513     tagToOrdinal[tags[i*tagSize]][tags[i*tagSize+1]][tags[i*tagSize+2]] = i;
00514   }
00515 }
00516 
00517 
00518 
00519 } // end namespace Intrepid
00520 
00521 #endif