|
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 } 00384 00385 // Pure virtual destructor (gives warnings if not included). 00386 // Following "Effective C++: 3rd Ed." item 7 the implementation 00387 // is included in the definition file. 00388 template<class ArrayScalar> 00389 DofCoordsInterface<ArrayScalar>::~DofCoordsInterface() {}
1.7.4