Intrepid
http://trilinos.sandia.gov/packages/docs/r10.10/packages/intrepid/src/Discretization/Basis/Intrepid_HGRAD_TRI_Cn_FEM_ORTHDef.hpp
Go to the documentation of this file.
00001 #ifndef INTREPID_HGRAD_TRI_CN_FEM_ORTHDEF_HPP
00002 #define INTREPID_HGRAD_TRI_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_TRI_Cn_FEM_ORTH<Scalar,ArrayScalar>::Basis_HGRAD_TRI_Cn_FEM_ORTH( int degree )
00056 {
00057     this -> basisCardinality_  = (degree+1)*(degree+2)/2;
00058     this -> basisDegree_       = degree;
00059     this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Triangle<3> >() );
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_TRI_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_TRI_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   // add more here and put in appropriate extra case statements below to enable higher derivatives.
00114   void (*tabulators[])(ArrayScalar &, const int, const ArrayScalar &)
00115     = { TabulatorTri<Scalar,ArrayScalar,0>::tabulate ,
00116         TabulatorTri<Scalar,ArrayScalar,1>::tabulate ,
00117         TabulatorTri<Scalar,ArrayScalar,2>::tabulate };
00118 
00119 
00120   switch (operatorType) {
00121   case OPERATOR_VALUE:
00122     tabulators[0]( outputValues , deg , inputPoints );
00123     break;
00124   case OPERATOR_GRAD:
00125     tabulators[1]( outputValues , deg , inputPoints );
00126     break;
00127   case OPERATOR_D1:
00128   case OPERATOR_D2:
00129     // add more case OPEATOR_Dn statements if you've added more items to the
00130     // array above.
00131     tabulators[operatorType-OPERATOR_D1+1]( outputValues , deg , inputPoints );
00132     break;
00133   default:
00134     TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument,
00135                         ">>> ERROR (Basis_HGRAD_TRI_Cn_FEM_ORTH): invalid or unsupported operator" );
00136 
00137   }
00138 
00139   return;
00140 }
00141 
00142 
00143 
00144 template<class Scalar, class ArrayScalar>
00145 void Basis_HGRAD_TRI_Cn_FEM_ORTH<Scalar, ArrayScalar>::getValues(ArrayScalar&           outputValues,
00146                                                                  const ArrayScalar &    inputPoints,
00147                                                                  const ArrayScalar &    cellVertices,
00148                                                                  const EOperator        operatorType) const {
00149   TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
00150                       ">>> ERROR (Basis_HGRAD_TRI_Cn_FEM): FEM Basis calling an FVD member function");
00151 }
00152 
00153 
00154 
00155 template<typename Scalar, typename ArrayScalar>
00156 void TabulatorTri<Scalar,ArrayScalar,0>::tabulate(ArrayScalar &outputValues ,
00157                                                   const int deg ,
00158                                                   const ArrayScalar &z )
00159 {
00160   const int np = z.dimension( 0 );
00161   
00162   // each point needs to be transformed from Pavel's element
00163   // z(i,0) --> (2.0 * z(i,0) - 1.0)
00164   // z(i,1) --> (2.0 * z(i,1) - 1.0)
00165   
00166   // set up constant term
00167   int idx_cur = TabulatorTri<Scalar,ArrayScalar,0>::idx(0,0);
00168   int idx_curp1,idx_curm1;
00169   
00170   // set D^{0,0} = 1.0
00171   for (int i=0;i<np;i++) {
00172     outputValues(idx_cur,i) = Scalar( 1.0 ) + z(i,0) - z(i,0) + z(i,1) - z(i,1);
00173   }
00174   
00175 
00176   if (deg > 0) {
00177     Teuchos::Array<Scalar> f1(np),f2(np),f3(np);
00178     
00179     for (int i=0;i<np;i++) {
00180       f1[i] = 0.5 * (1.0+2.0*(2.0*z(i,0)-1.0)+(2.0*z(i,1)-1.0));
00181       f2[i] = 0.5 * (1.0-(2.0*z(i,1)-1.0));
00182       f3[i] = f2[i] * f2[i];
00183     }
00184     
00185     // set D^{1,0} = f1
00186     idx_cur = TabulatorTri<Scalar,ArrayScalar,0>::idx(1,0);
00187     for (int i=0;i<np;i++) {
00188       outputValues(idx_cur,i) = f1[i];
00189     }
00190     
00191     // recurrence in p
00192     for (int p=1;p<deg;p++) {
00193       idx_cur = TabulatorTri<Scalar,ArrayScalar,0>::idx(p,0);
00194       idx_curp1 = TabulatorTri<Scalar,ArrayScalar,0>::idx(p+1,0);
00195       idx_curm1 = TabulatorTri<Scalar,ArrayScalar,0>::idx(p-1,0);
00196       Scalar a = (2.0*p+1.0)/(1.0+p);
00197       Scalar b = p / (p+1.0);
00198       
00199       for (int i=0;i<np;i++) {
00200         outputValues(idx_curp1,i) = a * f1[i] * outputValues(idx_cur,i)
00201           - b * f3[i] * outputValues(idx_curm1,i);
00202       }
00203     }
00204     
00205     // D^{p,1}
00206     for (int p=0;p<deg;p++) {
00207       int idxp0 = TabulatorTri<Scalar,ArrayScalar,0>::idx(p,0);
00208       int idxp1 = TabulatorTri<Scalar,ArrayScalar,0>::idx(p,1);
00209       for (int i=0;i<np;i++) {
00210         outputValues(idxp1,i) = outputValues(idxp0,i)
00211           *0.5*(1.0+2.0*p+(3.0+2.0*p)*(2.0*z(i,1)-1.0));
00212       }
00213     }
00214     
00215     
00216     // recurrence in q
00217     for (int p=0;p<deg-1;p++) {
00218       for (int q=1;q<deg-p;q++) {
00219         int idxpqp1=TabulatorTri<Scalar,ArrayScalar,0>::idx(p,q+1);
00220         int idxpq=TabulatorTri<Scalar,ArrayScalar,0>::idx(p,q);
00221         int idxpqm1=TabulatorTri<Scalar,ArrayScalar,0>::idx(p,q-1);
00222         Scalar a,b,c;
00223         TabulatorTri<Scalar,ArrayScalar,0>::jrc((Scalar)(2*p+1),(Scalar)0,q,a,b,c);
00224         for (int i=0;i<np;i++) {
00225           outputValues(idxpqp1,i)
00226             = (a*(2.0*z(i,1)-1.0)+b)*outputValues(idxpq,i)
00227             - c*outputValues(idxpqm1,i);
00228         }
00229       }
00230     }
00231   }    
00232   
00233   // orthogonalize
00234   for (int p=0;p<=deg;p++) {
00235     for (int q=0;q<=deg-p;q++) {
00236       for (int i=0;i<np;i++) {
00237         outputValues(TabulatorTri<Scalar,ArrayScalar,0>::idx(p,q),i) *= sqrt( (p+0.5)*(p+q+1.0));
00238       }
00239     }
00240   }
00241   
00242   return;
00243 }
00244  
00245 
00246  
00247 template<typename Scalar, typename ArrayScalar>
00248 void TabulatorTri<Scalar,ArrayScalar,1>::tabulate(ArrayScalar &outputValues ,
00249                                                   const int deg ,
00250                                                   const ArrayScalar &z ) 
00251 {
00252   const int np = z.dimension(0);
00253   const int card = outputValues.dimension(0);
00254   FieldContainer<Sacado::Fad::DFad<Scalar> > dZ( z.dimension(0) , z.dimension(1) );
00255   for (int i=0;i<np;i++) {
00256     for (int j=0;j<2;j++) {
00257       dZ(i,j) = Sacado::Fad::DFad<Scalar>( z(i,j) );
00258       dZ(i,j).diff(j,2);
00259     }
00260   }
00261   FieldContainer<Sacado::Fad::DFad<Scalar> > dResult(card,np);
00262 
00263   TabulatorTri<Sacado::Fad::DFad<Scalar>,FieldContainer<Sacado::Fad::DFad<Scalar> >,0>::tabulate( dResult ,
00264                                                                                                   deg ,
00265                                                                                                   dZ );
00266 
00267   for (int i=0;i<card;i++) {
00268     for (int j=0;j<np;j++) {
00269       for (int k=0;k<2;k++) {
00270         outputValues(i,j,k) = dResult(i,j).dx(k);
00271       }
00272     }
00273   }
00274 
00275   return;
00276 
00277 }
00278 
00279 
00280 
00281 template<typename Scalar, typename ArrayScalar, unsigned derivOrder>
00282 void TabulatorTri<Scalar,ArrayScalar,derivOrder>::tabulate( ArrayScalar &outputValues ,
00283                                                             const int deg ,
00284                                                             const ArrayScalar &z ) 
00285 {
00286   const int np = z.dimension(0);
00287   const int card = outputValues.dimension(0);
00288   FieldContainer<Sacado::Fad::DFad<Scalar> > dZ( z.dimension(0) , z.dimension(1) );
00289   for (int i=0;i<np;i++) {
00290     for (int j=0;j<2;j++) {
00291       dZ(i,j) = Sacado::Fad::DFad<Scalar>( z(i,j) );
00292       dZ(i,j).diff(j,2);
00293     }
00294   }
00295   FieldContainer<Sacado::Fad::DFad<Scalar> > dResult(card,np,derivOrder+1);
00296 
00297   TabulatorTri<Sacado::Fad::DFad<Scalar>,FieldContainer<Sacado::Fad::DFad<Scalar> >,derivOrder-1>::tabulate(dResult ,
00298                                                                                                             deg ,
00299                                                                                                             dZ );
00300 
00301   for (int i=0;i<card;i++) {
00302     for (int j=0;j<np;j++) {
00303       outputValues(i,j,0) = dResult(i,j,0).dx(0);
00304       for (unsigned k=0;k<derivOrder;k++) {
00305         outputValues(i,j,k+1) = dResult(i,j,k).dx(1);
00306       }
00307     }
00308   }
00309 
00310   return;
00311 
00312 
00313 }
00314 
00315 
00316 }// namespace Intrepid
00317 #endif