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