|
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 00049 namespace Intrepid { 00050 00051 template<class Scalar, class ArrayScalar> 00052 Basis_HCURL_HEX_I1_FEM<Scalar,ArrayScalar>::Basis_HCURL_HEX_I1_FEM() 00053 { 00054 this -> basisCardinality_ = 12; 00055 this -> basisDegree_ = 1; 00056 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >() ); 00057 this -> basisType_ = BASIS_FEM_DEFAULT; 00058 this -> basisCoordinates_ = COORDINATES_CARTESIAN; 00059 this -> basisTagsAreSet_ = false; 00060 } 00061 00062 template<class Scalar, class ArrayScalar> 00063 void Basis_HCURL_HEX_I1_FEM<Scalar, ArrayScalar>::initializeTags() { 00064 00065 // Basis-dependent intializations 00066 int tagSize = 4; // size of DoF tag 00067 int posScDim = 0; // position in the tag, counting from 0, of the subcell dim 00068 int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal 00069 int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell 00070 00071 // An array with local DoF tags assigned to basis functions, in the order of their local enumeration 00072 int tags[] = { 00073 1, 0, 0, 1, 00074 1, 1, 0, 1, 00075 1, 2, 0, 1, 00076 1, 3, 0, 1, 00077 1, 4, 0, 1, 00078 1, 5, 0, 1, 00079 1, 6, 0, 1, 00080 1, 7, 0, 1, 00081 1, 8, 0, 1, 00082 1, 9, 0, 1, 00083 1, 10, 0, 1, 00084 1, 11, 0, 1 }; 00085 00086 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays: 00087 Intrepid::setOrdinalTagData(this -> tagToOrdinal_, 00088 this -> ordinalToTag_, 00089 tags, 00090 this -> basisCardinality_, 00091 tagSize, 00092 posScDim, 00093 posScOrd, 00094 posDfOrd); 00095 } 00096 00097 00098 00099 template<class Scalar, class ArrayScalar> 00100 void Basis_HCURL_HEX_I1_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar & outputValues, 00101 const ArrayScalar & inputPoints, 00102 const EOperator operatorType) const { 00103 00104 // Verify arguments 00105 #ifdef HAVE_INTREPID_DEBUG 00106 Intrepid::getValues_HCURL_Args<Scalar, ArrayScalar>(outputValues, 00107 inputPoints, 00108 operatorType, 00109 this -> getBaseCellTopology(), 00110 this -> getCardinality() ); 00111 #endif 00112 00113 // Number of evaluation points = dim 0 of inputPoints 00114 int dim0 = inputPoints.dimension(0); 00115 00116 // Temporaries: (x,y,z) coordinates of the evaluation point 00117 Scalar x = 0.0; 00118 Scalar y = 0.0; 00119 Scalar z = 0.0; 00120 00121 switch (operatorType) { 00122 case OPERATOR_VALUE: 00123 for (int i0 = 0; i0 < dim0; i0++) { 00124 x = inputPoints(i0, 0); 00125 y = inputPoints(i0, 1); 00126 z = inputPoints(i0, 2); 00127 00128 // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim) 00129 outputValues(0, i0, 0) = (1.0 - y)*(1.0 - z)/8.0; 00130 outputValues(0, i0, 1) = 0.0; 00131 outputValues(0, i0, 2) = 0.0; 00132 00133 outputValues(1, i0, 0) = 0.0; 00134 outputValues(1, i0, 1) = (1.0 + x)*(1.0 - z)/8.0; 00135 outputValues(1, i0, 2) = 0.0; 00136 00137 outputValues(2, i0, 0) = -(1.0 + y)*(1.0 - z)/8.0; 00138 outputValues(2, i0, 1) = 0.0; 00139 outputValues(2, i0, 2) = 0.0; 00140 00141 outputValues(3, i0, 0) = 0.0; 00142 outputValues(3, i0, 1) = -(1.0 - x)*(1.0 - z)/8.0; 00143 outputValues(3, i0, 2) = 0.0; 00144 00145 outputValues(4, i0, 0) = (1.0 - y)*(1.0 + z)/8.0; 00146 outputValues(4, i0, 1) = 0.0; 00147 outputValues(4, i0, 2) = 0.0; 00148 00149 outputValues(5, i0, 0) = 0.0; 00150 outputValues(5, i0, 1) = (1.0 + x)*(1.0 + z)/8.0; 00151 outputValues(5, i0, 2) = 0.0; 00152 00153 outputValues(6, i0, 0) = -(1.0 + y)*(1.0 + z)/8.0; 00154 outputValues(6, i0, 1) = 0.0; 00155 outputValues(6, i0, 2) = 0.0; 00156 00157 outputValues(7, i0, 0) = 0.0; 00158 outputValues(7, i0, 1) = -(1.0 - x)*(1.0 + z)/8.0; 00159 outputValues(7, i0, 2) = 0.0; 00160 00161 outputValues(8, i0, 0) = 0.0; 00162 outputValues(8, i0, 1) = 0.0; 00163 outputValues(8, i0, 2) = (1.0 - x)*(1.0 - y)/8.0; 00164 00165 outputValues(9, i0, 0) = 0.0; 00166 outputValues(9, i0, 1) = 0.0; 00167 outputValues(9, i0, 2) = (1.0 + x)*(1.0 - y)/8.0; 00168 00169 outputValues(10, i0, 0) = 0.0; 00170 outputValues(10, i0, 1) = 0.0; 00171 outputValues(10, i0, 2) = (1.0 + x)*(1.0 + y)/8.0; 00172 00173 outputValues(11, i0, 0) = 0.0; 00174 outputValues(11, i0, 1) = 0.0; 00175 outputValues(11, i0, 2) = (1.0 - x)*(1.0 + y)/8.0; 00176 } 00177 break; 00178 00179 case OPERATOR_CURL: 00180 for (int i0 = 0; i0 < dim0; i0++) { 00181 x = inputPoints(i0, 0); 00182 y = inputPoints(i0, 1); 00183 z = inputPoints(i0, 2); 00184 00185 // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim) 00186 outputValues(0, i0, 0) = 0.0; 00187 outputValues(0, i0, 1) = -(1.0 - y)/8.0; 00188 outputValues(0, i0, 2) = (1.0 - z)/8.0; 00189 00190 outputValues(1, i0, 0) = (1.0 + x)/8.0; 00191 outputValues(1, i0, 1) = 0.0; 00192 outputValues(1, i0, 2) = (1.0 - z)/8.0; 00193 00194 outputValues(2, i0, 0) = 0.0; 00195 outputValues(2, i0, 1) = (1.0 + y)/8.0; 00196 outputValues(2, i0, 2) = (1.0 - z)/8.0; 00197 00198 outputValues(3, i0, 0) = -(1.0 - x)/8.0; 00199 outputValues(3, i0, 1) = 0.0; 00200 outputValues(3, i0, 2) = (1.0 - z)/8.0; 00201 00202 outputValues(4, i0, 0) = 0.0; 00203 outputValues(4, i0, 1) = (1.0 - y)/8.0; 00204 outputValues(4, i0, 2) = (1.0 + z)/8.0; 00205 00206 outputValues(5, i0, 0) = -(1.0 + x)/8.0; 00207 outputValues(5, i0, 1) = 0.0; 00208 outputValues(5, i0, 2) = (1.0 + z)/8.0; 00209 00210 outputValues(6, i0, 0) = 0.0; 00211 outputValues(6, i0, 1) = -(1.0 + y)/8.0; 00212 outputValues(6, i0, 2) = (1.0 + z)/8.0; 00213 00214 outputValues(7, i0, 0) = (1.0 - x)/8.0; 00215 outputValues(7, i0, 1) = 0.0; 00216 outputValues(7, i0, 2) = (1.0 + z)/8.0; 00217 00218 outputValues(8, i0, 0) = -(1.0 - x)/8.0; 00219 outputValues(8, i0, 1) = (1.0 - y)/8.0; 00220 outputValues(8, i0, 2) = 0.0; 00221 00222 outputValues(9, i0, 0) = -(1.0 + x)/8.0; 00223 outputValues(9, i0, 1) = -(1.0 - y)/8.0; 00224 outputValues(9, i0, 2) = 0.0; 00225 00226 outputValues(10, i0, 0) = (1.0 + x)/8.0; 00227 outputValues(10, i0, 1) = -(1.0 + y)/8.0; 00228 outputValues(10, i0, 2) = 0.0; 00229 00230 outputValues(11, i0, 0) = (1.0 - x)/8.0; 00231 outputValues(11, i0, 1) = (1.0 + y)/8.0; 00232 outputValues(11, i0, 2) = 0.0; 00233 } 00234 break; 00235 00236 case OPERATOR_DIV: 00237 TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument, 00238 ">>> ERROR (Basis_HCURL_HEX_I1_FEM): DIV is invalid operator for HCURL Basis Functions"); 00239 break; 00240 00241 case OPERATOR_GRAD: 00242 TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_GRAD), std::invalid_argument, 00243 ">>> ERROR (Basis_HCURL_HEX_I1_FEM): GRAD is invalid operator for HCURL Basis Functions"); 00244 break; 00245 00246 case OPERATOR_D1: 00247 case OPERATOR_D2: 00248 case OPERATOR_D3: 00249 case OPERATOR_D4: 00250 case OPERATOR_D5: 00251 case OPERATOR_D6: 00252 case OPERATOR_D7: 00253 case OPERATOR_D8: 00254 case OPERATOR_D9: 00255 case OPERATOR_D10: 00256 TEUCHOS_TEST_FOR_EXCEPTION( ( (operatorType == OPERATOR_D1) || 00257 (operatorType == OPERATOR_D2) || 00258 (operatorType == OPERATOR_D3) || 00259 (operatorType == OPERATOR_D4) || 00260 (operatorType == OPERATOR_D5) || 00261 (operatorType == OPERATOR_D6) || 00262 (operatorType == OPERATOR_D7) || 00263 (operatorType == OPERATOR_D8) || 00264 (operatorType == OPERATOR_D9) || 00265 (operatorType == OPERATOR_D10) ), 00266 std::invalid_argument, 00267 ">>> ERROR (Basis_HCURL_HEX_I1_FEM): Invalid operator type"); 00268 break; 00269 00270 default: 00271 TEUCHOS_TEST_FOR_EXCEPTION( ( (operatorType != OPERATOR_VALUE) && 00272 (operatorType != OPERATOR_GRAD) && 00273 (operatorType != OPERATOR_CURL) && 00274 (operatorType != OPERATOR_DIV) && 00275 (operatorType != OPERATOR_D1) && 00276 (operatorType != OPERATOR_D2) && 00277 (operatorType != OPERATOR_D3) && 00278 (operatorType != OPERATOR_D4) && 00279 (operatorType != OPERATOR_D5) && 00280 (operatorType != OPERATOR_D6) && 00281 (operatorType != OPERATOR_D7) && 00282 (operatorType != OPERATOR_D8) && 00283 (operatorType != OPERATOR_D9) && 00284 (operatorType != OPERATOR_D10) ), 00285 std::invalid_argument, 00286 ">>> ERROR (Basis_HCURL_HEX_I1_FEM): Invalid operator type"); 00287 } 00288 } 00289 00290 00291 00292 template<class Scalar, class ArrayScalar> 00293 void Basis_HCURL_HEX_I1_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar& outputValues, 00294 const ArrayScalar & inputPoints, 00295 const ArrayScalar & cellVertices, 00296 const EOperator operatorType) const { 00297 TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error, 00298 ">>> ERROR (Basis_HCURL_HEX_I1_FEM): FEM Basis calling an FVD member function"); 00299 } 00300 00301 template<class Scalar, class ArrayScalar> 00302 void Basis_HCURL_HEX_I1_FEM<Scalar, ArrayScalar>::getDofCoords(ArrayScalar & DofCoords) const { 00303 #ifdef HAVE_INTREPID_DEBUG 00304 // Verify rank of output array. 00305 TEUCHOS_TEST_FOR_EXCEPTION( !(DofCoords.rank() == 2), std::invalid_argument, 00306 ">>> ERROR: (Intrepid::Basis_HCURL_QUAD_I1_FEM::getDofCoords) rank = 2 required for DofCoords array"); 00307 // Verify 0th dimension of output array. 00308 TEUCHOS_TEST_FOR_EXCEPTION( !( DofCoords.dimension(0) == this -> basisCardinality_ ), std::invalid_argument, 00309 ">>> ERROR: (Intrepid::Basis_HCURL_QUAD_I1_FEM::getDofCoords) mismatch in number of DoF and 0th dimension of DofCoords array"); 00310 // Verify 1st dimension of output array. 00311 TEUCHOS_TEST_FOR_EXCEPTION( !( DofCoords.dimension(1) == (int)(this -> basisCellTopology_.getDimension()) ), std::invalid_argument, 00312 ">>> ERROR: (Intrepid::Basis_HCURL_QUAD_I1_FEM::getDofCoords) incorrect reference cell (1st) dimension in DofCoords array"); 00313 #endif 00314 00315 DofCoords(0,0) = 0.0; DofCoords(0,1) = -1.0; DofCoords(0,2) = -1.0; 00316 DofCoords(1,0) = 1.0; DofCoords(1,1) = 0.0; DofCoords(1,2) = -1.0; 00317 DofCoords(2,0) = 0.0; DofCoords(2,1) = 1.0; DofCoords(2,2) = -1.0; 00318 DofCoords(3,0) = -1.0; DofCoords(3,1) = 0.0; DofCoords(3,2) = -1.0; 00319 DofCoords(4,0) = 0.0; DofCoords(4,1) = -1.0; DofCoords(4,2) = 1.0; 00320 DofCoords(5,0) = 1.0; DofCoords(5,1) = 0.0; DofCoords(5,2) = 1.0; 00321 DofCoords(6,0) = 0.0; DofCoords(6,1) = 1.0; DofCoords(6,2) = 1.0; 00322 DofCoords(7,0) = -1.0; DofCoords(7,1) = 0.0; DofCoords(7,2) = 1.0; 00323 DofCoords(8,0) = -1.0; DofCoords(8,1) = -1.0; DofCoords(8,2) = 0.0; 00324 DofCoords(9,0) = 1.0; DofCoords(9,1) = -1.0; DofCoords(9,2) = 0.0; 00325 DofCoords(10,0) = 1.0; DofCoords(10,1) = 1.0; DofCoords(10,2) = 0.0; 00326 DofCoords(11,0) = -1.0; DofCoords(11,1) = 1.0; DofCoords(11,2) = 0.0; 00327 } 00328 00329 }// namespace Intrepid
1.7.6.1