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