|
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 00130 switch (operatorType) { 00131 00132 case OPERATOR_VALUE: 00133 for (int i0 = 0; i0 < dim0; i0++) { 00134 x = inputPoints(i0, 0); 00135 y = inputPoints(i0, 1); 00136 z = inputPoints(i0, 2); 00137 00138 pt_tets = getLocalSubTetrahedra(x,y,z); 00139 00140 // idependent verification shows that a given point will produce 00141 // the same shape functions for each tet that contains it, so 00142 // we only need to use the first one returned. 00143 //for (int pt = 0; pt < pt_tets.size(); ++pt) 00144 if (pt_tets[0] != -1) { 00145 int subtet = pt_tets[0]; 00146 aux = 0.0; 00147 // outputValues is a rank-2 array with dimensions (basisCardinality_, dim0) 00148 switch (subtet) { 00149 case 0: 00150 outputValues(0, i0) = 1. - 2. * (x + y + z); 00151 outputValues(4, i0) = 2. * x; 00152 outputValues(6, i0) = 2. * y; 00153 outputValues(7, i0) = 2. * z; 00154 break; 00155 case 1: 00156 outputValues(1, i0) = 2. * x - 1.; 00157 outputValues(4, i0) = 2. - 2. * (x + y + z); 00158 outputValues(5, i0) = 2. * y; 00159 outputValues(8, i0) = 2. * z; 00160 break; 00161 case 2: 00162 outputValues(2, i0) = 2. * y - 1.; 00163 outputValues(5, i0) = 2. * x; 00164 outputValues(6, i0) = 2. - 2. * (x + y + z); 00165 outputValues(9, i0) = 2. * z; 00166 break; 00167 case 3: 00168 outputValues(3, i0) = 2. * z - 1.; 00169 outputValues(7, i0) = 2. - 2. * (x + y + z); 00170 outputValues(8, i0) = 2. * x; 00171 outputValues(9, i0) = 2. * y; 00172 break; 00173 case 4: 00174 outputValues(4, i0) = 1. - 2. * (y + z); 00175 outputValues(5, i0) = 2. * (x + y) - 1.; 00176 outputValues(8, i0) = 2. * (x + z) - 1.; 00177 aux = 2. - 4. * x; 00178 break; 00179 case 5: 00180 outputValues(5, i0) = 1. - 2. * (x + y); 00181 outputValues(8, i0) = 1. - 2. * (x + z); 00182 outputValues(9, i0) = 2. * (y + z) - 1.; 00183 aux = 4. - 4. * (x + y + z); 00184 break; 00185 case 6: 00186 outputValues(7, i0) = 1. - 2. * (x + y); 00187 outputValues(8, i0) = 2. * (x + z) - 1.; 00188 outputValues(9, i0) = 2. * (y + z) - 1.; 00189 aux = 2. - 4. * z; 00190 break; 00191 case 7: 00192 outputValues(4, i0) = 1. - 2. * (y + z); 00193 outputValues(7, i0) = 1. - 2. * (x + y); 00194 outputValues(8, i0) = 2. * (x + z) - 1.; 00195 aux = 4. * y; 00196 break; 00197 case 8: 00198 outputValues(4, i0) = 1. - 2. * (y + z); 00199 outputValues(5, i0) = 2. * (x + y) - 1.; 00200 outputValues(6, i0) = 1. - 2. * (x + z); 00201 aux = 4. * z; 00202 break; 00203 case 9: 00204 outputValues(5, i0) = 2. * (x + y) - 1.; 00205 outputValues(6, i0) = 1. - 2. * (x + z); 00206 outputValues(9, i0) = 2. * (y + z) - 1.; 00207 aux = 2. - 4. * y; 00208 break; 00209 case 10: 00210 outputValues(6, i0) = 1. - 2. * (x + z); 00211 outputValues(7, i0) = 1. - 2. * (x + y); 00212 outputValues(9, i0) = 2. * (y + z) - 1.; 00213 aux = 4. * x; 00214 break; 00215 case 11: 00216 outputValues(4, i0) = 1. - 2. * (y + z); 00217 outputValues(6, i0) = 1. - 2. * (x + z); 00218 outputValues(7, i0) = 1. - 2. * (x + y); 00219 aux = 4. * (x + y + z) - 2.; 00220 break; 00221 } 00222 outputValues(4, i0) += aux/6.0; 00223 outputValues(5, i0) += aux/6.0; 00224 outputValues(6, i0) += aux/6.0; 00225 outputValues(7, i0) += aux/6.0; 00226 outputValues(8, i0) += aux/6.0; 00227 outputValues(9, i0) += aux/6.0; 00228 } 00229 } 00230 break; 00231 00232 case OPERATOR_GRAD: 00233 case OPERATOR_D1: 00234 { 00235 // initialize to 0.0 since we will be accumulating 00236 outputValues.initialize(0.0); 00237 00238 FieldContainer<Scalar> Lopt(10,3); 00239 for (int pt=0; pt < dim0; ++pt) { 00240 00241 Scalar r = inputPoints(pt, 0); 00242 Scalar s = inputPoints(pt, 1); 00243 Scalar t = inputPoints(pt, 2); 00244 00245 Lopt(0,0) = (-17 + 20*r + 20*s + 20*t)/8.; 00246 Lopt(0,1) = (-17 + 20*r + 20*s + 20*t)/8.; 00247 Lopt(0,2) = (-17 + 20*r + 20*s + 20*t)/8.; 00248 Lopt(1,0) = -0.375 + (5*r)/2.; 00249 Lopt(1,1) = 0.; 00250 Lopt(1,2) = 0.; 00251 Lopt(2,0) = 0.; 00252 Lopt(2,1) = -0.375 + (5*s)/2.; 00253 Lopt(2,2) = 0.; 00254 Lopt(3,0) = 0.; 00255 Lopt(3,1) = 0.; 00256 Lopt(3,2) = -0.375 + (5*t)/2.; 00257 Lopt(4,0) = (-35*(-1 + 2*r + s + t))/12.; 00258 Lopt(4,1) = (-4 - 35*r + 5*s + 10*t)/12.; 00259 Lopt(4,2) = (-4 - 35*r + 10*s + 5*t)/12.; 00260 Lopt(5,0) = (-1 + 5*r + 40*s - 5*t)/12.; 00261 Lopt(5,1) = (-1 + 40*r + 5*s - 5*t)/12.; 00262 Lopt(5,2) = (-5*(-1 + r + s + 2*t))/12.; 00263 Lopt(6,0) = (-4 + 5*r - 35*s + 10*t)/12.; 00264 Lopt(6,1) = (-35*(-1 + r + 2*s + t))/12.; 00265 Lopt(6,2) = (-4 + 10*r - 35*s + 5*t)/12.; 00266 Lopt(7,0) = (-4 + 5*r + 10*s - 35*t)/12.; 00267 Lopt(7,1) = (-4 + 10*r + 5*s - 35*t)/12.; 00268 Lopt(7,2) = (-35*(-1 + r + s + 2*t))/12.; 00269 Lopt(8,0) = (-1 + 5*r - 5*s + 40*t)/12.; 00270 Lopt(8,1) = (-5*(-1 + r + 2*s + t))/12.; 00271 Lopt(8,2) = (-1 + 40*r - 5*s + 5*t)/12.; 00272 Lopt(9,0) = (-5*(-1 + 2*r + s + t))/12.; 00273 Lopt(9,1) = (-1 - 5*r + 5*s + 40*t)/12.; 00274 Lopt(9,2) = (-1 - 5*r + 40*s + 5*t)/12.; 00275 00276 for (int node=0; node < 10; ++node) { 00277 for (int dim=0; dim < 3; ++dim) { 00278 outputValues(node,pt,dim) = Lopt(node,dim); 00279 } 00280 } 00281 } 00282 } 00283 break; 00284 00285 case OPERATOR_CURL: 00286 TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_CURL), std::invalid_argument, 00287 ">>> ERROR (Basis_HGRAD_TET_COMP12_FEM): CURL is invalid operator for rank-0 (scalar) functions in 3D"); 00288 break; 00289 00290 case OPERATOR_DIV: 00291 TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument, 00292 ">>> ERROR (Basis_HGRAD_TET_COMP12_FEM): DIV is invalid operator for rank-0 (scalar) functions in 3D"); 00293 break; 00294 00295 case OPERATOR_D2: 00296 case OPERATOR_D3: 00297 case OPERATOR_D4: 00298 case OPERATOR_D5: 00299 case OPERATOR_D6: 00300 case OPERATOR_D7: 00301 case OPERATOR_D8: 00302 case OPERATOR_D9: 00303 case OPERATOR_D10: 00304 { 00305 // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, DkCardinality) 00306 int DkCardinality = Intrepid::getDkCardinality(operatorType, 00307 this -> basisCellTopology_.getDimension() ); 00308 for(int dofOrd = 0; dofOrd < this -> basisCardinality_; dofOrd++) { 00309 for (int i0 = 0; i0 < dim0; i0++) { 00310 for(int dkOrd = 0; dkOrd < DkCardinality; dkOrd++){ 00311 outputValues(dofOrd, i0, dkOrd) = 0.0; 00312 } 00313 } 00314 } 00315 } 00316 break; 00317 default: 00318 TEUCHOS_TEST_FOR_EXCEPTION( !( Intrepid::isValidOperator(operatorType) ), std::invalid_argument, 00319 ">>> ERROR (Basis_HGRAD_TET_COMP12_FEM): Invalid operator type"); 00320 } 00321 } 00322 00323 00324 00325 template<class Scalar, class ArrayScalar> 00326 void Basis_HGRAD_TET_COMP12_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar& outputValues, 00327 const ArrayScalar & inputPoints, 00328 const ArrayScalar & cellVertices, 00329 const EOperator operatorType) const { 00330 TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error, 00331 ">>> ERROR (Basis_HGRAD_TET_COMP12_FEM): FEM Basis calling an FVD member function"); 00332 00333 } 00334 00335 template<class Scalar, class ArrayScalar> 00336 void Basis_HGRAD_TET_COMP12_FEM<Scalar, ArrayScalar>::getDofCoords(ArrayScalar & DofCoords) const { 00337 #ifdef HAVE_INTREPID_DEBUG 00338 // Verify rank of output array. 00339 TEUCHOS_TEST_FOR_EXCEPTION( !(DofCoords.rank() == 2), std::invalid_argument, 00340 ">>> ERROR: (Intrepid::Basis_HGRAD_TET_COMP12_FEM::getDofCoords) rank = 2 required for DofCoords array"); 00341 // Verify 0th dimension of output array. 00342 TEUCHOS_TEST_FOR_EXCEPTION( !( DofCoords.dimension(0) == this -> basisCardinality_ ), std::invalid_argument, 00343 ">>> ERROR: (Intrepid::Basis_HGRAD_TET_COMP12_FEM::getDofCoords) mismatch in number of DoF and 0th dimension of DofCoords array"); 00344 // Verify 1st dimension of output array. 00345 TEUCHOS_TEST_FOR_EXCEPTION( !( DofCoords.dimension(1) == (int)(this -> basisCellTopology_.getDimension()) ), std::invalid_argument, 00346 ">>> ERROR: (Intrepid::Basis_HGRAD_TET_COMP12_FEM::getDofCoords) incorrect reference cell (1st) dimension in DofCoords array"); 00347 #endif 00348 00349 DofCoords(0,0) = 0.0; DofCoords(0,1) = 0.0; DofCoords(0,2) = 0.0; 00350 DofCoords(1,0) = 1.0; DofCoords(1,1) = 0.0; DofCoords(1,2) = 0.0; 00351 DofCoords(2,0) = 0.0; DofCoords(2,1) = 1.0; DofCoords(2,2) = 0.0; 00352 DofCoords(3,0) = 0.0; DofCoords(3,1) = 0.0; DofCoords(3,2) = 1.0; 00353 DofCoords(4,0) = 0.5; DofCoords(4,1) = 0.0; DofCoords(4,2) = 0.0; 00354 DofCoords(5,0) = 0.5; DofCoords(5,1) = 0.5; DofCoords(5,2) = 0.0; 00355 DofCoords(6,0) = 0.0; DofCoords(6,1) = 0.5; DofCoords(6,2) = 0.0; 00356 DofCoords(7,0) = 0.0; DofCoords(7,1) = 0.0; DofCoords(7,2) = 0.5; 00357 DofCoords(8,0) = 0.5; DofCoords(8,1) = 0.0; DofCoords(8,2) = 0.5; 00358 DofCoords(9,0) = 0.0; DofCoords(9,1) = 0.5; DofCoords(9,2) = 0.5; 00359 } 00360 00361 template<class Scalar, class ArrayScalar> 00362 Teuchos::Array<int> 00363 Basis_HGRAD_TET_COMP12_FEM<Scalar, ArrayScalar>::getLocalSubTetrahedra(Scalar x, Scalar y, Scalar z) const 00364 { 00365 00366 Teuchos::Array<int> subTets; 00367 int count(0); 00368 00369 // local coords 00370 Scalar xyz = x + y + z; 00371 Scalar xy = x + y; 00372 Scalar xz = x + z; 00373 Scalar yz = y + z; 00374 00375 // cycle through each subdomain and push back if the point lies within 00376 00377 // 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 00378 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) ) { 00379 count++; 00380 subTets.push_back(0); 00381 } 00382 00383 // 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 00384 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) ) { 00385 count++; 00386 subTets.push_back(1); 00387 } 00388 00389 // 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 00390 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) ) { 00391 count++; 00392 subTets.push_back(2); 00393 } 00394 00395 // 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 00396 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) ) { 00397 count++; 00398 subTets.push_back(3); 00399 } 00400 00401 // 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 00402 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) ) { 00403 count++; 00404 subTets.push_back(4); 00405 } 00406 00407 // 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 00408 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) ) { 00409 count++; 00410 subTets.push_back(5); 00411 } 00412 00413 // 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 00414 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) ) { 00415 count++; 00416 subTets.push_back(6); 00417 } 00418 00419 // 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 00420 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) ) { 00421 count++; 00422 subTets.push_back(7); 00423 } 00424 00425 // 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 00426 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) ) { 00427 count++; 00428 subTets.push_back(8); 00429 } 00430 00431 // 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 00432 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) ) { 00433 count++; 00434 subTets.push_back(9); 00435 } 00436 00437 // 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 00438 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) ) { 00439 count++; 00440 subTets.push_back(10); 00441 } 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 count++; 00446 subTets.push_back(11); 00447 } 00448 00449 // if the point doesn't lie in the parent domain return -1 00450 if (count == 0) { 00451 subTets.push_back(-1); 00452 } 00453 00454 return subTets; 00455 } 00456 00457 00458 template<class Scalar, class ArrayScalar> 00459 Intrepid::FieldContainer<Scalar> 00460 Basis_HGRAD_TET_COMP12_FEM<Scalar, ArrayScalar>::getWeights(const ArrayScalar & inPts) const 00461 { 00462 int numPoints = inPts.dimension(0); 00463 Intrepid::FieldContainer<Scalar> w(numPoints, 12); 00464 w.initialize(0.0); 00465 Teuchos::Array< Teuchos::Array<int> > pt_tets; 00466 00467 for (int pt = 0; pt < numPoints; ++pt) 00468 pt_tets.push_back(getLocalSubTetrahedra(inPts(pt,0), inPts(pt,1), inPts(pt,2))); 00469 00470 Teuchos::Array<int> flat; 00471 FieldContainer<int> count(12); 00472 00473 for (int pt = 0; pt < numPoints; ++pt) 00474 for (int i = 0; i < pt_tets[pt].size(); ++i) 00475 flat.push_back(pt_tets[pt][i]); 00476 00477 for (int i = 0; i < flat.size(); ++i) 00478 count(flat[i])++; 00479 00480 for (int pt = 0; pt < numPoints; ++pt) 00481 for (int i = 0; i < pt_tets[pt].size(); ++i) 00482 w(pt, pt_tets[pt][i]) = 1.0/count(pt_tets[pt][i]); 00483 00484 return w; 00485 } 00486 00487 00488 00489 template<class Scalar, class ArrayScalar> 00490 Intrepid::FieldContainer<Scalar> 00491 Basis_HGRAD_TET_COMP12_FEM<Scalar, ArrayScalar>::getBarycentricCoords(const ArrayScalar & inPts) const 00492 { 00493 int numPoints = inPts.dimension(0); 00494 Intrepid::FieldContainer<Scalar> lambda(numPoints, 4); 00495 00496 for (int pt = 0; pt < numPoints; ++pt) 00497 { 00498 lambda(pt,0) = 1. - inPts(pt,0) - inPts(pt,1) - inPts(pt,2); 00499 lambda(pt,1) = inPts(pt,0); 00500 lambda(pt,2) = inPts(pt,1); 00501 lambda(pt,3) = inPts(pt,2); 00502 } 00503 00504 return lambda; 00505 } 00506 00507 template<class Scalar, class ArrayScalar> 00508 Scalar 00509 Basis_HGRAD_TET_COMP12_FEM<Scalar, ArrayScalar>::det44(const Intrepid::FieldContainer<Scalar> a) const 00510 { 00511 Scalar det = a(0,3) * a(1,2) * a(2,1) * a(3,0) 00512 - a(0,2) * a(1,3) * a(2,1) * a(3,0) 00513 - a(0,3) * a(1,1) * a(2,2) * a(3,0) 00514 + a(0,1) * a(1,3) * a(2,2) * a(3,0) 00515 + a(0,2) * a(1,1) * a(2,3) * a(3,0) 00516 - a(0,1) * a(1,2) * a(2,3) * a(3,0) 00517 - a(0,3) * a(1,2) * a(2,0) * a(3,1) 00518 + a(0,2) * a(1,3) * a(2,0) * a(3,1) 00519 + a(0,3) * a(1,0) * a(2,2) * a(3,1) 00520 - a(0,0) * a(1,3) * a(2,2) * a(3,1) 00521 - a(0,2) * a(1,0) * a(2,3) * a(3,1) 00522 + a(0,0) * a(1,2) * a(2,3) * a(3,1) 00523 + a(0,3) * a(1,1) * a(2,0) * a(3,2) 00524 - a(0,1) * a(1,3) * a(2,0) * a(3,2) 00525 - a(0,3) * a(1,0) * a(2,1) * a(3,2) 00526 + a(0,0) * a(1,3) * a(2,1) * a(3,2) 00527 + a(0,1) * a(1,0) * a(2,3) * a(3,2) 00528 - a(0,0) * a(1,1) * a(2,3) * a(3,2) 00529 - a(0,2) * a(1,1) * a(2,0) * a(3,3) 00530 + a(0,1) * a(1,2) * a(2,0) * a(3,3) 00531 + a(0,2) * a(1,0) * a(2,1) * a(3,3) 00532 - a(0,0) * a(1,2) * a(2,1) * a(3,3) 00533 - a(0,1) * a(1,0) * a(2,2) * a(3,3) 00534 + a(0,0) * a(1,1) * a(2,2) * a(3,3); 00535 00536 return det; 00537 } 00538 00539 template<class Scalar, class ArrayScalar> 00540 Intrepid::FieldContainer<Scalar> 00541 Basis_HGRAD_TET_COMP12_FEM<Scalar, ArrayScalar>::inverse44(const Intrepid::FieldContainer<Scalar> a) const 00542 { 00543 Intrepid::FieldContainer<Scalar> ai(4,4); 00544 Scalar xj = det44(a); 00545 00546 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)); 00547 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)); 00548 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)); 00549 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)); 00550 00551 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)); 00552 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)); 00553 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)); 00554 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)); 00555 00556 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)); 00557 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)); 00558 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)); 00559 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)); 00560 00561 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)); 00562 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)); 00563 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)); 00564 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)); 00565 00566 return ai; 00567 } 00568 00569 template<class Scalar, class ArrayScalar> 00570 Intrepid::FieldContainer<Scalar> 00571 Basis_HGRAD_TET_COMP12_FEM<Scalar, ArrayScalar>::getSubTetGrads() const 00572 { 00573 Intrepid::FieldContainer<Scalar> dx(3,11,12); 00574 dx.initialize(0.0); 00575 // fill in dx 00576 dx(0,0,0) = -2.0; 00577 dx(0,1,1) = 2.0; 00578 dx(0,4,0) = 2.0; 00579 dx(0,4,1) = -2.0; 00580 dx(0,5,2) = 2.0; 00581 dx(0,5,4) = 2.0; 00582 dx(0,5,5) = 2.0; 00583 dx(0,5,8) = 2.0; 00584 dx(0,5,9) = 2.0; 00585 dx(0,6,2) = -2.0; 00586 dx(0,6,8) = -2.0; 00587 dx(0,6,9) = -2.0; 00588 dx(0,6,10) = -2.0; 00589 dx(0,6,11) = -2.0; 00590 dx(0,7,3) = -2.0; 00591 dx(0,7,6) = -2.0; 00592 dx(0,7,7) = -2.0; 00593 dx(0,7,10) = -2.0; 00594 dx(0,7,11) = -2.0; 00595 dx(0,8,3) = 2.0; 00596 dx(0,8,4) = 2.0; 00597 dx(0,8,5) = 2.0; 00598 dx(0,8,6) = 2.0; 00599 dx(0,8,7) = 2.0; 00600 dx(0,10,4) = -4.0; 00601 dx(0,10,5) = -4.0; 00602 dx(0,10,10) = 4.0; 00603 dx(0,10,11) = 4.0; 00604 00605 // fill in dy 00606 dx(1,0,0) = -2.0; 00607 dx(1,2,2) = 2.0; 00608 dx(1,4,1) = -2.0; 00609 dx(1,4,4) = -2.0; 00610 dx(1,4,7) = -2.0; 00611 dx(1,4,8) = -2.0; 00612 dx(1,4,11) = -2.0; 00613 dx(1,5,1) = 2.0; 00614 dx(1,5,4) = 2.0; 00615 dx(1,5,5) = 2.0; 00616 dx(1,5,8) = 2.0; 00617 dx(1,5,9) = 2.0; 00618 dx(1,6,0) = 2.0; 00619 dx(1,6,2) = -2.0; 00620 dx(1,7,3) = -2.0; 00621 dx(1,7,6) = -2.0; 00622 dx(1,7,7) = -2.0; 00623 dx(1,7,10) = -2.0; 00624 dx(1,7,11) = -2.0; 00625 dx(1,9,3) = 2.0; 00626 dx(1,9,5) = 2.0; 00627 dx(1,9,6) = 2.0; 00628 dx(1,9,9) = 2.0; 00629 dx(1,9,10) = 2.0; 00630 dx(1,10,5) = -4.0; 00631 dx(1,10,7) = 4.0; 00632 dx(1,10,9) = -4.0; 00633 dx(1,10,11) = 4.0; 00634 00635 // fill in dz 00636 dx(2,0,0) = -2.0; 00637 dx(2,3,3) = 2.0; 00638 dx(2,4,1) = -2.0; 00639 dx(2,4,4) = -2.0; 00640 dx(2,4,7) = -2.0; 00641 dx(2,4,8) = -2.0; 00642 dx(2,4,11) = -2.0; 00643 dx(2,6,2) = -2.0; 00644 dx(2,6,8) = -2.0; 00645 dx(2,6,9) = -2.0; 00646 dx(2,6,10) = -2.0; 00647 dx(2,6,11) = -2.0; 00648 dx(2,7,0) = 2.0; 00649 dx(2,7,3) = -2.0; 00650 dx(2,8,1) = 2.0; 00651 dx(2,8,4) = 2.0; 00652 dx(2,8,5) = 2.0; 00653 dx(2,8,6) = 2.0; 00654 dx(2,8,7) = 2.0; 00655 dx(2,9,2) = 2.0; 00656 dx(2,9,5) = 2.0; 00657 dx(2,9,6) = 2.0; 00658 dx(2,9,9) = 2.0; 00659 dx(2,9,10) = 2.0; 00660 dx(2,10,5) = -4.0; 00661 dx(2,10,6) = -4.0; 00662 dx(2,10,8) = 4.0; 00663 dx(2,10,11) = 4.0; 00664 00665 return dx; 00666 } 00667 00668 template<class Scalar, class ArrayScalar> 00669 Intrepid::FieldContainer<Scalar> 00670 Basis_HGRAD_TET_COMP12_FEM<Scalar, ArrayScalar>::getSubTetDetF() const 00671 { 00672 Intrepid::FieldContainer<Scalar> xJ(12); 00673 // set sub-elem jacobians 00674 xJ(0) = 1./48.; xJ(1) = 1./48.; xJ(2) = 1./48.; xJ(3) = 1./48.; 00675 xJ(4) = 1./96.; xJ(5) = 1./96.; xJ(6) = 1./96.; xJ(7) = 1./96.; 00676 xJ(8) = 1./96.; xJ(9) = 1./96.; xJ(10) = 1./96.; xJ(11) = 1./96.; 00677 00678 return xJ; 00679 } 00680 00681 00682 }// namespace Intrepid 00683 #endif
1.7.6.1