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