|
Intrepid
|
00001 #ifndef INTREPID_HGRAD_TET_COMP12_FEMDEF_HPP 00002 #define INTREPID_HGRAD_TET_COMP12_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 00052 namespace Intrepid { 00053 00054 template<class Scalar, class ArrayScalar> 00055 Basis_HGRAD_TET_COMP12_FEM<Scalar, ArrayScalar>::Basis_HGRAD_TET_COMP12_FEM() 00056 { 00057 this -> basisCardinality_ = 10; 00058 this -> basisDegree_ = 1; 00059 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<11> >() ); 00060 this -> basisType_ = BASIS_FEM_DEFAULT; 00061 this -> basisCoordinates_ = COORDINATES_CARTESIAN; 00062 this -> basisTagsAreSet_ = false; 00063 } 00064 00065 00066 template<class Scalar, class ArrayScalar> 00067 void Basis_HGRAD_TET_COMP12_FEM<Scalar, ArrayScalar>::initializeTags() { 00068 00069 // Basis-dependent intializations 00070 int tagSize = 4; // size of DoF tag 00071 int posScDim = 0; // position in the tag, counting from 0, of the subcell dim 00072 int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal 00073 int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell 00074 00075 // An array with local DoF tags assigned to basis functions, in the order of their local enumeration 00076 int tags[] = { 0, 0, 0, 1, 00077 0, 1, 0, 1, 00078 0, 2, 0, 1, 00079 0, 3, 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 }; 00087 00088 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays: 00089 Intrepid::setOrdinalTagData(this -> tagToOrdinal_, 00090 this -> ordinalToTag_, 00091 tags, 00092 this -> basisCardinality_, 00093 tagSize, 00094 posScDim, 00095 posScOrd, 00096 posDfOrd); 00097 } 00098 00099 00100 00101 template<class Scalar, class ArrayScalar> 00102 void Basis_HGRAD_TET_COMP12_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar & outputValues, 00103 const ArrayScalar & inputPoints, 00104 const EOperator operatorType) const { 00105 00106 // Verify arguments 00107 #ifdef HAVE_INTREPID_DEBUG 00108 Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues, 00109 inputPoints, 00110 operatorType, 00111 this -> getBaseCellTopology(), 00112 this -> getCardinality() ); 00113 #endif 00114 00115 // Number of evaluation points = dim 0 of inputPoints 00116 int dim0 = inputPoints.dimension(0); 00117 00118 // Temporaries: (x,y,z) coordinates of the evaluation point 00119 Scalar x = 0.0; 00120 Scalar y = 0.0; 00121 Scalar z = 0.0; 00122 00123 // Temporary for the auxiliary node basis function 00124 Scalar aux = 0.0; 00125 00126 // Array to store all the subtets containing the given pt 00127 Teuchos::Array<int> pt_tets; 00128 00129 switch (operatorType) { 00130 00131 case OPERATOR_VALUE: 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 pt_tets = getLocalSubTetrahedra(x,y,z); 00138 00139 // idependent verification shows that a given point will produce 00140 // the same shape functions for each tet that contains it, so 00141 // we only need to use the first one returned. 00142 //for (int pt = 0; pt < pt_tets.size(); ++pt) 00143 int subtet = pt_tets[0]; //pt_tets[pt]; 00144 aux = 0.0; 00145 // outputValues is a rank-2 array with dimensions (basisCardinality_, dim0) 00146 switch (subtet) { 00147 case 0: 00148 outputValues(0, i0) = 1. - 2. * (x + y + z); 00149 outputValues(4, i0) = 2. * x; 00150 outputValues(6, i0) = 2. * y; 00151 outputValues(7, i0) = 2. * z; 00152 break; 00153 case 1: 00154 outputValues(1, i0) = 2. * x - 1.; 00155 outputValues(4, i0) = 2. - 2. * (x + y + z); 00156 outputValues(5, i0) = 2. * y; 00157 outputValues(8, i0) = 2. * z; 00158 break; 00159 case 2: 00160 outputValues(2, i0) = 2. * y - 1.; 00161 outputValues(5, i0) = 2. * x; 00162 outputValues(6, i0) = 2. - 2. * (x + y + z); 00163 outputValues(9, i0) = 2. * z; 00164 break; 00165 case 3: 00166 outputValues(3, i0) = 2. * z - 1.; 00167 outputValues(7, i0) = 2. - 2. * (x + y + z); 00168 outputValues(8, i0) = 2. * x; 00169 outputValues(9, i0) = 2. * y; 00170 break; 00171 case 4: 00172 outputValues(4, i0) = 1. - 2. * (y + z); 00173 outputValues(5, i0) = 2. * (x + y) - 1.; 00174 outputValues(8, i0) = 2. * (x + z) - 1.; 00175 aux = 2. - 4. * x; 00176 break; 00177 case 5: 00178 outputValues(5, i0) = 1. - 2. * (x + y); 00179 outputValues(8, i0) = 1. - 2. * (x + z); 00180 outputValues(9, i0) = 2. * (y + z) - 1.; 00181 aux = 4. - 4. * (x + y + z); 00182 break; 00183 case 6: 00184 outputValues(7, i0) = 1. - 2. * (x + y); 00185 outputValues(8, i0) = 2. * (x + z) - 1.; 00186 outputValues(9, i0) = 2. * (y + z) - 1.; 00187 aux = 2. - 4. * z; 00188 break; 00189 case 7: 00190 outputValues(4, i0) = 1. - 2. * (y + z); 00191 outputValues(7, i0) = 1. - 2. * (x + y); 00192 outputValues(8, i0) = 2. * (x + z) - 1.; 00193 aux = 4. * y; 00194 break; 00195 case 8: 00196 outputValues(4, i0) = 1. - 2. * (y + z); 00197 outputValues(5, i0) = 2. * (x + y) - 1.; 00198 outputValues(6, i0) = 1. - 2. * (x + z); 00199 aux = 4. * z; 00200 break; 00201 case 9: 00202 outputValues(5, i0) = 2. * (x + y) - 1.; 00203 outputValues(6, i0) = 1. - 2. * (x + z); 00204 outputValues(9, i0) = 2. * (y + z) - 1.; 00205 aux = 2. - 4. * y; 00206 break; 00207 case 10: 00208 outputValues(6, i0) = 1. - 2. * (x + z); 00209 outputValues(7, i0) = 1. - 2. * (x + y); 00210 outputValues(9, i0) = 2. * (y + z) - 1.; 00211 aux = 4. * x; 00212 break; 00213 case 11: 00214 outputValues(4, i0) = 1. - 2. * (y + z); 00215 outputValues(6, i0) = 1. - 2. * (x + z); 00216 outputValues(7, i0) = 1. - 2. * (x + y); 00217 aux = 4. * (x + y + z) - 2.; 00218 break; 00219 } 00220 outputValues(4, i0) += aux/6.0; 00221 outputValues(5, i0) += aux/6.0; 00222 outputValues(6, i0) += aux/6.0; 00223 outputValues(7, i0) += aux/6.0; 00224 outputValues(8, i0) += aux/6.0; 00225 outputValues(9, i0) += aux/6.0; 00226 //} 00227 } 00228 break; 00229 00230 case OPERATOR_GRAD: 00231 case OPERATOR_D1: 00232 { 00233 // initialize to 0.0 since we will be accumulating 00234 outputValues.initialize(0.0); 00235 00236 // local looping indices 00237 // NOTE: bc is barycentric coord 00238 int node, tet, pt, dim, bc1, bc2; 00239 00240 // get 12 subtet gradients sized as 3, 11, 12 00241 Intrepid::FieldContainer<Scalar> dx = getSubTetGrads(); 00242 // get subtet jacobian determinants (12) 00243 Intrepid::FieldContainer<Scalar> xJ = getSubTetDetF(); 00244 // get integration weights (numPoints, 12) 00245 Intrepid::FieldContainer<Scalar> w = getWeights(inputPoints); 00246 // get Barycentric Coordinates of points 00247 Intrepid::FieldContainer<Scalar> lambda = getBarycentricCoords(inputPoints); 00248 00249 // declare FieldContainer for intermediate calcs 00250 // a - for the 4x4 projection matrix 00251 FieldContainer<Scalar> a(4, 4); 00252 FieldContainer<Scalar> ai; 00253 00254 // b - for the 3 x 4 x 10 integrated gradients 00255 FieldContainer<Scalar> b(3,4,10); 00256 00257 // c - product of inv(a) and b 00258 FieldContainer<Scalar> c(4,3,10); 00259 00260 // again initialize for safety 00261 a.initialize(0.0); 00262 b.initialize(0.0); 00263 c.initialize(0.0); 00264 00265 for (pt=0; pt < dim0; ++pt) 00266 { 00267 // compute a 00268 for (tet=0; tet < 12; ++tet) 00269 for (bc1=0; bc1 < 4; ++bc1) 00270 for (bc2=0; bc2 < 4; ++bc2) 00271 a(bc1,bc2) += xJ(tet) * w(pt,tet) * lambda(pt,bc1) * lambda(pt,bc2); 00272 00273 // compute b 00274 for (tet=0; tet < 12; ++tet) 00275 { 00276 for (node=0; node < 10; ++node) 00277 for (bc1=0; bc1 < 4; ++bc1) 00278 for (dim=0; dim < 3; ++dim) 00279 b(dim,bc1,node) += xJ(tet) * w(pt,tet) * lambda(pt,bc1) * dx(dim,node,tet); 00280 00281 // add in contribution from the auxiliary node, averaged over the 6 mid-edge ndoes 00282 for (node=4; node < 10; ++node) 00283 for (bc1=0; bc1 < 4; ++bc1) 00284 for (dim=0; dim < 3; ++dim) 00285 b(dim,bc1,node) += xJ(tet) * w(pt,tet) * lambda(pt,bc1) * dx(dim,10,tet)/6; 00286 } 00287 } // loop over pts 00288 00289 // compute inverse of a (4x4 inverse) 00290 ai = inverse44(a); 00291 00292 // compute c 00293 for (node=0; node < 10; ++node) 00294 for (bc1=0; bc1 < 4; ++bc1) 00295 for (bc2=0; bc2 < 4; ++bc2) 00296 for (dim=0; dim < 3; ++dim) 00297 c(bc1,dim,node) += ai(bc1,bc2) * b(dim,bc2,node); 00298 00299 // fill in shape function derivatives 00300 for (pt=0; pt < dim0; ++pt) 00301 for (node=0; node < 10; ++node) 00302 for (dim=0; dim < 3; ++dim) 00303 for (bc1=0; bc1 < 4; ++bc1) 00304 outputValues(node,pt,dim) += lambda(pt,bc1) * c(bc1,dim,node); 00305 } 00306 break; 00307 00308 case OPERATOR_CURL: 00309 TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_CURL), std::invalid_argument, 00310 ">>> ERROR (Basis_HGRAD_TET_COMP12_FEM): CURL is invalid operator for rank-0 (scalar) functions in 3D"); 00311 break; 00312 00313 case OPERATOR_DIV: 00314 TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument, 00315 ">>> ERROR (Basis_HGRAD_TET_COMP12_FEM): DIV is invalid operator for rank-0 (scalar) functions in 3D"); 00316 break; 00317 00318 case OPERATOR_D2: 00319 case OPERATOR_D3: 00320 case OPERATOR_D4: 00321 case OPERATOR_D5: 00322 case OPERATOR_D6: 00323 case OPERATOR_D7: 00324 case OPERATOR_D8: 00325 case OPERATOR_D9: 00326 case OPERATOR_D10: 00327 { 00328 // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, DkCardinality) 00329 int DkCardinality = Intrepid::getDkCardinality(operatorType, 00330 this -> basisCellTopology_.getDimension() ); 00331 for(int dofOrd = 0; dofOrd < this -> basisCardinality_; dofOrd++) { 00332 for (int i0 = 0; i0 < dim0; i0++) { 00333 for(int dkOrd = 0; dkOrd < DkCardinality; dkOrd++){ 00334 outputValues(dofOrd, i0, dkOrd) = 0.0; 00335 } 00336 } 00337 } 00338 } 00339 break; 00340 default: 00341 TEUCHOS_TEST_FOR_EXCEPTION( !( Intrepid::isValidOperator(operatorType) ), std::invalid_argument, 00342 ">>> ERROR (Basis_HGRAD_TET_COMP12_FEM): Invalid operator type"); 00343 } 00344 } 00345 00346 00347 00348 template<class Scalar, class ArrayScalar> 00349 void Basis_HGRAD_TET_COMP12_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar& outputValues, 00350 const ArrayScalar & inputPoints, 00351 const ArrayScalar & cellVertices, 00352 const EOperator operatorType) const { 00353 TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error, 00354 ">>> ERROR (Basis_HGRAD_TET_COMP12_FEM): FEM Basis calling an FVD member function"); 00355 00356 } 00357 00358 template<class Scalar, class ArrayScalar> 00359 void Basis_HGRAD_TET_COMP12_FEM<Scalar, ArrayScalar>::getDofCoords(ArrayScalar & DofCoords) const { 00360 #ifdef HAVE_INTREPID_DEBUG 00361 // Verify rank of output array. 00362 TEUCHOS_TEST_FOR_EXCEPTION( !(DofCoords.rank() == 2), std::invalid_argument, 00363 ">>> ERROR: (Intrepid::Basis_HGRAD_TET_COMP12_FEM::getDofCoords) rank = 2 required for DofCoords array"); 00364 // Verify 0th dimension of output array. 00365 TEUCHOS_TEST_FOR_EXCEPTION( !( DofCoords.dimension(0) == this -> basisCardinality_ ), std::invalid_argument, 00366 ">>> ERROR: (Intrepid::Basis_HGRAD_TET_COMP12_FEM::getDofCoords) mismatch in number of DoF and 0th dimension of DofCoords array"); 00367 // Verify 1st dimension of output array. 00368 TEUCHOS_TEST_FOR_EXCEPTION( !( DofCoords.dimension(1) == (int)(this -> basisCellTopology_.getDimension()) ), std::invalid_argument, 00369 ">>> ERROR: (Intrepid::Basis_HGRAD_TET_COMP12_FEM::getDofCoords) incorrect reference cell (1st) dimension in DofCoords array"); 00370 #endif 00371 00372 DofCoords(0,0) = 0.0; DofCoords(0,1) = 0.0; DofCoords(0,2) = 0.0; 00373 DofCoords(1,0) = 1.0; DofCoords(1,1) = 0.0; DofCoords(1,2) = 0.0; 00374 DofCoords(2,0) = 0.0; DofCoords(2,1) = 1.0; DofCoords(2,2) = 0.0; 00375 DofCoords(3,0) = 0.0; DofCoords(3,1) = 0.0; DofCoords(3,2) = 1.0; 00376 DofCoords(4,0) = 0.5; DofCoords(4,1) = 0.0; DofCoords(4,2) = 0.0; 00377 DofCoords(5,0) = 0.5; DofCoords(5,1) = 0.5; DofCoords(5,2) = 0.0; 00378 DofCoords(6,0) = 0.0; DofCoords(6,1) = 0.5; DofCoords(6,2) = 0.0; 00379 DofCoords(7,0) = 0.0; DofCoords(7,1) = 0.0; DofCoords(7,2) = 0.5; 00380 DofCoords(8,0) = 0.5; DofCoords(8,1) = 0.0; DofCoords(8,2) = 0.5; 00381 DofCoords(9,0) = 0.0; DofCoords(9,1) = 0.5; DofCoords(9,2) = 0.5; 00382 } 00383 00384 template<class Scalar, class ArrayScalar> 00385 Teuchos::Array<int> 00386 Basis_HGRAD_TET_COMP12_FEM<Scalar, ArrayScalar>::getLocalSubTetrahedra(Scalar x, Scalar y, Scalar z) const 00387 { 00388 00389 Teuchos::Array<int> subTets; 00390 00391 // local coords 00392 Scalar xyz = x + y + z; 00393 Scalar xy = x + y; 00394 Scalar xz = x + z; 00395 Scalar yz = y + z; 00396 00397 // cycle through each subdomain and push back if the point lies within 00398 00399 // subtet #0 E0 := 0.0 <= r + s + t <= 0.5 && 0.0 <= r <= 0.5 && 0.0 <= s <= 0.5 && 0.0 <= t <= 0.5 00400 if ( (0.0 <= xyz && xyz <= 0.5) && (0.0 <= x && x <= 0.5) && (0.0 <= y && y <= 0.5) && (0.0 <= z && z <= 0.5) ) 00401 subTets.push_back(0); 00402 00403 // subtet #1 E1 := 0.5 <= r + s + t <= 1.0 && 0.5 <= r <= 1.0 && 0.0 <= s <= 0.5 && 0.0 <= t <= 0.5 00404 if ( (0.5 <= xyz && xyz <= 1.0) && (0.5 <= x && x <= 1.0) && (0.0 <= y && y <= 0.5) && (0.0 <= z && z <= 0.5) ) 00405 subTets.push_back(1); 00406 00407 // subtet #2 E2 := 0.5 <= r + s + t <= 1.0 && 0.0 <= r <= 0.5 && 0.5 <= s <= 1.0 && 0.0 <= t <= 0.5 00408 if ( (0.5 <= xyz && xyz <= 1.0) && (0.0 <= x && x <= 0.5) && (0.5 <= y && y <= 1.0) && (0.0 <= z && z<= 0.5) ) 00409 subTets.push_back(2); 00410 00411 // subtet #3 E3 := 0.5 <= r + s + t <= 1.0 && 0.0 <= r <= 0.5 && 0.0 <= s <= 0.5 && 0.5 <= t <= 1.0 00412 if ( (0.5 <= xyz && xyz <= 1.0) && (0.0 <= x && x <= 0.5) && (0.0 <= y && y <= 0.5) && (0.5 <= z && z <= 1.0) ) 00413 subTets.push_back(3); 00414 00415 // subtet #4 E4 := 0.0 <= s + t <= 0.5 && 0.5 <= r + s <= 1.0 && 0.5 <= r + t <= 1.0 && 0.0 <= r <= 0.5 00416 if ( (0.0 <= yz && yz <= 0.5) && (0.5 <= xy && xy <= 1.0) && (0.5 <= xz && xz <= 1.0) && (0.0 <= x && x <= 0.5) ) 00417 subTets.push_back(4); 00418 00419 // subtet #5 E5 := 0.5 <= r + s <= 1.0 && 0.5 <= s + t <= 1.0 && 0.5 <= r + t <= 1.0 && 0.75 <= r + s + t <= 1.0 00420 if ( (0.5 <= xy && xy <= 1.0) && (0.5 <= yz && yz <= 1.0) && (0.5 <= xz && xz <= 1.0) && (0.75 <= xyz && xyz <= 1.0) ) 00421 subTets.push_back(5); 00422 00423 // subtet #6 E6 := 0.5 <= s + t <= 1.0 && 0.0 <= r + s <= 0.5 && 0.5 <= r + t <= 1.0 && 0.0 <= t <= 0.5 00424 if ( (0.5 <= yz && yz <= 1.0) && (0.0 <= xy && xy <= 0.5) && (0.5 <= xz && xz <= 1.0) && (0.0 <= z && z <= 0.5) ) 00425 subTets.push_back(6); 00426 00427 // subtet #7 E7 := 0.0 <= s + t <= 0.5 && 0.0 <= r + s <= 0.5 && 0.5 <= r + t <= 1.0 && 0.0 <= s <= 0.25 00428 if ( (0.0 <= yz && yz <= 0.5) && (0.0 <= xy && xy <= 0.5) && (0.5 <= xz && xz <= 1.0) && (0.0 <= y && y <= 0.25) ) 00429 subTets.push_back(7); 00430 00431 // subtet #8 E8 := 0.0 <= r + t <= 0.5 && 0.0 <= s + t <= 0.5 && 0.5 <= r + s <= 1.0 && 0.0 <= t <= 0.25 00432 if ( (0.0 <= xz && xz <= 0.5) && (0.0 <= yz && yz <= 0.5) && (0.5 <= xy && xy <= 1.0) && (0.0 <= z && z <= 0.5) ) 00433 subTets.push_back(8); 00434 00435 // subtet #9 E9 := 0.0 <= r + t <= 0.5 && 0.5 <= r + s <= 1.0 && 0.5 <= s + t <= 1.0 && 0.0 <= s <= 0.5 00436 if ( (0.0 <= xz && xz <= 0.5) && (0.5 <= xy && xy <= 1.0) && (0.5 <= yz && yz <= 1.0) && (0.0 <= y && y <= 0.5) ) 00437 subTets.push_back(9); 00438 00439 // subtet #10 E10 := 0.0 <= r + t <= 0.5 && 0.5 <= s + t <= 1.0 && 0.0 <= r + s <= 0.5 && 0.0 <= r <= 0.25 00440 if ( (0.0 <= xz && xz <= 0.5) && (0.5 <= yz && yz <= 1.0) && (0.0 <= xy && xy <= 0.5) && (0.0 <= x && x <= 0.25) ) 00441 subTets.push_back(10); 00442 00443 // subtet #11 E11 := 0.5 <= r + s + t <= 0.75 && 0.0 <= r + t <= 0.5 && 0.0 <= s + t <= 0.5 && 0.0 <= r + s <= 0.5 00444 if ( (0.5 <= xyz && xyz <= 0.75) && (0.0 <= xz && xz <= 0.5) && (0.0 <= yz && yz <= 0.5) && (0.0 <= xy && xy <= 0.5) ) 00445 subTets.push_back(11); 00446 00447 return subTets; 00448 } 00449 00450 00451 template<class Scalar, class ArrayScalar> 00452 Intrepid::FieldContainer<Scalar> 00453 Basis_HGRAD_TET_COMP12_FEM<Scalar, ArrayScalar>::getWeights(const ArrayScalar & inPts) const 00454 { 00455 int numPoints = inPts.dimension(0); 00456 Intrepid::FieldContainer<Scalar> w(numPoints, 12); 00457 w.initialize(0.0); 00458 Teuchos::Array< Teuchos::Array<int> > pt_tets; 00459 00460 for (int pt = 0; pt < numPoints; ++pt) 00461 pt_tets.push_back(getLocalSubTetrahedra(inPts(pt,0), inPts(pt,1), inPts(pt,2))); 00462 00463 Teuchos::Array<int> flat; 00464 FieldContainer<int> count(12); 00465 00466 for (int pt = 0; pt < numPoints; ++pt) 00467 for (int i = 0; i < pt_tets[pt].size(); ++i) 00468 flat.push_back(pt_tets[pt][i]); 00469 00470 for (int i = 0; i < flat.size(); ++i) 00471 count(flat[i])++; 00472 00473 for (int pt = 0; pt < numPoints; ++pt) 00474 for (int i = 0; i < pt_tets[pt].size(); ++i) 00475 w(pt, pt_tets[pt][i]) = 1.0/count(pt_tets[pt][i]); 00476 00477 return w; 00478 } 00479 00480 00481 00482 template<class Scalar, class ArrayScalar> 00483 Intrepid::FieldContainer<Scalar> 00484 Basis_HGRAD_TET_COMP12_FEM<Scalar, ArrayScalar>::getBarycentricCoords(const ArrayScalar & inPts) const 00485 { 00486 int numPoints = inPts.dimension(0); 00487 Intrepid::FieldContainer<Scalar> lambda(numPoints, 4); 00488 00489 for (int pt = 0; pt < numPoints; ++pt) 00490 { 00491 lambda(pt,0) = 1. - inPts(pt,0) - inPts(pt,1) - inPts(pt,2); 00492 lambda(pt,1) = inPts(pt,0); 00493 lambda(pt,2) = inPts(pt,1); 00494 lambda(pt,3) = inPts(pt,2); 00495 } 00496 00497 return lambda; 00498 } 00499 00500 template<class Scalar, class ArrayScalar> 00501 Scalar 00502 Basis_HGRAD_TET_COMP12_FEM<Scalar, ArrayScalar>::det44(const Intrepid::FieldContainer<Scalar> a) const 00503 { 00504 Scalar det = a(0,3) * a(1,2) * a(2,1) * a(3,0) 00505 - a(0,2) * a(1,3) * a(2,1) * a(3,0) 00506 - a(0,3) * a(1,1) * a(2,2) * a(3,0) 00507 + a(0,1) * a(1,3) * a(2,2) * a(3,0) 00508 + a(0,2) * a(1,1) * a(2,3) * a(3,0) 00509 - a(0,1) * a(1,2) * a(2,3) * a(3,0) 00510 - a(0,3) * a(1,2) * a(2,0) * a(3,1) 00511 + a(0,2) * a(1,3) * a(2,0) * a(3,1) 00512 + a(0,3) * a(1,0) * a(2,2) * a(3,1) 00513 - a(0,0) * a(1,3) * a(2,2) * a(3,1) 00514 - a(0,2) * a(1,0) * a(2,3) * a(3,1) 00515 + a(0,0) * a(1,2) * a(2,3) * a(3,1) 00516 + a(0,3) * a(1,1) * a(2,0) * a(3,2) 00517 - a(0,1) * a(1,3) * a(2,0) * a(3,2) 00518 - a(0,3) * a(1,0) * a(2,1) * a(3,2) 00519 + a(0,0) * a(1,3) * a(2,1) * a(3,2) 00520 + a(0,1) * a(1,0) * a(2,3) * a(3,2) 00521 - a(0,0) * a(1,1) * a(2,3) * a(3,2) 00522 - a(0,2) * a(1,1) * a(2,0) * a(3,3) 00523 + a(0,1) * a(1,2) * a(2,0) * a(3,3) 00524 + a(0,2) * a(1,0) * a(2,1) * a(3,3) 00525 - a(0,0) * a(1,2) * a(2,1) * a(3,3) 00526 - a(0,1) * a(1,0) * a(2,2) * a(3,3) 00527 + a(0,0) * a(1,1) * a(2,2) * a(3,3); 00528 00529 return det; 00530 } 00531 00532 template<class Scalar, class ArrayScalar> 00533 Intrepid::FieldContainer<Scalar> 00534 Basis_HGRAD_TET_COMP12_FEM<Scalar, ArrayScalar>::inverse44(const Intrepid::FieldContainer<Scalar> a) const 00535 { 00536 Intrepid::FieldContainer<Scalar> ai(4,4); 00537 Scalar xj = det44(a); 00538 00539 ai(0,0) = (1/xj) * (-a(1,3) * a(2,2) * a(3,1) + a(1,2) * a(2,3) * a(3,1) + a(1,3) * a(2,1) * a(3,2) - a(1,1) * a(2,3) * a(3,2) - a(1,2) * a(2,1) * a(3,3) + a(1,1) * a(2,2) * a(3,3)); 00540 ai(0,1) = (1/xj) * ( a(0,3) * a(2,2) * a(3,1) - a(0,2) * a(2,3) * a(3,1) - a(0,3) * a(2,1) * a(3,2) + a(0,1) * a(2,3) * a(3,2) + a(0,2) * a(2,1) * a(3,3) - a(0,1) * a(2,2) * a(3,3)); 00541 ai(0,2) = (1/xj) * (-a(0,3) * a(1,2) * a(3,1) + a(0,2) * a(1,3) * a(3,1) + a(0,3) * a(1,1) * a(3,2) - a(0,1) * a(1,3) * a(3,2) - a(0,2) * a(1,1) * a(3,3) + a(0,1) * a(1,2) * a(3,3)); 00542 ai(0,3) = (1/xj) * ( a(0,3) * a(1,2) * a(2,1) - a(0,2) * a(1,3) * a(2,1) - a(0,3) * a(1,1) * a(2,2) + a(0,1) * a(1,3) * a(2,2) + a(0,2) * a(1,1) * a(2,3) - a(0,1) * a(1,2) * a(2,3)); 00543 00544 ai(1,0) = (1/xj) * ( a(1,3) * a(2,2) * a(3,0) - a(1,2) * a(2,3) * a(3,0) - a(1,3) * a(2,0) * a(3,2) + a(1,0) * a(2,3) * a(3,2) + a(1,2) * a(2,0) * a(3,3) - a(1,0) * a(2,2) * a(3,3)); 00545 ai(1,1) = (1/xj) * (-a(0,3) * a(2,2) * a(3,0) + a(0,2) * a(2,3) * a(3,0) + a(0,3) * a(2,0) * a(3,2) - a(0,0) * a(2,3) * a(3,2) - a(0,2) * a(2,0) * a(3,3) + a(0,0) * a(2,2) * a(3,3)); 00546 ai(1,2) = (1/xj) * ( a(0,3) * a(1,2) * a(3,0) - a(0,2) * a(1,3) * a(3,0) - a(0,3) * a(1,0) * a(3,2) + a(0,0) * a(1,3) * a(3,2) + a(0,2) * a(1,0) * a(3,3) - a(0,0) * a(1,2) * a(3,3)); 00547 ai(1,3) = (1/xj) * (-a(0,3) * a(1,2) * a(2,0) + a(0,2) * a(1,3) * a(2,0) + a(0,3) * a(1,0) * a(2,2) - a(0,0) * a(1,3) * a(2,2) - a(0,2) * a(1,0) * a(2,3) + a(0,0) * a(1,2) * a(2,3)); 00548 00549 ai(2,0) = (1/xj) * (-a(1,3) * a(2,1) * a(3,0) + a(1,1) * a(2,3) * a(3,0) + a(1,3) * a(2,0) * a(3,1) - a(1,0) * a(2,3) * a(3,1) - a(1,1) * a(2,0) * a(3,3) + a(1,0) * a(2,1) * a(3,3)); 00550 ai(2,1) = (1/xj) * ( a(0,3) * a(2,1) * a(3,0) - a(0,1) * a(2,3) * a(3,0) - a(0,3) * a(2,0) * a(3,1) + a(0,0) * a(2,3) * a(3,1) + a(0,1) * a(2,0) * a(3,3) - a(0,0) * a(2,1) * a(3,3)); 00551 ai(2,2) = (1/xj) * (-a(0,3) * a(1,1) * a(3,0) + a(0,1) * a(1,3) * a(3,0) + a(0,3) * a(1,0) * a(3,1) - a(0,0) * a(1,3) * a(3,1) - a(0,1) * a(1,0) * a(3,3) + a(0,0) * a(1,1) * a(3,3)); 00552 ai(2,3) = (1/xj) * ( a(0,3) * a(1,1) * a(2,0) - a(0,1) * a(1,3) * a(2,0) - a(0,3) * a(1,0) * a(2,1) + a(0,0) * a(1,3) * a(2,1) + a(0,1) * a(1,0) * a(2,3) - a(0,0) * a(1,1) * a(2,3)); 00553 00554 ai(3,0) = (1/xj) * ( a(1,2) * a(2,1) * a(3,0) - a(1,1) * a(2,2) * a(3,0) - a(1,2) * a(2,0) * a(3,1) + a(1,0) * a(2,2) * a(3,1) + a(1,1) * a(2,0) * a(3,2) - a(1,0) * a(2,1) * a(3,2)); 00555 ai(3,1) = (1/xj) * (-a(0,2) * a(2,1) * a(3,0) + a(0,1) * a(2,2) * a(3,0) + a(0,2) * a(2,0) * a(3,1) - a(0,0) * a(2,2) * a(3,1) - a(0,1) * a(2,0) * a(3,2) + a(0,0) * a(2,1) * a(3,2)); 00556 ai(3,2) = (1/xj) * ( a(0,2) * a(1,1) * a(3,0) - a(0,1) * a(1,2) * a(3,0) - a(0,2) * a(1,0) * a(3,1) + a(0,0) * a(1,2) * a(3,1) + a(0,1) * a(1,0) * a(3,2) - a(0,0) * a(1,1) * a(3,2)); 00557 ai(3,3) = (1/xj) * (-a(0,2) * a(1,1) * a(2,0) + a(0,1) * a(1,2) * a(2,0) + a(0,2) * a(1,0) * a(2,1) - a(0,0) * a(1,2) * a(2,1) - a(0,1) * a(1,0) * a(2,2) + a(0,0) * a(1,1) * a(2,2)); 00558 00559 return ai; 00560 } 00561 00562 template<class Scalar, class ArrayScalar> 00563 Intrepid::FieldContainer<Scalar> 00564 Basis_HGRAD_TET_COMP12_FEM<Scalar, ArrayScalar>::getSubTetGrads() const 00565 { 00566 Intrepid::FieldContainer<Scalar> dx(3,11,12); 00567 dx.initialize(0.0); 00568 // fill in dx 00569 dx(0,0,0) = -2.0; 00570 dx(0,1,1) = 2.0; 00571 dx(0,4,0) = 2.0; 00572 dx(0,4,1) = -2.0; 00573 dx(0,5,2) = 2.0; 00574 dx(0,5,4) = 2.0; 00575 dx(0,5,5) = 2.0; 00576 dx(0,5,8) = 2.0; 00577 dx(0,5,9) = 2.0; 00578 dx(0,6,2) = -2.0; 00579 dx(0,6,8) = -2.0; 00580 dx(0,6,9) = -2.0; 00581 dx(0,6,10) = -2.0; 00582 dx(0,6,11) = -2.0; 00583 dx(0,7,3) = -2.0; 00584 dx(0,7,6) = -2.0; 00585 dx(0,7,7) = -2.0; 00586 dx(0,7,10) = -2.0; 00587 dx(0,7,11) = -2.0; 00588 dx(0,8,3) = 2.0; 00589 dx(0,8,4) = 2.0; 00590 dx(0,8,5) = 2.0; 00591 dx(0,8,6) = 2.0; 00592 dx(0,8,7) = 2.0; 00593 dx(0,10,4) = -4.0; 00594 dx(0,10,5) = -4.0; 00595 dx(0,10,10) = 4.0; 00596 dx(0,10,11) = 4.0; 00597 00598 // fill in dy 00599 dx(1,0,0) = -2.0; 00600 dx(1,2,2) = 2.0; 00601 dx(1,4,1) = -2.0; 00602 dx(1,4,4) = -2.0; 00603 dx(1,4,7) = -2.0; 00604 dx(1,4,8) = -2.0; 00605 dx(1,4,11) = -2.0; 00606 dx(1,5,1) = 2.0; 00607 dx(1,5,4) = 2.0; 00608 dx(1,5,5) = 2.0; 00609 dx(1,5,8) = 2.0; 00610 dx(1,5,9) = 2.0; 00611 dx(1,6,0) = 2.0; 00612 dx(1,6,2) = -2.0; 00613 dx(1,7,3) = -2.0; 00614 dx(1,7,6) = -2.0; 00615 dx(1,7,7) = -2.0; 00616 dx(1,7,10) = -2.0; 00617 dx(1,7,11) = -2.0; 00618 dx(1,9,3) = 2.0; 00619 dx(1,9,5) = 2.0; 00620 dx(1,9,6) = 2.0; 00621 dx(1,9,9) = 2.0; 00622 dx(1,9,10) = 2.0; 00623 dx(1,10,5) = -4.0; 00624 dx(1,10,7) = 4.0; 00625 dx(1,10,9) = -4.0; 00626 dx(1,10,11) = 4.0; 00627 00628 // fill in dz 00629 dx(2,0,0) = -2.0; 00630 dx(2,3,3) = 2.0; 00631 dx(2,4,1) = -2.0; 00632 dx(2,4,4) = -2.0; 00633 dx(2,4,7) = -2.0; 00634 dx(2,4,8) = -2.0; 00635 dx(2,4,11) = -2.0; 00636 dx(2,6,2) = -2.0; 00637 dx(2,6,8) = -2.0; 00638 dx(2,6,9) = -2.0; 00639 dx(2,6,10) = -2.0; 00640 dx(2,6,11) = -2.0; 00641 dx(2,7,0) = 2.0; 00642 dx(2,7,3) = -2.0; 00643 dx(2,8,1) = 2.0; 00644 dx(2,8,4) = 2.0; 00645 dx(2,8,5) = 2.0; 00646 dx(2,8,6) = 2.0; 00647 dx(2,8,7) = 2.0; 00648 dx(2,9,2) = 2.0; 00649 dx(2,9,5) = 2.0; 00650 dx(2,9,6) = 2.0; 00651 dx(2,9,9) = 2.0; 00652 dx(2,9,10) = 2.0; 00653 dx(2,10,5) = -4.0; 00654 dx(2,10,6) = -4.0; 00655 dx(2,10,8) = 4.0; 00656 dx(2,10,11) = 4.0; 00657 00658 return dx; 00659 } 00660 00661 template<class Scalar, class ArrayScalar> 00662 Intrepid::FieldContainer<Scalar> 00663 Basis_HGRAD_TET_COMP12_FEM<Scalar, ArrayScalar>::getSubTetDetF() const 00664 { 00665 Intrepid::FieldContainer<Scalar> xJ(12); 00666 // set sub-elem jacobians 00667 xJ(0) = 1./48.; xJ(1) = 1./48.; xJ(2) = 1./48.; xJ(3) = 1./48.; 00668 xJ(4) = 1./96.; xJ(5) = 1./96.; xJ(6) = 1./96.; xJ(7) = 1./96.; 00669 xJ(8) = 1./96.; xJ(9) = 1./96.; xJ(10) = 1./96.; xJ(11) = 1./96.; 00670 00671 return xJ; 00672 } 00673 00674 00675 }// namespace Intrepid 00676 #endif
1.7.6.1