Intrepid
http://trilinos.sandia.gov/packages/docs/r10.8/packages/intrepid/src/Discretization/Basis/Intrepid_HGRAD_TET_COMP12_FEMDef.hpp
Go to the documentation of this file.
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       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       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       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   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   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   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   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<Scalar> 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