|
Intrepid
|
00001 #ifndef INTREPID_HGRAD_PYR_I2_FEMDEF_HPP 00002 #define INTREPID_HGRAD_PYR_I2_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_PYR_I2_FEM<Scalar, ArrayScalar>::Basis_HGRAD_PYR_I2_FEM() 00055 { 00056 this -> basisCardinality_ = 13; 00057 this -> basisDegree_ = 2; 00058 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Pyramid<5> >() ); 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_PYR_I2_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 0, 4, 0, 1, 00080 1, 0, 0, 1, 00081 1, 1, 0, 1, 00082 1, 2, 0, 1, 00083 1, 3, 0, 1, 00084 1, 4, 0, 1, 00085 1, 5, 0, 1, 00086 1, 6, 0, 1, 00087 1, 7, 0, 1, 00088 }; 00089 00090 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays: 00091 Intrepid::setOrdinalTagData(this -> tagToOrdinal_, 00092 this -> ordinalToTag_, 00093 tags, 00094 this -> basisCardinality_, 00095 tagSize, 00096 posScDim, 00097 posScOrd, 00098 posDfOrd); 00099 } 00100 00101 00102 00103 template<class Scalar, class ArrayScalar> 00104 void Basis_HGRAD_PYR_I2_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar & outputValues, 00105 const ArrayScalar & inputPoints, 00106 const EOperator operatorType) const { 00107 00108 // Verify arguments 00109 #ifdef HAVE_INTREPID_DEBUG 00110 Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues, 00111 inputPoints, 00112 operatorType, 00113 this -> getBaseCellTopology(), 00114 this -> getCardinality() ); 00115 #endif 00116 00117 // Number of evaluation points = dim 0 of inputPoints 00118 int dim0 = inputPoints.dimension(0); 00119 00120 // Temporaries: (x,y,z) coordinates of the evaluation point 00121 Scalar x = 0.0; 00122 Scalar y = 0.0; 00123 Scalar z = 0.0; 00124 const Scalar eps = std::numeric_limits<Scalar>::epsilon( ); 00125 00126 switch (operatorType) { 00127 00128 case OPERATOR_VALUE: 00129 for (int i0 = 0; i0 < dim0; i0++) { 00130 x = inputPoints(i0, 0); 00131 y = inputPoints(i0, 1); 00132 z = inputPoints(i0, 2); 00133 00134 //be sure that the basis functions are defined when z is very close to 1. 00135 00136 if(fabs(z-1.0) < eps) { 00137 if(z <= 1.0) z = 1.0-eps; 00138 else z = 1.0+eps; 00139 //z = 1.0-eps; 00140 } 00141 00142 Scalar w = 1.0/(1.0 - z); 00143 00144 // outputValues is a rank-2 array with dimensions (basisCardinality_, dim0) 00145 outputValues(0, i0) = 0.25 * (-x - y - 1.0)*((1.0-x)*(1.0-y) - z + x*y*z*w); 00146 outputValues(1, i0) = 0.25 * ( x - y - 1.0)*((1.0+x)*(1.0-y) - z - x*y*z*w); 00147 outputValues(2, i0) = 0.25 * ( x + y - 1.0)*((1.0+x)*(1.0+y) - z + x*y*z*w); 00148 outputValues(3, i0) = 0.25 * (-x + y - 1.0)*((1.0-x)*(1.0+y) - z - x*y*z*w); 00149 00150 outputValues(4, i0) = z * (2.0*z - 1.0); 00151 00152 outputValues(5, i0) = 0.5 * (1.0 + x - z)*(1.0 - x - z)*(1.0 - y - z)*w; 00153 outputValues(6, i0) = 0.5 * (1.0 + y - z)*(1.0 - y - z)*(1.0 + x - z)*w; 00154 outputValues(7, i0) = 0.5 * (1.0 + x - z)*(1.0 - x - z)*(1.0 + y - z)*w; 00155 outputValues(8, i0) = 0.5 * (1.0 + y - z)*(1.0 - y - z)*(1.0 - x - z)*w; 00156 00157 outputValues(9, i0) = z*(1.0 - x - z)*(1.0 - y - z)*w; 00158 outputValues(10,i0) = z*(1.0 + x - z)*(1.0 - y - z)*w; 00159 outputValues(11,i0) = z*(1.0 + x - z)*(1.0 + y - z)*w; 00160 outputValues(12,i0) = z*(1.0 - x - z)*(1.0 + y - z)*w; 00161 } 00162 break; 00163 00164 case OPERATOR_GRAD: 00165 case OPERATOR_D1: 00166 for (int i0 = 0; i0 < dim0; i0++) { 00167 x = inputPoints(i0,0); 00168 y = inputPoints(i0,1); 00169 z = inputPoints(i0,2); 00170 00171 //be sure that the basis functions are defined when z is very close to 1. 00172 00173 if(fabs(z-1.0) < eps) { 00174 if(z <= 1.0) z = 1.0-eps; 00175 else z = 1.0+eps; 00176 } 00177 00178 Scalar w = 1.0/(1.0 - z); 00179 00180 // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim) 00181 outputValues(0, i0, 0) = 0.25*(-1.0-x-y)*(-1.0+y + y*z*w) - 0.25*((1.0-x)*(1.0-y)-z + x*y*z*w); 00182 outputValues(0, i0, 1) = 0.25*(-1.0-x-y)*(-1.0+x + x*z*w) - 0.25*((1.0-x)*(1.0-y)-z + x*y*z*w); 00183 outputValues(0, i0, 2) = 0.25*(-1.0-x-y)*(-1.0 + x*y*w + x*y*z*w*w); 00184 00185 outputValues(1, i0, 0) = 0.25*(-1.0+x-y)*( 1.0-y - y*z*w) + 0.25*((1.0+x)*(1.0-y)-z - x*y*z*w); 00186 outputValues(1, i0, 1) = 0.25*(-1.0+x-y)*(-1.0-x - x*z*w) - 0.25*((1.0+x)*(1.0-y)-z - x*y*z*w); 00187 outputValues(1, i0, 2) = 0.25*(-1.0+x-y)*(-1.0 - x*y*w - x*y*z*w*w); 00188 00189 outputValues(2, i0, 0) = 0.25*(-1.0+x+y)*(1.0+y + y*z*w) + 0.25*((1.0+x)*(1.0+y)-z + x*y*z*w); 00190 outputValues(2, i0, 1) = 0.25*(-1.0+x+y)*(1.0+x + x*z*w) + 0.25*((1.0+x)*(1.0+y)-z + x*y*z*w); 00191 outputValues(2, i0, 2) = 0.25*(-1.0+x+y)*(-1.0 + x*y*w + x*y*z*w*w); 00192 00193 outputValues(3, i0, 0) = 0.25*(-1.0-x+y)*(-1.0-y - y*z*w) - 0.25*((1.0-x)*(1.0+y)-z - x*y*z*w); 00194 outputValues(3, i0, 1) = 0.25*(-1.0-x+y)*( 1.0-x - x*z*w) + 0.25*((1.0-x)*(1.0+y)-z - x*y*z*w); 00195 outputValues(3, i0, 2) = 0.25*(-1.0-x+y)*(-1.0 - x*y*w - x*y*z*w*w); 00196 00197 outputValues(4, i0, 0) = 0.0; 00198 outputValues(4, i0, 1) = 0.0; 00199 outputValues(4, i0, 2) = -1.0 + 4.0*z; 00200 00201 outputValues(5, i0, 0) = -x*w*(1.0-y-z); 00202 outputValues(5, i0, 1) = -0.5*(1.0-x-z)*(1.0+x-z)*w; 00203 outputValues(5, i0, 2) = 0.5*y*x*x*w*w + 0.5*y - 1.0+z; 00204 00205 outputValues(6, i0, 0) = 0.5*(1.0-y-z)*(1.0+y-z)*w; 00206 outputValues(6, i0, 1) = -y*w*(1.0+x-z); 00207 outputValues(6, i0, 2) = -0.5*x*y*y*w*w - 0.5*x - 1.0+z; 00208 00209 outputValues(7, i0, 0) = -x*w*(1.0+y-z); 00210 outputValues(7, i0, 1) = 0.5*(1.0-x-z)*(1.0+x-z)*w; 00211 outputValues(7, i0, 2) = -0.5*y*x*x*w*w - 0.5*y - 1.0+z; 00212 00213 outputValues(8, i0, 0) = -0.5*(1.0-y-z)*(1.0+y-z)*w; 00214 outputValues(8, i0, 1) = -y*w*(1.0-x-z); 00215 outputValues(8, i0, 2) = 0.5*x*y*y*w*w + 0.5*x - 1.0+z; 00216 00217 outputValues(9, i0, 0) = -(1.0-y-z)*z*w; 00218 outputValues(9, i0, 1) = -(1.0-x-z)*z*w; 00219 outputValues(9, i0, 2) = x*y*w*w + 1.0 - x - y - 2.0*z; 00220 00221 outputValues(10,i0, 0) = (1.0-y-z)*z*w; 00222 outputValues(10,i0, 1) = -(1.0+x-z)*z*w; 00223 outputValues(10,i0, 2) = -x*y*w*w + 1.0 + x - y - 2.0*z; 00224 00225 outputValues(11,i0, 0) = (1.0+y-z)*z*w; 00226 outputValues(11,i0, 1) = (1.0+x-z)*z*w; 00227 outputValues(11,i0, 2) = x*y*w*w + 1.0 + x + y - 2.0*z; 00228 00229 outputValues(12,i0, 0) = -(1.0+y-z)*z*w; 00230 outputValues(12,i0, 1) = (1.0-x-z)*z*w; 00231 outputValues(12,i0, 2) = -x*y*w*w + 1.0 - x + y - 2.0*z; 00232 00233 } 00234 break; 00235 00236 case OPERATOR_CURL: 00237 TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_CURL), std::invalid_argument, 00238 ">>> ERROR (Basis_HGRAD_PYR_I2_FEM): CURL is invalid operator for rank-0 (scalar) functions in 3D"); 00239 break; 00240 00241 case OPERATOR_DIV: 00242 TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument, 00243 ">>> ERROR (Basis_HGRAD_PYR_I2_FEM): DIV is invalid operator for rank-0 (scalar) functions in 3D"); 00244 break; 00245 00246 case OPERATOR_D2: 00247 for (int i0 = 0; i0 < dim0; i0++) { 00248 x = inputPoints(i0,0); 00249 y = inputPoints(i0,1); 00250 z = inputPoints(i0,2); 00251 00252 if(fabs(z-1.0) < eps) { 00253 if(z <= 1.0) z = 1.0-eps; 00254 else z = 1.0+eps; 00255 } 00256 Scalar w = 1.0/(1.0 - z); 00257 00258 outputValues(0, i0, 0) = -0.5*(-1.0+y+y*z*w); 00259 outputValues(0, i0, 1) = -(-0.25 + 0.5*x + 0.5*y + 0.5*z)*w; 00260 outputValues(0, i0, 2) = 0.25 + (-0.25*y-0.5*x*y-0.25*y*y)*w*w; 00261 outputValues(0, i0, 3) = -0.5*(-1.0+x+x*z*w); 00262 outputValues(0, i0, 4) = 0.25 + (-0.25*x*x-0.25*x-0.5*x*y)*w*w; 00263 outputValues(0, i0, 5) = 0.5*x*y*(-1.0-x-y)*w*w*w; 00264 00265 outputValues(1, i0, 0) = 0.5*(1.0-y-y*z*w); 00266 outputValues(1, i0, 1) =-(0.25 + 0.5*x - 0.5*y - 0.5*z)*w; 00267 outputValues(1, i0, 2) = -0.25 + (0.25*y-0.5*x*y+0.25*y*y)*w*w; 00268 outputValues(1, i0, 3) = -0.5*(-1.0-x-x*z*w); 00269 outputValues(1, i0, 4) = 0.25 + (-0.25*x*x + 0.25*x + 0.5*x*y)*w*w; 00270 outputValues(1, i0, 5) = -0.5*x*y*(-1.0+x-y)*w*w*w; 00271 00272 outputValues(2, i0, 0) = 0.5*(1.0+y+y*z*w); 00273 outputValues(2, i0, 1) =-(-0.25 - 0.5*x - 0.5*y + 0.5*z)*w; 00274 outputValues(2, i0, 2) = -0.25 + (-0.25*y+0.5*x*y+0.25*y*y)*w*w; 00275 outputValues(2, i0, 3) = 0.5*(1.0+x+x*z*w); 00276 outputValues(2, i0, 4) = -0.25 + (0.25*x*x -0.25*x + 0.5*x*y)*w*w; 00277 outputValues(2, i0, 5) = 0.5*x*y*(-1.0+x+y)*w*w*w; 00278 00279 outputValues(3, i0, 0) = -0.5*(-1.0-y-y*z*w); 00280 outputValues(3, i0, 1) =-(0.25 - 0.5*x + 0.5*y - 0.5*z)*w; 00281 outputValues(3, i0, 2) = 0.25 + (0.25*y+0.5*x*y-0.25*y*y)*w*w; 00282 outputValues(3, i0, 3) = 0.5*(1.0-x-x*z*w); 00283 outputValues(3, i0, 4) = -0.25 + (0.25*x + 0.25*x*x - 0.5*x*y)*w*w; 00284 outputValues(3, i0, 5) = -0.5*x*y*(-1.0-x+y)*w*w*w; 00285 00286 outputValues(4, i0, 0) = 0.0; 00287 outputValues(4, i0, 1) = 0.0; 00288 outputValues(4, i0, 2) = 0.0; 00289 outputValues(4, i0, 3) = 0.0; 00290 outputValues(4, i0, 4) = 0.0; 00291 outputValues(4, i0, 5) = 4.0; 00292 00293 outputValues(5, i0, 0) = -(1.0-y-z)*w; 00294 outputValues(5, i0, 1) = x*w; 00295 outputValues(5, i0, 2) = x*y*w*w; 00296 outputValues(5, i0, 3) = 0.0; 00297 outputValues(5, i0, 4) = 0.5*x*x*w*w + 0.5; 00298 outputValues(5, i0, 5) = x*x*y*w*w*w + 1.0; 00299 00300 outputValues(6, i0, 0) = 0.0; 00301 outputValues(6, i0, 1) = -y*w; 00302 outputValues(6, i0, 2) = -0.5*y*y*w*w - 0.5; 00303 outputValues(6, i0, 3) =-(1.0+x-z)*w; 00304 outputValues(6, i0, 4) = -x*y*w*w; 00305 outputValues(6, i0, 5) = -x*y*y*w*w*w + 1.0; 00306 00307 outputValues(7, i0, 0) = -(1.0+y-z)*w; 00308 outputValues(7, i0, 1) = -x*w; 00309 outputValues(7, i0, 2) = -x*y*w*w; 00310 outputValues(7, i0, 3) = 0.0; 00311 outputValues(7, i0, 4) = -0.5*x*x*w*w - 0.5; 00312 outputValues(7, i0, 5) = -x*x*y*w*w*w + 1.0; 00313 00314 outputValues(8, i0, 0) = 0.0; 00315 outputValues(8, i0, 1) = y*w; 00316 outputValues(8, i0, 2) = 0.5*y*y*w*w + 0.5; 00317 outputValues(8, i0, 3) = -(1.0-x-z)*w; 00318 outputValues(8, i0, 4) = x*y*w*w; 00319 outputValues(8, i0, 5) = x*y*y*w*w*w + 1.0; 00320 00321 outputValues(9, i0, 0) = 0.0; 00322 outputValues(9, i0, 1) = z*w; 00323 outputValues(9, i0, 2) = y*w*w - 1.0; 00324 outputValues(9, i0, 3) = 0.0; 00325 outputValues(9, i0, 4) = x*w*w - 1.0; 00326 outputValues(9, i0, 5) = 2.0*x*y*w*w*w - 2.0; 00327 00328 outputValues(10,i0, 0) = 0.0; 00329 outputValues(10,i0, 1) = -z*w; 00330 outputValues(10,i0, 2) = -y*w*w + 1.0; 00331 outputValues(10,i0, 3) = 0.0; 00332 outputValues(10,i0, 4) = -x*w*w - 1.0; 00333 outputValues(10,i0, 5) = -2.0*x*y*w*w*w - 2.0; 00334 00335 outputValues(11,i0, 0) = 0.0; 00336 outputValues(11,i0, 1) = z*w; 00337 outputValues(11,i0, 2) = y*w*w + 1.0; 00338 outputValues(11,i0, 3) = 0.0; 00339 outputValues(11,i0, 4) = x*w*w + 1.0; 00340 outputValues(11,i0, 5) = 2.0*x*y*w*w*w - 2.0; 00341 00342 outputValues(12,i0, 0) = 0.0; 00343 outputValues(12,i0, 1) = -z*w; 00344 outputValues(12,i0, 2) = -y*w*w - 1.0; 00345 outputValues(12,i0, 3) = 0.0; 00346 outputValues(12,i0, 4) = -x*w*w + 1.0; 00347 outputValues(12,i0, 5) = -2.0*x*y*w*w*w - 2.0; 00348 00349 } 00350 break; 00351 00352 case OPERATOR_D3: 00353 for (int i0 = 0; i0 < dim0; i0++) { 00354 x = inputPoints(i0,0); 00355 y = inputPoints(i0,1); 00356 z = inputPoints(i0,2); 00357 00358 if(fabs(z-1.0) < eps) { 00359 if(z <= 1.0) z = 1.0-eps; 00360 else z = 1.0+eps; 00361 } 00362 00363 Scalar w = 1.0/(1.0 - z); 00364 00365 outputValues(0, i0, 0) = 0.0; 00366 outputValues(0, i0, 1) = -0.5*w; 00367 outputValues(0, i0, 2) = -0.5*y*w*w; 00368 outputValues(0, i0, 3) = -0.5*w; 00369 outputValues(0, i0, 4) =(-0.25 - 0.5*x - 0.5*y)*w*w; 00370 outputValues(0, i0, 5) =-(0.5*y + x*y + 0.5*y*y)*w*w*w; 00371 outputValues(0, i0, 6) = 0.0; 00372 outputValues(0, i0, 7) = -0.5*x*w*w; 00373 outputValues(0, i0, 8) =-(0.5*x + x*y + 0.5*x*x)*w*w*w; 00374 outputValues(0, i0, 9) = 1.5*x*y*(-1.0 - x - y)*w*w*w*w; 00375 00376 outputValues(1, i0, 0) = 0.0; 00377 outputValues(1, i0, 1) = -0.5*w; 00378 outputValues(1, i0, 2) = -0.5*y*w*w; 00379 outputValues(1, i0, 3) = 0.5*w; 00380 outputValues(1, i0, 4) = (0.25 - 0.5*x - 0.5*y)*w*w; 00381 outputValues(1, i0, 5) =-(-0.5*y + x*y - 0.5*y*y)*w*w*w; 00382 outputValues(1, i0, 6) = 0.0; 00383 outputValues(1, i0, 7) = 0.5*x*w*w; 00384 outputValues(1, i0, 8) =-(-0.5*x - x*y + 0.5*x*x)*w*w*w; 00385 outputValues(1, i0, 9) = -1.5*x*y*(-1.0 + x - y)*w*w*w*w; 00386 00387 outputValues(2, i0, 0) = 0.0; 00388 outputValues(2, i0, 1) = 0.5*w; 00389 outputValues(2, i0, 2) = 0.5*y*w*w; 00390 outputValues(2, i0, 3) = 0.5*w; 00391 outputValues(2, i0, 4) = (-0.25 + 0.5*x + 0.5*y)*w*w; 00392 outputValues(2, i0, 5) =-(0.5*y - x*y - 0.5*y*y)*w*w*w; 00393 outputValues(2, i0, 6) = 0.0; 00394 outputValues(2, i0, 7) = 0.5*x*w*w; 00395 outputValues(2, i0, 8) =-(0.5*x - x*y - 0.5*x*x)*w*w*w; 00396 outputValues(2, i0, 9) = 1.5*x*y*(-1.0 + x + y)*w*w*w*w; 00397 00398 outputValues(3, i0, 0) = 0.0; 00399 outputValues(3, i0, 1) = 0.5*w; 00400 outputValues(3, i0, 2) = 0.5*y*w*w; 00401 outputValues(3, i0, 3) = -0.5*w; 00402 outputValues(3, i0, 4) = (0.25 + 0.5*x - 0.5*y)*w*w; 00403 outputValues(3, i0, 5) =-(-0.5*y - x*y + 0.5*y*y)*w*w*w; 00404 outputValues(3, i0, 6) = 0.0; 00405 outputValues(3, i0, 7) = -0.5*x*w*w; 00406 outputValues(3, i0, 8) =-(-0.5*x + x*y - 0.5*x*x)*w*w*w; 00407 outputValues(3, i0, 9) = -1.5*x*y*(-1.0 - x + y)*w*w*w*w; 00408 00409 outputValues(4, i0, 0) = 0.0; 00410 outputValues(4, i0, 1) = 0.0; 00411 outputValues(4, i0, 2) = 0.0; 00412 outputValues(4, i0, 3) = 0.0; 00413 outputValues(4, i0, 4) = 0.0; 00414 outputValues(4, i0, 5) = 0.0; 00415 outputValues(4, i0, 6) = 0.0; 00416 outputValues(4, i0, 7) = 0.0; 00417 outputValues(4, i0, 8) = 0.0; 00418 outputValues(4, i0, 9) = 0.0; 00419 00420 outputValues(5, i0, 0) = 0.0; 00421 outputValues(5, i0, 1) = w; 00422 outputValues(5, i0, 2) = y*w*w; 00423 outputValues(5, i0, 3) = 0.0; 00424 outputValues(5, i0, 4) = x*w*w; 00425 outputValues(5, i0, 5) = 2.0*y*x*w*w*w; 00426 outputValues(5, i0, 6) = 0.0; 00427 outputValues(5, i0, 7) = 0.0; 00428 outputValues(5, i0, 8) = x*x*w*w*w; 00429 outputValues(5, i0, 9) = 3.0*x*x*y*w*w*w*w; 00430 00431 outputValues(6, i0, 0) = 0.0; 00432 outputValues(6, i0, 1) = 0.0; 00433 outputValues(6, i0, 2) = 0.0; 00434 outputValues(6, i0, 3) = -w; 00435 outputValues(6, i0, 4) = -y*w*w; 00436 outputValues(6, i0, 5) = -y*y*w*w*w; 00437 outputValues(6, i0, 6) = 0.0; 00438 outputValues(6, i0, 7) = -x*w*w ; 00439 outputValues(6, i0, 8) = -2.0*x*y*w*w*w; 00440 outputValues(6, i0, 9) = -3.0*x*y*y*w*w*w*w; 00441 00442 outputValues(7, i0, 0) = 0.0; 00443 outputValues(7, i0, 1) = -w; 00444 outputValues(7, i0, 2) = -y*w*w; 00445 outputValues(7, i0, 3) = 0.0; 00446 outputValues(7, i0, 4) = -x*w*w; 00447 outputValues(7, i0, 5) = -2.0*x*y*w*w*w; 00448 outputValues(7, i0, 6) = 0.0; 00449 outputValues(7, i0, 7) = 0.0; 00450 outputValues(7, i0, 8) = -x*x*w*w*w; 00451 outputValues(7, i0, 9) = -3.0*x*x*y*w*w*w*w; 00452 00453 outputValues(8, i0, 0) = 0.0; 00454 outputValues(8, i0, 1) = 0.0; 00455 outputValues(8, i0, 2) = 0.0; 00456 outputValues(8, i0, 3) = w; 00457 outputValues(8, i0, 4) = y*w*w; 00458 outputValues(8, i0, 5) = y*y*w*w*w; 00459 outputValues(8, i0, 6) = 0.0; 00460 outputValues(8, i0, 7) = x*w*w; 00461 outputValues(8, i0, 8) = 2.0*x*y*w*w*w; 00462 outputValues(8, i0, 9) = 3.0*x*y*y*w*w*w*w ; 00463 00464 outputValues(9, i0, 0) = 0.0; 00465 outputValues(9, i0, 1) = 0.0; 00466 outputValues(9, i0, 2) = 0.0; 00467 outputValues(9, i0, 3) = 0.0; 00468 outputValues(9, i0, 4) = w*w; 00469 outputValues(9, i0, 5) = 2.0*y*w*w*w; 00470 outputValues(9, i0, 6) = 0.0; 00471 outputValues(9, i0, 7) = 0.0; 00472 outputValues(9, i0, 8) = 2.0*x*w*w*w; 00473 outputValues(9, i0, 9) = 6.0*x*y*w*w*w*w; 00474 00475 outputValues(10,i0, 0) = 0.0; 00476 outputValues(10,i0, 1) = 0.0; 00477 outputValues(10,i0, 2) = 0.0; 00478 outputValues(10,i0, 3) = 0.0; 00479 outputValues(10,i0, 4) = -w*w; 00480 outputValues(10,i0, 5) = -2.0*y*w*w*w; 00481 outputValues(10,i0, 6) = 0.0; 00482 outputValues(10,i0, 7) = 0.0; 00483 outputValues(10,i0, 8) = -2.0*x*w*w*w; 00484 outputValues(10,i0, 9) = -6.0*x*y*w*w*w*w; 00485 00486 outputValues(11,i0, 0) = 0.0; 00487 outputValues(11,i0, 1) = 0.0; 00488 outputValues(11,i0, 2) = 0.0; 00489 outputValues(11,i0, 3) = 0.0; 00490 outputValues(11,i0, 4) = w*w; 00491 outputValues(11,i0, 5) = 2.0*y*w*w*w; 00492 outputValues(11,i0, 6) = 0.0; 00493 outputValues(11,i0, 7) = 0.0; 00494 outputValues(11,i0, 8) = 2.0*x*w*w*w; 00495 outputValues(11,i0, 9) = 6.0*x*y*w*w*w*w; 00496 00497 outputValues(12,i0, 0) = 0.0; 00498 outputValues(12,i0, 1) = 0.0; 00499 outputValues(12,i0, 2) = 0.0; 00500 outputValues(12,i0, 3) = 0.0; 00501 outputValues(12,i0, 4) = -w*w; 00502 outputValues(12,i0, 5) = -2.0*y*w*w*w; 00503 outputValues(12,i0, 6) = 0.0; 00504 outputValues(12,i0, 7) = 0.0; 00505 outputValues(12,i0, 8) = -2.0*x*w*w*w; 00506 outputValues(12,i0, 9) = -6.0*x*y*w*w*w*w; 00507 00508 } 00509 break; 00510 00511 case OPERATOR_D4: 00512 case OPERATOR_D5: 00513 case OPERATOR_D6: 00514 case OPERATOR_D7: 00515 case OPERATOR_D8: 00516 case OPERATOR_D9: 00517 case OPERATOR_D10: 00518 { 00519 // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, DkCardinality) 00520 int DkCardinality = Intrepid::getDkCardinality(operatorType, 00521 this -> basisCellTopology_.getDimension() ); 00522 for(int dofOrd = 0; dofOrd < this -> basisCardinality_; dofOrd++) { 00523 for (int i0 = 0; i0 < dim0; i0++) { 00524 for(int dkOrd = 0; dkOrd < DkCardinality; dkOrd++){ 00525 outputValues(dofOrd, i0, dkOrd) = 0.0; 00526 } 00527 } 00528 } 00529 } 00530 break; 00531 00532 default: 00533 TEUCHOS_TEST_FOR_EXCEPTION( !( Intrepid::isValidOperator(operatorType) ), std::invalid_argument, 00534 ">>> ERROR (Basis_HGRAD_PYR_I2_FEM): Invalid operator type"); 00535 } 00536 } 00537 00538 00539 00540 template<class Scalar, class ArrayScalar> 00541 void Basis_HGRAD_PYR_I2_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar& outputValues, 00542 const ArrayScalar & inputPoints, 00543 const ArrayScalar & cellVertices, 00544 const EOperator operatorType) const { 00545 TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error, 00546 ">>> ERROR (Basis_HGRAD_PYR_I2_FEM): FEM Basis calling an FVD member function"); 00547 } 00548 }// namespace Intrepid 00549 #endif
1.7.6.1