Intrepid
http://trilinos.sandia.gov/packages/docs/dev/packages/intrepid/src/Discretization/Basis/Intrepid_HGRAD_TET_Cn_FEM_ORTHDef.hpp
Go to the documentation of this file.
00001 #ifndef INTREPID_HGRAD_TET_CN_FEM_ORTHDEF_HPP
00002 #define INTREPID_HGRAD_TET_CN_FEM_ORTHDEF_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_Cn_FEM_ORTH<Scalar,ArrayScalar>::Basis_HGRAD_TET_Cn_FEM_ORTH( int degree )
00056   {
00057     this -> basisCardinality_  = (degree+1)*(degree+2)*(degree+3)/6;
00058     this -> basisDegree_       = degree;
00059     this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<4> >() );
00060     this -> basisType_         = BASIS_FEM_HIERARCHICAL;
00061     this -> basisCoordinates_  = COORDINATES_CARTESIAN;
00062     this -> basisTagsAreSet_   = false;
00063   }
00064   
00065   
00066   
00067   template<class Scalar, class ArrayScalar>
00068   void Basis_HGRAD_TET_Cn_FEM_ORTH<Scalar, ArrayScalar>::initializeTags() {
00069   
00070     // Basis-dependent initializations
00071     int tagSize  = 4;        // size of DoF tag, i.e., number of fields in the tag
00072     int posScDim = 0;        // position in the tag, counting from 0, of the subcell dim 
00073     int posScOrd = 1;        // position in the tag, counting from 0, of the subcell ordinal
00074     int posDfOrd = 2;        // position in the tag, counting from 0, of DoF ordinal relative to the subcell
00075   
00076     // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration 
00077     int *tags = new int[tagSize * this->getCardinality()];
00078     for (int i=0;i<this->getCardinality();i++) {
00079       tags[4*i] = 2;
00080       tags[4*i+1] = 0;
00081       tags[4*i+2] = i;
00082       tags[4*i+3] = this->getCardinality();
00083     }
00084 
00085     // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
00086     Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
00087                                 this -> ordinalToTag_,
00088                                 tags,
00089                                 this -> basisCardinality_,
00090                                 tagSize,
00091                                 posScDim,
00092                                 posScOrd,
00093                                 posDfOrd);
00094   }  
00095   
00096 
00097 
00098   template<class Scalar, class ArrayScalar> 
00099   void Basis_HGRAD_TET_Cn_FEM_ORTH<Scalar, ArrayScalar>::getValues(ArrayScalar &        outputValues,
00100                                                                   const ArrayScalar &  inputPoints,
00101                                                                   const EOperator      operatorType) const {
00102   
00103     // Verify arguments
00104 #ifdef HAVE_INTREPID_DEBUG
00105     Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
00106                                                         inputPoints,
00107                                                         operatorType,
00108                                                         this -> getBaseCellTopology(),
00109                                                         this -> getCardinality() );
00110 #endif
00111     const int deg = this->getDegree();
00112   
00113     switch (operatorType) {
00114     case OPERATOR_VALUE:
00115       {
00116         TabulatorTet<Scalar,ArrayScalar,0>::tabulate( outputValues ,
00117                                                       deg ,
00118                                                       inputPoints );
00119       }
00120       break;
00121     case OPERATOR_GRAD:
00122     case OPERATOR_D1:
00123       {
00124         TabulatorTet<Scalar,ArrayScalar,1>::tabulate( outputValues ,
00125                                                       deg ,
00126                                                       inputPoints );
00127       }
00128       break;
00129     default:
00130       TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument,
00131                           ">>> ERROR (Basis_HGRAD_TET_Cn_FEM_ORTH): invalid or unsupported operator" );
00132     }
00133 
00134     return;
00135   }
00136   
00137   template<class Scalar, class ArrayScalar>
00138   void Basis_HGRAD_TET_Cn_FEM_ORTH<Scalar, ArrayScalar>::getValues(ArrayScalar&           outputValues,
00139                                                                   const ArrayScalar &    inputPoints,
00140                                                                   const ArrayScalar &    cellVertices,
00141                                                                   const EOperator        operatorType) const {
00142     TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
00143                         ">>> ERROR (Basis_HGRAD_TET_Cn_FEM_ORTH): FEM Basis calling an FVD member function");
00144   }
00145 
00146   template<class Scalar, class ArrayScalar>
00147   void TabulatorTet<Scalar,ArrayScalar,0>::tabulate( ArrayScalar &outputValues ,
00148                                                     const int deg ,
00149                                                     const ArrayScalar &z )
00150   {
00151     const int np = z.dimension( 0 );
00152     int idxcur;
00153   
00154     // each point needs to be transformed from Pavel's element
00155     // z(i,0) --> (2.0 * z(i,0) - 1.0)
00156     // z(i,1) --> (2.0 * z(i,1) - 1.0)
00157     // z(i,2) --> (2.0 * z(i,2) - 1.0)
00158   
00159     Teuchos::Array<Scalar> f1(np),f2(np),f3(np),f4(np),f5(np);
00160   
00161     for (int i=0;i<np;i++) {
00162       f1[i] = 0.5 * ( 2.0 + 2.0*(2.0*z(i,0)-1.0) + (2.0*z(i,1)-1.0) + (2.0*z(i,2)-1.0) );
00163       Scalar foo =  0.5 * ( (2.0*z(i,1)-1.0) + (2.0*z(i,2)-1.0) );
00164       f2[i] = foo * foo;
00165       f3[i] = 0.5 * ( 1.0 + 2.0 * (2.0*z(i,1)-1.0) + (2.0*z(i,2)-1.0) );
00166       f4[i] = 0.5 * ( 1.0 - (2.0*z(i,2)-1.0) );
00167       f5[i] = f4[i] * f4[i];
00168     }
00169 
00170     // constant term
00171     idxcur = TabulatorTet<Scalar,ArrayScalar,0>::idx(0,0,0);
00172     for (int i=0;i<np;i++) {
00173       outputValues(idxcur,i) = 1.0 + z(i,0) - z(i,0) + z(i,1) - z(i,1) + z(i,2) - z(i,2);
00174     }
00175   
00176     if (deg > 0) {
00177 
00178       // D^{1,0,0}
00179       idxcur = TabulatorTet<Scalar,ArrayScalar,0>::idx(1,0,0);
00180       for (int i=0;i<np;i++) {
00181         outputValues(idxcur,i) = f1[i];
00182       }
00183   
00184       // p recurrence
00185       for (int p=1;p<deg;p++) {
00186         Scalar a1 = (2.0 * p + 1.0) / ( p + 1.0);
00187         Scalar a2 = p / ( p + 1.0 );
00188         int idxp = TabulatorTet<Scalar,ArrayScalar,0>::idx(p,0,0);
00189         int idxpp1 = TabulatorTet<Scalar,ArrayScalar,0>::idx(p+1,0,0);
00190         int idxpm1 = TabulatorTet<Scalar,ArrayScalar,0>::idx(p-1,0,0);
00191         for (int i=0;i<np;i++) {
00192           outputValues(idxpp1,i) = a1 * f1[i] * outputValues(idxp,i) - a2 * f2[i] * outputValues(idxpm1,i);
00193         }
00194       }
00195       // q = 1
00196       for (int p=0;p<deg;p++) {
00197         int idx0 = TabulatorTet<Scalar,ArrayScalar,0>::idx(p,0,0);
00198         int idx1 = TabulatorTet<Scalar,ArrayScalar,0>::idx(p,1,0);
00199         for (int i=0;i<np;i++) {
00200           outputValues(idx1,i) = outputValues(idx0,i) * ( p * ( 1.0 + (2.0*z(i,1)-1.0) ) +
00201                                                           0.5 * ( 2.0 + 3.0 * (2.0*z(i,1)-1.0) + (2.0*z(i,2)-1.0) ) );
00202         }
00203       }
00204   
00205       // q recurrence
00206       for (int p=0;p<deg-1;p++) {
00207         for (int q=1;q<deg-p;q++) {
00208           Scalar aq,bq,cq;
00209 
00210           TabulatorTet<Scalar,ArrayScalar,0>::jrc(2.0*p+1.0 ,0 ,q, aq, bq, cq);
00211           int idxpqp1 = TabulatorTet<Scalar,ArrayScalar,0>::idx(p,q+1,0);
00212           int idxpq = TabulatorTet<Scalar,ArrayScalar,0>::idx(p,q,0);
00213           int idxpqm1 = TabulatorTet<Scalar,ArrayScalar,0>::idx(p,q-1,0);
00214           for (int i=0;i<np;i++) {
00215             outputValues(idxpqp1,i) = ( aq * f3[i] + bq * f4[i] ) * outputValues(idxpq,i) 
00216               - ( cq * f5[i] ) * outputValues(idxpqm1,i);
00217           }
00218         }
00219       }
00220   
00221       // r = 1
00222       for (int p=0;p<deg;p++) {
00223         for (int q=0;q<deg-p;q++) {
00224           int idxpq1 = TabulatorTet<Scalar,ArrayScalar,0>::idx(p,q,1);
00225           int idxpq0 = TabulatorTet<Scalar,ArrayScalar,0>::idx(p,q,0);
00226           for (int i=0;i<np;i++) {
00227             outputValues(idxpq1,i) = outputValues(idxpq0,i) * ( 1.0 + p + q + ( 2.0 + q + 
00228                                                                                 p ) * (2.0*z(i,2)-1.0) );
00229           }
00230         }
00231       }
00232       // general r recurrence
00233       for (int p=0;p<deg-1;p++) {
00234         for (int q=0;q<deg-p-1;q++) {
00235           for (int r=1;r<deg-p-q;r++) {
00236             Scalar ar,br,cr;
00237             int idxpqrp1 = TabulatorTet<Scalar,ArrayScalar,0>::idx(p,q,r+1);
00238             int idxpqr = TabulatorTet<Scalar,ArrayScalar,0>::idx(p,q,r);
00239             int idxpqrm1 = TabulatorTet<Scalar,ArrayScalar,0>::idx(p,q,r-1);
00240             jrc(2.0*p+2.0*q+2.0, 0.0, r, ar, br, cr);
00241             for (int i=0;i<np;i++) {
00242               outputValues(idxpqrp1,i) = (ar * (2.0*z(i,2)-1.0) + br) * outputValues( idxpqr , i ) - cr * outputValues(idxpqrm1,i);
00243             }
00244           }
00245         }
00246       }
00247 
00248     }  
00249     // normalize
00250     for (int p=0;p<=deg;p++) {
00251       for (int q=0;q<=deg-p;q++) {
00252         for (int r=0;r<=deg-p-q;r++) {
00253           int idxcur = TabulatorTet<Scalar,ArrayScalar,0>::idx(p,q,r);
00254           Scalar scal = sqrt( (p+0.5)*(p+q+1.0)*(p+q+r+1.5) );
00255           for (int i=0;i<np;i++) {
00256             outputValues(idxcur,i) *= scal;
00257           }
00258         }
00259       }
00260     }
00261   
00262     return;
00263   
00264   }
00265 
00266 
00267   template<typename Scalar, typename ArrayScalar>
00268   void TabulatorTet<Scalar,ArrayScalar,1>::tabulate( ArrayScalar &outputValues ,
00269                                                     const int deg ,
00270                                                     const ArrayScalar &z ) 
00271   {
00272     const int np = z.dimension(0);
00273     const int card = outputValues.dimension(0);
00274     FieldContainer<Sacado::Fad::DFad<Scalar> > dZ( z.dimension(0) , z.dimension(1) );
00275     for (int i=0;i<np;i++) {
00276       for (int j=0;j<3;j++) {
00277         dZ(i,j) = Sacado::Fad::DFad<Scalar>( z(i,j) );
00278         dZ(i,j).diff(j,3);
00279       }
00280     }
00281     FieldContainer<Sacado::Fad::DFad<Scalar> > dResult(card,np);
00282 
00283     TabulatorTet<Sacado::Fad::DFad<Scalar>,FieldContainer<Sacado::Fad::DFad<Scalar> >,0>::tabulate( dResult ,
00284                                                                                                     deg ,
00285                                                                                                     dZ );
00286 
00287     for (int i=0;i<card;i++) {
00288       for (int j=0;j<np;j++) {
00289         for (int k=0;k<3;k++) {
00290           outputValues(i,j,k) = dResult(i,j).dx(k);
00291         }
00292       }
00293     }
00294 
00295     return;
00296 
00297   }
00298 
00299 
00300 }// namespace Intrepid
00301 
00302 
00303 #endif