|
Intrepid
|
00001 #ifndef INTREPID_HGRAD_HEX_C1_FEMDEF_HPP 00002 #define INTREPID_HGRAD_HEX_C1_FEMDEF_HPP 00003 // @HEADER 00004 // ************************************************************************ 00005 // 00006 // Intrepid Package 00007 // Copyright (2007) Sandia Corporation 00008 // 00009 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive 00010 // license for use of this work by or on behalf of the U.S. Government. 00011 // 00012 // Redistribution and use in source and binary forms, with or without 00013 // modification, are permitted provided that the following conditions are 00014 // met: 00015 // 00016 // 1. Redistributions of source code must retain the above copyright 00017 // notice, this list of conditions and the following disclaimer. 00018 // 00019 // 2. Redistributions in binary form must reproduce the above copyright 00020 // notice, this list of conditions and the following disclaimer in the 00021 // documentation and/or other materials provided with the distribution. 00022 // 00023 // 3. Neither the name of the Corporation nor the names of the 00024 // contributors may be used to endorse or promote products derived from 00025 // this software without specific prior written permission. 00026 // 00027 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY 00028 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 00029 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 00030 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE 00031 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 00032 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 00033 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 00034 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 00035 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 00036 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 00037 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00038 // 00039 // Questions? Contact Pavel Bochev (pbboche@sandia.gov) 00040 // Denis Ridzal (dridzal@sandia.gov), or 00041 // Kara Peterson (kjpeter@sandia.gov) 00042 // 00043 // ************************************************************************ 00044 // @HEADER 00045 00051 namespace Intrepid { 00052 00053 00054 template<class Scalar, class ArrayScalar> 00055 Basis_HGRAD_HEX_C1_FEM<Scalar,ArrayScalar>::Basis_HGRAD_HEX_C1_FEM() 00056 { 00057 this -> basisCardinality_ = 8; 00058 this -> basisDegree_ = 1; 00059 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >() ); 00060 this -> basisType_ = BASIS_FEM_DEFAULT; 00061 this -> basisCoordinates_ = COORDINATES_CARTESIAN; 00062 this -> basisTagsAreSet_ = false; 00063 } 00064 00065 00066 00067 template<class Scalar, class ArrayScalar> 00068 void Basis_HGRAD_HEX_C1_FEM<Scalar, ArrayScalar>::initializeTags() { 00069 00070 // Basis-dependent intializations 00071 int tagSize = 4; // size of DoF tag, i.e., number of fields in the tag 00072 int posScDim = 0; // position in the tag, counting from 0, of the subcell dim 00073 int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal 00074 int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell 00075 00076 // An array with local DoF tags assigned to basis functions, in the order of their local enumeration 00077 int tags[] = { 0, 0, 0, 1, 00078 0, 1, 0, 1, 00079 0, 2, 0, 1, 00080 0, 3, 0, 1, 00081 0, 4, 0, 1, 00082 0, 5, 0, 1, 00083 0, 6, 0, 1, 00084 0, 7, 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_HGRAD_HEX_C1_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_HGRAD_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 00123 case OPERATOR_VALUE: 00124 for (int i0 = 0; i0 < dim0; i0++) { 00125 x = inputPoints(i0, 0); 00126 y = inputPoints(i0, 1); 00127 z = inputPoints(i0, 2); 00128 00129 // outputValues is a rank-2 array with dimensions (basisCardinality_, dim0) 00130 outputValues(0, i0) = (1.0 - x)*(1.0 - y)*(1.0 - z)/8.0; 00131 outputValues(1, i0) = (1.0 + x)*(1.0 - y)*(1.0 - z)/8.0; 00132 outputValues(2, i0) = (1.0 + x)*(1.0 + y)*(1.0 - z)/8.0; 00133 outputValues(3, i0) = (1.0 - x)*(1.0 + y)*(1.0 - z)/8.0; 00134 00135 outputValues(4, i0) = (1.0 - x)*(1.0 - y)*(1.0 + z)/8.0; 00136 outputValues(5, i0) = (1.0 + x)*(1.0 - y)*(1.0 + z)/8.0; 00137 outputValues(6, i0) = (1.0 + x)*(1.0 + y)*(1.0 + z)/8.0; 00138 outputValues(7, i0) = (1.0 - x)*(1.0 + y)*(1.0 + z)/8.0; 00139 } 00140 break; 00141 00142 case OPERATOR_GRAD: 00143 case OPERATOR_D1: 00144 for (int i0 = 0; i0 < dim0; i0++) { 00145 x = inputPoints(i0,0); 00146 y = inputPoints(i0,1); 00147 z = inputPoints(i0,2); 00148 00149 // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim) 00150 outputValues(0, i0, 0) = -(1.0 - y)*(1.0 - z)/8.0; 00151 outputValues(0, i0, 1) = -(1.0 - x)*(1.0 - z)/8.0; 00152 outputValues(0, i0, 2) = -(1.0 - x)*(1.0 - y)/8.0; 00153 00154 outputValues(1, i0, 0) = (1.0 - y)*(1.0 - z)/8.0; 00155 outputValues(1, i0, 1) = -(1.0 + x)*(1.0 - z)/8.0; 00156 outputValues(1, i0, 2) = -(1.0 + x)*(1.0 - y)/8.0; 00157 00158 outputValues(2, i0, 0) = (1.0 + y)*(1.0 - z)/8.0; 00159 outputValues(2, i0, 1) = (1.0 + x)*(1.0 - z)/8.0; 00160 outputValues(2, i0, 2) = -(1.0 + x)*(1.0 + y)/8.0; 00161 00162 outputValues(3, i0, 0) = -(1.0 + y)*(1.0 - z)/8.0; 00163 outputValues(3, i0, 1) = (1.0 - x)*(1.0 - z)/8.0; 00164 outputValues(3, i0, 2) = -(1.0 - x)*(1.0 + y)/8.0; 00165 00166 outputValues(4, i0, 0) = -(1.0 - y)*(1.0 + z)/8.0; 00167 outputValues(4, i0, 1) = -(1.0 - x)*(1.0 + z)/8.0; 00168 outputValues(4, i0, 2) = (1.0 - x)*(1.0 - y)/8.0; 00169 00170 outputValues(5, i0, 0) = (1.0 - y)*(1.0 + z)/8.0; 00171 outputValues(5, i0, 1) = -(1.0 + x)*(1.0 + z)/8.0; 00172 outputValues(5, i0, 2) = (1.0 + x)*(1.0 - y)/8.0; 00173 00174 outputValues(6, i0, 0) = (1.0 + y)*(1.0 + z)/8.0; 00175 outputValues(6, i0, 1) = (1.0 + x)*(1.0 + z)/8.0; 00176 outputValues(6, i0, 2) = (1.0 + x)*(1.0 + y)/8.0; 00177 00178 outputValues(7, i0, 0) = -(1.0 + y)*(1.0 + z)/8.0; 00179 outputValues(7, i0, 1) = (1.0 - x)*(1.0 + z)/8.0; 00180 outputValues(7, i0, 2) = (1.0 - x)*(1.0 + y)/8.0; 00181 } 00182 break; 00183 00184 case OPERATOR_CURL: 00185 TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_CURL), std::invalid_argument, 00186 ">>> ERROR (Basis_HGRAD_HEX_C1_FEM): CURL is invalid operator for rank-0 (scalar) functions in 3D"); 00187 break; 00188 00189 case OPERATOR_DIV: 00190 TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument, 00191 ">>> ERROR (Basis_HGRAD_HEX_C1_FEM): DIV is invalid operator for rank-0 (scalar) functions in 3D"); 00192 break; 00193 00194 case OPERATOR_D2: 00195 for (int i0 = 0; i0 < dim0; i0++) { 00196 x = inputPoints(i0,0); 00197 y = inputPoints(i0,1); 00198 z = inputPoints(i0,2); 00199 00200 // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, D2Cardinality = 6) 00201 outputValues(0, i0, 0) = 0.0; // {2, 0, 0} 00202 outputValues(0, i0, 1) = (1.0 - z)/8.0; // {1, 1, 0} 00203 outputValues(0, i0, 2) = (1.0 - y)/8.0; // {1, 0, 1} 00204 outputValues(0, i0, 3) = 0.0; // {0, 2, 0} 00205 outputValues(0, i0, 4) = (1.0 - x)/8.0; // {0, 1, 1} 00206 outputValues(0, i0, 5) = 0.0; // {0, 0, 2} 00207 00208 outputValues(1, i0, 0) = 0.0; // {2, 0, 0} 00209 outputValues(1, i0, 1) = -(1.0 - z)/8.0; // {1, 1, 0} 00210 outputValues(1, i0, 2) = -(1.0 - y)/8.0; // {1, 0, 1} 00211 outputValues(1, i0, 3) = 0.0; // {0, 2, 0} 00212 outputValues(1, i0, 4) = (1.0 + x)/8.0; // {0, 1, 1} 00213 outputValues(1, i0, 5) = 0.0; // {0, 0, 2} 00214 00215 outputValues(2, i0, 0) = 0.0; // {2, 0, 0} 00216 outputValues(2, i0, 1) = (1.0 - z)/8.0; // {1, 1, 0} 00217 outputValues(2, i0, 2) = -(1.0 + y)/8.0; // {1, 0, 1} 00218 outputValues(2, i0, 3) = 0.0; // {0, 2, 0} 00219 outputValues(2, i0, 4) = -(1.0 + x)/8.0; // {0, 1, 1} 00220 outputValues(2, i0, 5) = 0.0; // {0, 0, 2} 00221 00222 outputValues(3, i0, 0) = 0.0; // {2, 0, 0} 00223 outputValues(3, i0, 1) = -(1.0 - z)/8.0; // {1, 1, 0} 00224 outputValues(3, i0, 2) = (1.0 + y)/8.0; // {1, 0, 1} 00225 outputValues(3, i0, 3) = 0.0; // {0, 2, 0} 00226 outputValues(3, i0, 4) = -(1.0 - x)/8.0; // {0, 1, 1} 00227 outputValues(3, i0, 5) = 0.0; // {0, 0, 2} 00228 00229 00230 outputValues(4, i0, 0) = 0.0; // {2, 0, 0} 00231 outputValues(4, i0, 1) = (1.0 + z)/8.0; // {1, 1, 0} 00232 outputValues(4, i0, 2) = -(1.0 - y)/8.0; // {1, 0, 1} 00233 outputValues(4, i0, 3) = 0.0; // {0, 2, 0} 00234 outputValues(4, i0, 4) = -(1.0 - x)/8.0; // {0, 1, 1} 00235 outputValues(4, i0, 5) = 0.0; // {0, 0, 2} 00236 00237 outputValues(5, i0, 0) = 0.0; // {2, 0, 0} 00238 outputValues(5, i0, 1) = -(1.0 + z)/8.0; // {1, 1, 0} 00239 outputValues(5, i0, 2) = (1.0 - y)/8.0; // {1, 0, 1} 00240 outputValues(5, i0, 3) = 0.0; // {0, 2, 0} 00241 outputValues(5, i0, 4) = -(1.0 + x)/8.0; // {0, 1, 1} 00242 outputValues(5, i0, 5) = 0.0; // {0, 0, 2} 00243 00244 outputValues(6, i0, 0) = 0.0; // {2, 0, 0} 00245 outputValues(6, i0, 1) = (1.0 + z)/8.0; // {1, 1, 0} 00246 outputValues(6, i0, 2) = (1.0 + y)/8.0; // {1, 0, 1} 00247 outputValues(6, i0, 3) = 0.0; // {0, 2, 0} 00248 outputValues(6, i0, 4) = (1.0 + x)/8.0; // {0, 1, 1} 00249 outputValues(6, i0, 5) = 0.0; // {0, 0, 2} 00250 00251 outputValues(7, i0, 0) = 0.0; // {2, 0, 0} 00252 outputValues(7, i0, 1) = -(1.0 + z)/8.0; // {1, 1, 0} 00253 outputValues(7, i0, 2) = -(1.0 + y)/8.0; // {1, 0, 1} 00254 outputValues(7, i0, 3) = 0.0; // {0, 2, 0} 00255 outputValues(7, i0, 4) = (1.0 - x)/8.0; // {0, 1, 1} 00256 outputValues(7, i0, 5) = 0.0; // {0, 0, 2} 00257 } 00258 break; 00259 00260 case OPERATOR_D3: 00261 case OPERATOR_D4: 00262 case OPERATOR_D5: 00263 case OPERATOR_D6: 00264 case OPERATOR_D7: 00265 case OPERATOR_D8: 00266 case OPERATOR_D9: 00267 case OPERATOR_D10: 00268 { 00269 // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, DkCardinality) 00270 int DkCardinality = Intrepid::getDkCardinality(operatorType, 00271 this -> basisCellTopology_.getDimension() ); 00272 for(int dofOrd = 0; dofOrd < this -> basisCardinality_; dofOrd++) { 00273 for (int i0 = 0; i0 < dim0; i0++) { 00274 for(int dkOrd = 0; dkOrd < DkCardinality; dkOrd++){ 00275 outputValues(dofOrd, i0, dkOrd) = 0.0; 00276 } 00277 } 00278 } 00279 } 00280 break; 00281 00282 default: 00283 TEUCHOS_TEST_FOR_EXCEPTION( !( Intrepid::isValidOperator(operatorType) ), std::invalid_argument, 00284 ">>> ERROR (Basis_HGRAD_HEX_C1_FEM): Invalid operator type"); 00285 } 00286 } 00287 00288 00289 00290 template<class Scalar, class ArrayScalar> 00291 void Basis_HGRAD_HEX_C1_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar& outputValues, 00292 const ArrayScalar & inputPoints, 00293 const ArrayScalar & cellVertices, 00294 const EOperator operatorType) const { 00295 TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error, 00296 ">>> ERROR (Basis_HGRAD_HEX_C1_FEM): FEM Basis calling an FVD member function"); 00297 } 00298 00299 template<class Scalar, class ArrayScalar> 00300 void Basis_HGRAD_HEX_C1_FEM<Scalar, ArrayScalar>::getDofCoords(ArrayScalar & DofCoords) const { 00301 #ifdef HAVE_INTREPID_DEBUG 00302 // Verify rank of output array. 00303 TEUCHOS_TEST_FOR_EXCEPTION( !(DofCoords.rank() == 2), std::invalid_argument, 00304 ">>> ERROR: (Intrepid::Basis_HGRAD_HEX_C1_FEM::getDofCoords) rank = 2 required for DofCoords array"); 00305 // Verify 0th dimension of output array. 00306 TEUCHOS_TEST_FOR_EXCEPTION( !( DofCoords.dimension(0) == this -> basisCardinality_ ), std::invalid_argument, 00307 ">>> ERROR: (Intrepid::Basis_HGRAD_HEX_C1_FEM::getDofCoords) mismatch in number of DoF and 0th dimension of DofCoords array"); 00308 // Verify 1st dimension of output array. 00309 TEUCHOS_TEST_FOR_EXCEPTION( !( DofCoords.dimension(1) == (int)(this -> basisCellTopology_.getDimension()) ), std::invalid_argument, 00310 ">>> ERROR: (Intrepid::Basis_HGRAD_HEX_C1_FEM::getDofCoords) incorrect reference cell (1st) dimension in DofCoords array"); 00311 #endif 00312 00313 DofCoords(0,0) = -1.0; DofCoords(0,1) = -1.0; DofCoords(0,2) = -1.0; 00314 DofCoords(1,0) = 1.0; DofCoords(1,1) = -1.0; DofCoords(1,2) = -1.0; 00315 DofCoords(2,0) = 1.0; DofCoords(2,1) = 1.0; DofCoords(2,2) = -1.0; 00316 DofCoords(3,0) = -1.0; DofCoords(3,1) = 1.0; DofCoords(3,2) = -1.0; 00317 DofCoords(4,0) = -1.0; DofCoords(4,1) = -1.0; DofCoords(4,2) = 1.0; 00318 DofCoords(5,0) = 1.0; DofCoords(5,1) = -1.0; DofCoords(5,2) = 1.0; 00319 DofCoords(6,0) = 1.0; DofCoords(6,1) = 1.0; DofCoords(6,2) = 1.0; 00320 DofCoords(7,0) = -1.0; DofCoords(7,1) = 1.0; DofCoords(7,2) = 1.0; 00321 } 00322 00323 }// namespace Intrepid 00324 00325 #endif 00326
1.7.6.1