|
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_QUAD_C2_FEM<Scalar, ArrayScalar>::Basis_HGRAD_QUAD_C2_FEM() 00039 { 00040 this -> basisCardinality_ = 9; 00041 this -> basisDegree_ = 2; 00042 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<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_QUAD_C2_FEM<Scalar, ArrayScalar>::initializeTags() { 00051 00052 // Basis-dependent intializations 00053 int tagSize = 4; // size of DoF tag, i.e., number of fields in the 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 // edge midpoints 00064 1, 0, 0, 1, 00065 1, 1, 0, 1, 00066 1, 2, 0, 1, 00067 1, 3, 0, 1, 00068 // quad center 00069 2, 0, 0, 1}; 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_QUAD_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) coordinates of the evaluation point 00102 Scalar x = 0.0; 00103 Scalar y = 0.0; 00104 00105 switch (operatorType) { 00106 00107 case OPERATOR_VALUE: 00108 for (int i0 = 0; i0 < dim0; i0++) { 00109 x = inputPoints(i0, 0); 00110 y = inputPoints(i0, 1); 00111 00112 // outputValues is a rank-2 array with dimensions (basisCardinality_, dim0) 00113 outputValues(0, i0) = x*(x - 1.0)*y*(y - 1.0)/4.0; 00114 outputValues(1, i0) = x*(x + 1.0)*y*(y - 1.0)/4.0; 00115 outputValues(2, i0) = x*(x + 1.0)*y*(y + 1.0)/4.0; 00116 outputValues(3, i0) = x*(x - 1.0)*y*(y + 1.0)/4.0; 00117 // edge midpoints basis functions 00118 outputValues(4, i0) = (1.0 - x)*(1.0 + x)*y*(y - 1.0)/2.0; 00119 outputValues(5, i0) = x*(x + 1.0)*(1.0 - y)*(1.0 + y)/2.0; 00120 outputValues(6, i0) = (1.0 - x)*(1.0 + x)*y*(y + 1.0)/2.0; 00121 outputValues(7, i0) = x*(x - 1.0)*(1.0 - y)*(1.0 + y)/2.0; 00122 // quad bubble basis function 00123 outputValues(8, i0) = (1.0 - x)*(1.0 + x)*(1.0 - y)*(1.0 + y); 00124 } 00125 break; 00126 00127 case OPERATOR_GRAD: 00128 case OPERATOR_D1: 00129 for (int i0 = 0; i0 < dim0; i0++) { 00130 x = inputPoints(i0,0); 00131 y = inputPoints(i0,1); 00132 00133 // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim) 00134 outputValues(0, i0, 0) = (-0.25 + 0.5*x)*(-1. + y)*y; 00135 outputValues(0, i0, 1) = (-1.0 + x)*x*(-0.25 + 0.5*y); 00136 00137 outputValues(1, i0, 0) = (0.25 + 0.5*x)*(-1. + y)*y; 00138 outputValues(1, i0, 1) = x*(1. + x)*(-0.25 + 0.5*y); 00139 00140 outputValues(2, i0, 0) = (0.25 + 0.5*x)*y*(1. + y); 00141 outputValues(2, i0, 1) = x*(1. + x)*(0.25 + 0.5*y); 00142 00143 outputValues(3, i0, 0) = (-0.25 + 0.5*x)*y*(1. + y); 00144 outputValues(3, i0, 1) = (-1. + x)*x*(0.25 + 0.5*y); 00145 00146 outputValues(4, i0, 0) = x*(1.0 - y)*y; 00147 outputValues(4, i0, 1) = 0.5*(1.0 - x)*(1.0 + x)*(-1.0 + 2.0*y); 00148 00149 outputValues(5, i0, 0) = 0.5*(1.0 - y)*(1.0 + y)*(1.0 + 2.0*x); 00150 outputValues(5, i0, 1) =-x*(1.0 + x)*y; 00151 00152 outputValues(6, i0, 0) =-y*(1.0 + y)*x; 00153 outputValues(6, i0, 1) = 0.5*(1.0 - x)*(1.0 + x)*(1.0 + 2.0*y); 00154 00155 outputValues(7, i0, 0) = 0.5*(1.0 - y)*(1.0+ y)*(-1.0 + 2.0*x); 00156 outputValues(7, i0, 1) = (1.0 - x)*x*y; 00157 00158 outputValues(8, i0, 0) =-2.0*(1.0 - y)*(1.0 + y)*x; 00159 outputValues(8, i0, 1) =-2.0*(1.0 - x)*(1.0 + x)*y; 00160 } 00161 break; 00162 00163 case OPERATOR_CURL: 00164 for (int i0 = 0; i0 < dim0; i0++) { 00165 x = inputPoints(i0,0); 00166 y = inputPoints(i0,1); 00167 00168 // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim) 00169 // CURL(u) = (u_y, -u_x), is rotated GRAD 00170 outputValues(0, i0, 1) =-(-0.25 + 0.5*x)*(-1. + y)*y; 00171 outputValues(0, i0, 0) = (-1.0 + x)*x*(-0.25 + 0.5*y); 00172 00173 outputValues(1, i0, 1) =-(0.25 + 0.5*x)*(-1. + y)*y; 00174 outputValues(1, i0, 0) = x*(1. + x)*(-0.25 + 0.5*y); 00175 00176 outputValues(2, i0, 1) =-(0.25 + 0.5*x)*y*(1. + y); 00177 outputValues(2, i0, 0) = x*(1. + x)*(0.25 + 0.5*y); 00178 00179 outputValues(3, i0, 1) =-(-0.25 + 0.5*x)*y*(1. + y); 00180 outputValues(3, i0, 0) = (-1. + x)*x*(0.25 + 0.5*y); 00181 00182 outputValues(4, i0, 1) =-x*(1.0 - y)*y; 00183 outputValues(4, i0, 0) = 0.5*(1.0 - x)*(1.0 + x)*(-1.0 + 2.0*y); 00184 00185 outputValues(5, i0, 1) =-0.5*(1.0 - y)*(1.0 + y)*(1.0 + 2.0*x); 00186 outputValues(5, i0, 0) =-x*(1.0 + x)*y; 00187 00188 outputValues(6, i0, 1) = y*(1.0 + y)*x; 00189 outputValues(6, i0, 0) = 0.5*(1.0 - x)*(1.0 + x)*(1.0 + 2.0*y); 00190 00191 outputValues(7, i0, 1) =-0.5*(1.0 - y)*(1.0 + y)*(-1.0 + 2.0*x); 00192 outputValues(7, i0, 0) = (1.0 - x)*x*y; 00193 00194 outputValues(8, i0, 1) = 2.0*(1.0 - y)*(1.0 + y)*x; 00195 outputValues(8, i0, 0) =-2.0*(1.0 - x)*(1.0 + x)*y; 00196 } 00197 break; 00198 00199 case OPERATOR_DIV: 00200 TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument, 00201 ">>> ERROR (Basis_HGRAD_QUAD_C2_FEM): DIV is invalid operator for rank-0 (scalar) functions in 2D"); 00202 break; 00203 00204 case OPERATOR_D2: 00205 for (int i0 = 0; i0 < dim0; i0++) { 00206 x = inputPoints(i0,0); 00207 y = inputPoints(i0,1); 00208 00209 // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, D2Cardinality=3) 00210 outputValues(0, i0, 0) = 0.5*(-1.0 + y)*y; 00211 outputValues(0, i0, 1) = 0.25 - 0.5*y + x*(-0.5 + 1.*y); 00212 outputValues(0, i0, 2) = 0.5*(-1.0 + x)*x; 00213 00214 outputValues(1, i0, 0) = 0.5*(-1.0 + y)*y; 00215 outputValues(1, i0, 1) =-0.25 + 0.5*y + x*(-0.5 + 1.*y); 00216 outputValues(1, i0, 2) = 0.5*x*(1.0 + x); 00217 00218 outputValues(2, i0, 0) = 0.5*y*(1.0 + y); 00219 outputValues(2, i0, 1) = 0.25 + 0.5*y + x*(0.5 + 1.*y); 00220 outputValues(2, i0, 2) = 0.5*x*(1.0 + x); 00221 00222 outputValues(3, i0, 0) = 0.5*y*(1.0 + y); 00223 outputValues(3, i0, 1) =-0.25 - 0.5*y + x*(0.5 + 1.*y); 00224 outputValues(3, i0, 2) = 0.5*(-1.0 + x)*x; 00225 00226 outputValues(4, i0, 0) = (1.0 - y)*y; 00227 outputValues(4, i0, 1) = x*(1. - 2.*y); 00228 outputValues(4, i0, 2) = (1.0 - x)*(1.0 + x); 00229 00230 outputValues(5, i0, 0) = (1.0 - y)*(1.0 + y); 00231 outputValues(5, i0, 1) = x*(0. - 2.*y) - 1.*y; 00232 outputValues(5, i0, 2) =-x*(1.0 + x); 00233 00234 outputValues(6, i0, 0) =-y*(1.0 + y); 00235 outputValues(6, i0, 1) = x*(-1. - 2.*y); 00236 outputValues(6, i0, 2) = (1.0 - x)*(1.0 + x); 00237 00238 outputValues(7, i0, 0) = (1.0 - y)*(1.0 + y); 00239 outputValues(7, i0, 1) = x*(0. - 2.*y) + 1.*y; 00240 outputValues(7, i0, 2) = (1.0 - x)*x; 00241 00242 outputValues(8, i0, 0) =-2.0 + 2.0*y*y; 00243 outputValues(8, i0, 1) = 4*x*y; 00244 outputValues(8, i0, 2) =-2.0 + 2.0*x*x; 00245 00246 } 00247 break; 00248 00249 case OPERATOR_D3: 00250 // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, D3Cardinality=4) 00251 for (int i0 = 0; i0 < dim0; i0++) { 00252 x = inputPoints(i0,0); 00253 y = inputPoints(i0,1); 00254 00255 outputValues(0, i0, 0) = 0.0; 00256 outputValues(0, i0, 1) =-0.5 + y; 00257 outputValues(0, i0, 2) =-0.5 + x; 00258 outputValues(0, i0, 3) = 0.0; 00259 00260 outputValues(1, i0, 0) = 0.0; 00261 outputValues(1, i0, 1) =-0.5 + y; 00262 outputValues(1, i0, 2) = 0.5 + x; 00263 outputValues(1, i0, 3) = 0.0; 00264 00265 outputValues(2, i0, 0) = 0.0; 00266 outputValues(2, i0, 1) = 0.5 + y; 00267 outputValues(2, i0, 2) = 0.5 + x; 00268 outputValues(2, i0, 3) = 0.0; 00269 00270 outputValues(3, i0, 0) = 0.0; 00271 outputValues(3, i0, 1) = 0.5 + y; 00272 outputValues(3, i0, 2) =-0.5 + x; 00273 outputValues(3, i0, 3) = 0.0; 00274 00275 outputValues(4, i0, 0) = 0.0; 00276 outputValues(4, i0, 1) = 1.0 - 2.0*y; 00277 outputValues(4, i0, 2) =-2.0*x; 00278 outputValues(4, i0, 3) = 0.0; 00279 00280 outputValues(5, i0, 0) = 0.0; 00281 outputValues(5, i0, 1) =-2.0*y; 00282 outputValues(5, i0, 2) =-1.0 - 2.0*x; 00283 outputValues(5, i0, 3) = 0.0; 00284 00285 outputValues(6, i0, 0) = 0.0; 00286 outputValues(6, i0, 1) =-1.0 - 2.0*y; 00287 outputValues(6, i0, 2) =-2.0*x; 00288 outputValues(6, i0, 3) = 0.0; 00289 00290 outputValues(7, i0, 0) = 0.0; 00291 outputValues(7, i0, 1) =-2.0*y; 00292 outputValues(7, i0, 2) = 1.0 - 2.0*x; 00293 outputValues(7, i0, 3) = 0.0; 00294 00295 outputValues(8, i0, 0) = 0.0; 00296 outputValues(8, i0, 1) = 4.0*y; 00297 outputValues(8, i0, 2) = 4.0*x; 00298 outputValues(8, i0, 3) = 0.0; 00299 } 00300 break; 00301 00302 case OPERATOR_D4: 00303 // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, D4Cardinality=5) 00304 for (int i0 = 0; i0 < dim0; i0++) { 00305 00306 outputValues(0, i0, 0) = 0.0; 00307 outputValues(0, i0, 1) = 0.0; 00308 outputValues(0, i0, 2) = 1.0; 00309 outputValues(0, i0, 3) = 0.0; 00310 outputValues(0, i0, 4) = 0.0; 00311 00312 outputValues(1, i0, 0) = 0.0; 00313 outputValues(1, i0, 1) = 0.0; 00314 outputValues(1, i0, 2) = 1.0; 00315 outputValues(1, i0, 3) = 0.0; 00316 outputValues(1, i0, 4) = 0.0; 00317 00318 outputValues(2, i0, 0) = 0.0; 00319 outputValues(2, i0, 1) = 0.0; 00320 outputValues(2, i0, 2) = 1.0; 00321 outputValues(2, i0, 3) = 0.0; 00322 outputValues(2, i0, 4) = 0.0; 00323 00324 outputValues(3, i0, 0) = 0.0; 00325 outputValues(3, i0, 1) = 0.0; 00326 outputValues(3, i0, 2) = 1.0; 00327 outputValues(3, i0, 3) = 0.0; 00328 outputValues(3, i0, 4) = 0.0; 00329 00330 outputValues(4, i0, 0) = 0.0; 00331 outputValues(4, i0, 1) = 0.0; 00332 outputValues(4, i0, 2) =-2.0; 00333 outputValues(4, i0, 3) = 0.0; 00334 outputValues(4, i0, 4) = 0.0; 00335 00336 outputValues(5, i0, 0) = 0.0; 00337 outputValues(5, i0, 1) = 0.0; 00338 outputValues(5, i0, 2) =-2.0; 00339 outputValues(5, i0, 3) = 0.0; 00340 outputValues(5, i0, 4) = 0.0; 00341 00342 outputValues(6, i0, 0) = 0.0; 00343 outputValues(6, i0, 1) = 0.0; 00344 outputValues(6, i0, 2) =-2.0; 00345 outputValues(6, i0, 3) = 0.0; 00346 outputValues(6, i0, 4) = 0.0; 00347 00348 outputValues(7, i0, 0) = 0.0; 00349 outputValues(7, i0, 1) = 0.0; 00350 outputValues(7, i0, 2) =-2.0; 00351 outputValues(7, i0, 3) = 0.0; 00352 outputValues(7, i0, 4) = 0.0; 00353 00354 outputValues(8, i0, 0) = 0.0; 00355 outputValues(8, i0, 1) = 0.0; 00356 outputValues(8, i0, 2) = 4.0; 00357 outputValues(8, i0, 3) = 0.0; 00358 outputValues(8, i0, 4) = 0.0; 00359 } 00360 break; 00361 00362 case OPERATOR_D5: 00363 case OPERATOR_D6: 00364 case OPERATOR_D7: 00365 case OPERATOR_D8: 00366 case OPERATOR_D9: 00367 case OPERATOR_D10: 00368 { 00369 // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, DkCardinality) 00370 int DkCardinality = Intrepid::getDkCardinality(operatorType, 00371 this -> basisCellTopology_.getDimension() ); 00372 for(int dofOrd = 0; dofOrd < this -> basisCardinality_; dofOrd++) { 00373 for (int i0 = 0; i0 < dim0; i0++) { 00374 for(int dkOrd = 0; dkOrd < DkCardinality; dkOrd++){ 00375 outputValues(dofOrd, i0, dkOrd) = 0.0; 00376 } 00377 } 00378 } 00379 } 00380 break; 00381 00382 default: 00383 TEST_FOR_EXCEPTION( !( Intrepid::isValidOperator(operatorType) ), std::invalid_argument, 00384 ">>> ERROR (Basis_HGRAD_QUAD_C2_FEM): Invalid operator type"); 00385 } 00386 } 00387 00388 00389 00390 template<class Scalar, class ArrayScalar> 00391 void Basis_HGRAD_QUAD_C2_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar& outputValues, 00392 const ArrayScalar & inputPoints, 00393 const ArrayScalar & cellVertices, 00394 const EOperator operatorType) const { 00395 TEST_FOR_EXCEPTION( (true), std::logic_error, 00396 ">>> ERROR (Basis_HGRAD_QUAD_C2_FEM): FEM Basis calling an FVD member function"); 00397 } 00398 00399 00400 00401 template<class Scalar, class ArrayScalar> 00402 void Basis_HGRAD_QUAD_C2_FEM<Scalar, ArrayScalar>::getDofCoords(ArrayScalar & DofCoords) const { 00403 #ifdef HAVE_INTREPID_DEBUG 00404 // Verify rank of output array. 00405 TEST_FOR_EXCEPTION( !(DofCoords.rank() == 2), std::invalid_argument, 00406 ">>> ERROR: (Intrepid::Basis_HGRAD_QUAD_C2_FEM::getDofCoords) rank = 2 required for DofCoords array"); 00407 // Verify 0th dimension of output array. 00408 TEST_FOR_EXCEPTION( !( DofCoords.dimension(0) == this -> basisCardinality_ ), std::invalid_argument, 00409 ">>> ERROR: (Intrepid::Basis_HGRAD_QUAD_C2_FEM::getDofCoords) mismatch in number of DoF and 0th dimension of DofCoords array"); 00410 // Verify 1st dimension of output array. 00411 TEST_FOR_EXCEPTION( !( DofCoords.dimension(1) == (int)(this -> basisCellTopology_.getDimension()) ), std::invalid_argument, 00412 ">>> ERROR: (Intrepid::Basis_HGRAD_QUAD_C2_FEM::getDofCoords) incorrect reference cell (1st) dimension in DofCoords array"); 00413 #endif 00414 00415 DofCoords(0,0) = -1.0; DofCoords(0,1) = -1.0; 00416 DofCoords(1,0) = 1.0; DofCoords(1,1) = -1.0; 00417 DofCoords(2,0) = 1.0; DofCoords(2,1) = 1.0; 00418 DofCoords(3,0) = -1.0; DofCoords(3,1) = 1.0; 00419 00420 DofCoords(4,0) = 0.0; DofCoords(4,1) = -1.0; 00421 DofCoords(5,0) = 1.0; DofCoords(5,1) = 0.0; 00422 DofCoords(6,0) = 0.0; DofCoords(6,1) = 1.0; 00423 DofCoords(7,0) = -1.0; DofCoords(7,1) = 0.0; 00424 00425 DofCoords(8,0) = 0.0; DofCoords(8,1) = 0.0; 00426 00427 } 00428 00429 }// namespace Intrepid
1.7.4