Intrepid: http://trilinos.sandia.gov/packages/docs/r10.0/packages/intrepid/src/Discretization/Basis/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 
00036 template<class Scalar, class ArrayScalar>
00037 int Basis<Scalar, ArrayScalar>::getDofOrdinal(const int subcDim,
00038                                               const int subcOrd,
00039                                               const int subcDofOrd) {
00040   if (!basisTagsAreSet_) {
00041     initializeTags();
00042     basisTagsAreSet_ = true;
00043   }
00044   // Use .at() for bounds checking
00045   int dofOrdinal = tagToOrdinal_.at(subcDim).at(subcOrd).at(subcDofOrd);
00046   TEST_FOR_EXCEPTION( (dofOrdinal == -1), std::invalid_argument, 
00047                       ">>> ERROR (Basis): Invalid DoF tag");
00048   return dofOrdinal;
00049 }
00050 
00051 
00052 template<class Scalar,class ArrayScalar>
00053 const std::vector<std::vector<std::vector<int> > > & Basis<Scalar, ArrayScalar>::getDofOrdinalData( ) 
00054 {
00055   if (!basisTagsAreSet_) {
00056     initializeTags();
00057     basisTagsAreSet_ = true;
00058   }
00059   return tagToOrdinal_;
00060 }
00061 
00062 
00063 template<class Scalar, class ArrayScalar>
00064 const std::vector<int>&  Basis<Scalar, ArrayScalar>::getDofTag(int dofOrd) {
00065   if (!basisTagsAreSet_) {
00066     initializeTags();
00067     basisTagsAreSet_ = true;
00068   }
00069   // Use .at() for bounds checking
00070   return ordinalToTag_.at(dofOrd);
00071 }
00072 
00073 
00074 template<class Scalar, class ArrayScalar>
00075 const std::vector<std::vector<int> > & Basis<Scalar, ArrayScalar>::getAllDofTags() {
00076   if (!basisTagsAreSet_) {
00077     initializeTags();
00078     basisTagsAreSet_ = true;
00079   }
00080   return ordinalToTag_;
00081 }
00082 
00083 
00084 
00085 template<class Scalar, class ArrayScalar>
00086 inline int Basis<Scalar, ArrayScalar>::getCardinality() const {
00087   return basisCardinality_;   
00088 }
00089 
00090 
00091 template<class Scalar, class ArrayScalar>
00092 inline EBasis Basis<Scalar, ArrayScalar>::getBasisType() const {
00093   return basisType_;
00094 }
00095 
00096 
00097 template<class Scalar, class ArrayScalar>
00098 inline const shards::CellTopology Basis<Scalar, ArrayScalar>::getBaseCellTopology() const {
00099   return basisCellTopology_;
00100 }
00101 
00102 
00103 template<class Scalar, class ArrayScalar>
00104 inline int Basis<Scalar,ArrayScalar>::getDegree() const {
00105   return basisDegree_;
00106 }
00107 
00108 
00109 template<class Scalar, class ArrayScalar>
00110 inline ECoordinates Basis<Scalar, ArrayScalar>::getCoordinateSystem() const {
00111   return basisCoordinates_;
00112 }
00113   
00114 
00115 //--------------------------------------------------------------------------------------------//
00116 //                                                                                            //            
00117 //                            Helper functions of the Basis class                             //
00118 //                                                                                            //            
00119 //--------------------------------------------------------------------------------------------//
00120 
00121 template<class Scalar, class ArrayScalar>
00122 void getValues_HGRAD_Args(ArrayScalar &                outputValues,
00123                           const ArrayScalar &          inputPoints,
00124                           const EOperator              operatorType,
00125                           const shards::CellTopology&  cellTopo,
00126                           const int                    basisCard){
00127   
00128   int spaceDim = cellTopo.getDimension();
00129   
00130   // Verify inputPoints array
00131   TEST_FOR_EXCEPTION( !(inputPoints.rank() == 2), std::invalid_argument, 
00132                       ">>> ERROR: (Intrepid::getValues_HGRAD_Args) rank = 2 required for inputPoints array");
00133   
00134   TEST_FOR_EXCEPTION(  (inputPoints.dimension(0) <= 0), std::invalid_argument,
00135                       ">>> ERROR (Intrepid::getValues_HGRAD_Args): dim 0 (number of points) > 0 required for inputPoints array");
00136 
00137   TEST_FOR_EXCEPTION( !(inputPoints.dimension(1) == spaceDim), std::invalid_argument,
00138                       ">>> ERROR: (Intrepid::getValues_HGRAD_Args) dim 1 (spatial dimension) of inputPoints array  does not match cell dimension");
00139 
00140   
00141   // Verify that all inputPoints are in the reference cell
00142   /*
00143    TEST_FOR_EXCEPTION( !CellTools<Scalar>::checkPointSetInclusion(inputPoints, cellTopo), std::invalid_argument,
00144                        ">>> ERROR: (Intrepid::getValues_HGRAD_Args) One or more points are outside the " 
00145                        << cellTopo <<" reference cell");
00146    */
00147   
00148   
00149   // Verify that operatorType is admissible for HGRAD fields
00150   TEST_FOR_EXCEPTION( ( (spaceDim == 2) && (operatorType == OPERATOR_DIV) ), std::invalid_argument,
00151                       ">>> ERROR: (Intrepid::getValues_HGRAD_Args) DIV is invalid operator for rank-0 (scalar) fields in 2D."); 
00152   
00153   TEST_FOR_EXCEPTION( ( (spaceDim == 3) && ( (operatorType == OPERATOR_DIV) ||
00154                                              (operatorType == OPERATOR_CURL) ) ), std::invalid_argument,
00155                       ">>> ERROR: (Intrepid::getValues_HGRAD_Args) DIV and CURL are invalid operators for rank-0 (scalar) fields in 3D."); 
00156   
00157   
00158   // Check rank of outputValues (all operators are admissible in 1D) and its dim 2 when operator is
00159   // GRAD, CURL (only in 2D), or Dk.
00160   
00161   if(spaceDim == 1) {
00162     switch(operatorType){
00163       case OPERATOR_VALUE:
00164         TEST_FOR_EXCEPTION( !(outputValues.rank() == 2), std::invalid_argument,
00165                             ">>> ERROR: (Intrepid::getValues_HGRAD_Args) rank = 2 required for outputValues when operator = VALUE.");
00166         break;
00167       case OPERATOR_GRAD:
00168       case OPERATOR_CURL:
00169       case OPERATOR_DIV:
00170       case OPERATOR_D1:
00171       case OPERATOR_D2:
00172       case OPERATOR_D3:
00173       case OPERATOR_D4:
00174       case OPERATOR_D5:
00175       case OPERATOR_D6:
00176       case OPERATOR_D7:
00177       case OPERATOR_D8:
00178       case OPERATOR_D9:
00179       case OPERATOR_D10:
00180         TEST_FOR_EXCEPTION( !(outputValues.rank() == 3), std::invalid_argument,
00181                             ">>> ERROR: (Intrepid::getValues_HGRAD_Args) rank = 3 required for outputValues in 1D when operator = GRAD, CURL, DIV, or Dk.");
00182         
00183         TEST_FOR_EXCEPTION( !(outputValues.dimension(2) == 1 ),
00184                             std::invalid_argument,
00185                             ">>> ERROR: (Intrepid::getValues_HGRAD_Args) dim 2 of outputValues must equal 1 when operator = GRAD, CURL, DIV, or Dk.");
00186         break;
00187       default:
00188         TEST_FOR_EXCEPTION( (true), std::invalid_argument, ">>> ERROR: (Intrepid::getValues_HGRAD_Args) Invalid operator");
00189     }
00190   }
00191   else if(spaceDim > 1) {
00192     switch(operatorType){
00193       case OPERATOR_VALUE:
00194         TEST_FOR_EXCEPTION( !(outputValues.rank() == 2), std::invalid_argument,
00195                             ">>> ERROR: (Intrepid::getValues_HGRAD_Args) rank = 2 required for outputValues when operator = VALUE.");
00196         break;
00197       case OPERATOR_GRAD:
00198       case OPERATOR_CURL:
00199       case OPERATOR_D1:
00200         TEST_FOR_EXCEPTION( !(outputValues.rank() == 3), std::invalid_argument,
00201                             ">>> ERROR: (Intrepid::getValues_HGRAD_Args) rank = 3 required for outputValues in 2D and 3D when operator = GRAD, CURL (in 2D), or Dk.");
00202         
00203         TEST_FOR_EXCEPTION( !(outputValues.dimension(2) == spaceDim ),
00204                             std::invalid_argument,
00205                             ">>> ERROR: (Intrepid::getValues_HGRAD_Args) dim 2 of outputValues must equal cell dimension when operator = GRAD, CURL (in 2D), or D1.");
00206         break;
00207       case OPERATOR_D2:
00208       case OPERATOR_D3:
00209       case OPERATOR_D4:
00210       case OPERATOR_D5:
00211       case OPERATOR_D6:
00212       case OPERATOR_D7:
00213       case OPERATOR_D8:
00214       case OPERATOR_D9:
00215       case OPERATOR_D10:
00216         TEST_FOR_EXCEPTION( !(outputValues.rank() == 3), std::invalid_argument,
00217                             ">>> ERROR: (Intrepid::getValues_HGRAD_Args) rank = 3 required for outputValues in 2D and 3D when operator = GRAD, CURL (in 2D), or Dk.");
00218         
00219         TEST_FOR_EXCEPTION( !(outputValues.dimension(2) == Intrepid::getDkCardinality(operatorType, spaceDim) ),
00220                             std::invalid_argument,
00221                             ">>> ERROR: (Intrepid::getValues_HGRAD_Args) dim 2 of outputValues must equal cardinality of the Dk multiset.");
00222         break;
00223       default:
00224         TEST_FOR_EXCEPTION( (true), std::invalid_argument, ">>> ERROR: (Intrepid::getValues_HGRAD_Args) Invalid operator");
00225     }
00226   }
00227 
00228   
00229   // Verify dim 0 and dim 1 of outputValues
00230   TEST_FOR_EXCEPTION( !(outputValues.dimension(1) == inputPoints.dimension(0) ), 
00231                       std::invalid_argument, 
00232                       ">>> ERROR: (Intrepid::getValues_HGRAD_Args) dim 1 (number of points) of outputValues must equal dim 0 of inputPoints.");
00233   
00234   TEST_FOR_EXCEPTION( !(outputValues.dimension(0) == basisCard ),
00235                       std::invalid_argument,
00236                       ">>> ERROR: (Intrepid::getValues_HGRAD_Args) dim 0 (number of basis functions) of outputValues must equal basis cardinality.");
00237 }
00238 
00239 
00240 
00241 template<class Scalar, class ArrayScalar>
00242 void getValues_HCURL_Args(ArrayScalar &                outputValues,
00243                           const ArrayScalar &          inputPoints,
00244                           const EOperator              operatorType,
00245                           const shards::CellTopology&  cellTopo,
00246                           const int                    basisCard){
00247   
00248   int spaceDim = cellTopo.getDimension();
00249   
00250   // Verify that cell is 2D or 3D (this is redundant for default bases where we use correct cells,
00251   //  but will force any user-defined basis for HCURL spaces to use 2D or 3D cells
00252   TEST_FOR_EXCEPTION( !( (spaceDim == 2) || (spaceDim == 3) ), std::invalid_argument, 
00253                       ">>> ERROR: (Intrepid::getValues_HCURL_Args) cell dimension = 2 or 3 required for HCURL basis"); 
00254 
00255   
00256   // Verify inputPoints array
00257   TEST_FOR_EXCEPTION( !(inputPoints.rank() == 2), std::invalid_argument, 
00258                       ">>> ERROR: (Intrepid::getValues_HCURL_Args) rank = 2 required for inputPoints array"); 
00259   TEST_FOR_EXCEPTION(  (inputPoints.dimension(0) <= 0), std::invalid_argument,
00260                       ">>> ERROR (Intrepid::getValues_HCURL_Args): dim 0 (number of points) > 0 required for inputPoints array");
00261 
00262   TEST_FOR_EXCEPTION( !(inputPoints.dimension(1) == spaceDim), std::invalid_argument,
00263                       ">>> ERROR: (Intrepid::getValues_HCURL_Args) dim 1 (spatial dimension) of inputPoints array  does not match cell dimension");
00264   
00265   // Verify that all inputPoints are in the reference cell
00266   /*
00267   TEST_FOR_EXCEPTION( !CellTools<Scalar>::checkPointSetInclusion(inputPoints, cellTopo), std::invalid_argument,
00268                       ">>> ERROR: (Intrepid::getValues_HCURL_Args) One or more points are outside the " 
00269                       << cellTopo <<" reference cell");
00270    */
00271   
00272   // Verify that operatorType is admissible for HCURL fields
00273   TEST_FOR_EXCEPTION( !( (operatorType == OPERATOR_VALUE) || (operatorType == OPERATOR_CURL) ), std::invalid_argument,
00274                       ">>> ERROR: (Intrepid::getValues_HCURL_Args) operator = VALUE or CURL required for HCURL fields."); 
00275   
00276   
00277   // Check rank of outputValues: for VALUE should always be rank-3 array with (F,P,D) layout 
00278   switch(operatorType) {
00279     
00280     case OPERATOR_VALUE:
00281       TEST_FOR_EXCEPTION( !(outputValues.rank() == 3), std::invalid_argument,
00282                           ">>> ERROR: (Intrepid::getValues_HCURL_Args) rank = 3 required for outputValues when operator is VALUE");
00283       TEST_FOR_EXCEPTION( !(outputValues.dimension(2) == spaceDim ),
00284                           std::invalid_argument,
00285                           ">>> ERROR: (Intrepid::getValues_HCURL_Args) dim 2 of outputValues must equal cell dimension when operator is VALUE.");
00286       break;
00287       
00288     case OPERATOR_CURL:
00289       
00290       // in 3D we need an (F,P,D) container because CURL gives a vector field:
00291       if(spaceDim == 3) {
00292         TEST_FOR_EXCEPTION( !(outputValues.rank() == 3 ) ,
00293                             std::invalid_argument,
00294                             ">>> ERROR: (Intrepid::getValues_HCURL_Args) rank = 3 required for outputValues in 3D when operator is CURL");
00295         TEST_FOR_EXCEPTION( !(outputValues.dimension(2) == spaceDim),
00296                             std::invalid_argument,
00297                             ">>> ERROR: (Intrepid::getValues_HCURL_Args) dim 2 of outputValues must equal cell dimension in 3D when operator is CURL.");
00298       }
00299       // In 2D we need an (F,P) container because CURL gives a scalar field
00300       else if(spaceDim == 2) {
00301         TEST_FOR_EXCEPTION( !(outputValues.rank() == 2 ) ,
00302                             std::invalid_argument,
00303                             ">>> ERROR: (Intrepid::getValues_HCURL_Args) rank = 2 required for outputValues in 2D when operator is CURL");
00304       }
00305       break;
00306       
00307     default:
00308       TEST_FOR_EXCEPTION( (true), std::invalid_argument, ">>> ERROR: (Intrepid::getValues_HCURL_Args) Invalid operator");
00309   }
00310   
00311   
00312   // Verify dim 0 and dim 1 of outputValues
00313   TEST_FOR_EXCEPTION( !(outputValues.dimension(1) == inputPoints.dimension(0) ), 
00314                       std::invalid_argument, 
00315                       ">>> ERROR: (Intrepid::getValues_HCURL_Args) dim 1 (number of points) of outputValues must equal dim 0 of inputPoints.");
00316   
00317   TEST_FOR_EXCEPTION( !(outputValues.dimension(0) == basisCard ),
00318                       std::invalid_argument,
00319                       ">>> ERROR: (Intrepid::getValues_HCURL_Args) dim 0 (number of basis functions) of outputValues must equal basis cardinality.");
00320 
00321 }
00322 
00323 
00324 
00325 template<class Scalar, class ArrayScalar>
00326 void getValues_HDIV_Args(ArrayScalar &                outputValues,
00327                           const ArrayScalar &          inputPoints,
00328                           const EOperator              operatorType,
00329                           const shards::CellTopology&  cellTopo,
00330                           const int                    basisCard){
00331   
00332   int spaceDim = cellTopo.getDimension();
00333   
00334   // Verify inputPoints array
00335   TEST_FOR_EXCEPTION( !(inputPoints.rank() == 2), std::invalid_argument, 
00336                       ">>> ERROR: (Intrepid::getValues_HDIV_Args) rank = 2 required for inputPoints array"); 
00337   TEST_FOR_EXCEPTION(  (inputPoints.dimension(0) <= 0), std::invalid_argument,
00338                       ">>> ERROR (Intrepid::getValues_HDIV_Args): dim 0 (number of points) > 0 required for inputPoints array");
00339 
00340   TEST_FOR_EXCEPTION( !(inputPoints.dimension(1) == spaceDim), std::invalid_argument,
00341                       ">>> ERROR: (Intrepid::getValues_HDIV_Args) dim 1 (spatial dimension) of inputPoints array  does not match cell dimension");
00342   
00343   // Verify that all inputPoints are in the reference cell
00344   /*
00345   TEST_FOR_EXCEPTION( !CellTools<Scalar>::checkPointSetInclusion(inputPoints, cellTopo), std::invalid_argument,
00346                       ">>> ERROR: (Intrepid::getValues_HDIV_Args) One or more points are outside the " 
00347                       << cellTopo <<" reference cell");
00348    */
00349   
00350   // Verify that operatorType is admissible for HDIV fields
00351   TEST_FOR_EXCEPTION( !( (operatorType == OPERATOR_VALUE) || (operatorType == OPERATOR_DIV) ), std::invalid_argument,
00352                       ">>> ERROR: (Intrepid::getValues_HDIV_Args) operator = VALUE or DIV required for HDIV fields."); 
00353   
00354   
00355   // Check rank of outputValues 
00356   switch(operatorType) {
00357     case OPERATOR_VALUE:
00358       TEST_FOR_EXCEPTION( !(outputValues.rank() == 3), std::invalid_argument,
00359                           ">>> ERROR: (Intrepid::getValues_HDIV_Args) rank = 3 required for outputValues when operator is VALUE.");
00360       
00361       TEST_FOR_EXCEPTION( !(outputValues.dimension(2) == spaceDim ),
00362                           std::invalid_argument,
00363                           ">>> ERROR: (Intrepid::getValues_HDIV_Args) dim 2 of outputValues must equal cell dimension for operator VALUE.");
00364       break;
00365     case OPERATOR_DIV:
00366       TEST_FOR_EXCEPTION( !(outputValues.rank() == 2), std::invalid_argument,
00367                           ">>> ERROR: (Intrepid::getValues_HDIV_Args) rank = 2 required for outputValues when operator is DIV.");
00368       break;
00369 
00370     default:
00371       TEST_FOR_EXCEPTION( (true), std::invalid_argument, ">>> ERROR: (Intrepid::getValues_HDIV_Args) Invalid operator");
00372   }
00373 
00374   
00375   // Verify dim 0 and dim 1 of outputValues
00376   TEST_FOR_EXCEPTION( !(outputValues.dimension(1) == inputPoints.dimension(0) ), 
00377                       std::invalid_argument, 
00378                       ">>> ERROR: (Intrepid::getValues_HDIV_Args) dim 1 (number of points) of outputValues must equal dim 0 of inputPoints.");
00379   
00380   TEST_FOR_EXCEPTION( !(outputValues.dimension(0) == basisCard ),
00381                       std::invalid_argument,
00382                       ">>> ERROR: (Intrepid::getValues_HDIV_Args) dim 0 (number of basis functions) of outputValues must equal basis cardinality.");
00383 }
http://trilinos.sandia.gov/packages/docs/r10.0/packages/intrepid/src/Discretization/Basis/Intrepid_BasisDef.hpp
