|
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 // 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() {}
1.7.6.1