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