Intrepid
http://trilinos.sandia.gov/packages/docs/dev/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 
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