Intrepid
http://trilinos.sandia.gov/packages/docs/r10.10/packages/intrepid/src/Discretization/Basis/Intrepid_HGRAD_QUAD_Cn_FEMDef.hpp
Go to the documentation of this file.
00001 #ifndef INTREPID_HGRAD_QUAD_CN_FEMDEF_HPP
00002 #define INTREPID_HGRAD_QUAD_CN_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 
00051 using Teuchos::Array;
00052 using Teuchos::RCP;
00053 
00054 namespace Intrepid {
00055   template<class Scalar, class ArrayScalar>
00056   Basis_HGRAD_QUAD_Cn_FEM<Scalar,ArrayScalar>::Basis_HGRAD_QUAD_Cn_FEM( const int orderx , const int ordery,
00057                                                                         const ArrayScalar &pts_x ,
00058                                                                         const ArrayScalar &pts_y ): 
00059     ptsx_( pts_x.dimension(0) , 1 ) ,
00060     ptsy_( pts_y.dimension(0) , 1 )
00061   {
00062     Array<Array<RCP<Basis<Scalar,ArrayScalar> > > > bases(1);
00063     bases[0].resize(2);
00064     bases[0][0] = Teuchos::rcp( new Basis_HGRAD_LINE_Cn_FEM<Scalar,ArrayScalar>( orderx , pts_x ) );
00065     bases[0][1] = Teuchos::rcp( new Basis_HGRAD_LINE_Cn_FEM<Scalar,ArrayScalar>( ordery , pts_y ) );
00066     this->setBases( bases );
00067 
00068     this->basisCardinality_ = (orderx+1)*(ordery+1);
00069     if (orderx > ordery) {
00070       this->basisDegree_ = orderx;
00071     }
00072     else {
00073       this->basisDegree_ = ordery;
00074     }
00075     this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
00076     this -> basisType_         = BASIS_FEM_FIAT;
00077     this -> basisCoordinates_  = COORDINATES_CARTESIAN;
00078     this -> basisTagsAreSet_   = false;
00079 
00080     for (int i=0;i<pts_x.dimension(0);i++)
00081       {
00082         ptsx_(i,0) = pts_x(i,0);
00083       }
00084 
00085     for (int i=0;i<pts_y.dimension(0);i++)
00086       {
00087         ptsy_(i,0) = pts_y(i,0);
00088       }
00089 
00090   }
00091 
00092   template<class Scalar, class ArrayScalar>
00093   Basis_HGRAD_QUAD_Cn_FEM<Scalar,ArrayScalar>::Basis_HGRAD_QUAD_Cn_FEM( const int order,
00094                                                                         const EPointType &pointType ):
00095     ptsx_( order+1 , 1 ) ,
00096     ptsy_( order+1 , 1 )
00097   {
00098     Array<Array<RCP<Basis<Scalar,ArrayScalar> > > > bases(1);
00099     bases[0].resize(2);
00100     bases[0][0] = Teuchos::rcp( new Basis_HGRAD_LINE_Cn_FEM<Scalar,ArrayScalar>( order , pointType ) );
00101     bases[0][1] = Teuchos::rcp( new Basis_HGRAD_LINE_Cn_FEM<Scalar,ArrayScalar>( order , pointType ) );
00102     this->setBases( bases );
00103 
00104     this->basisCardinality_ = (order+1)*(order+1);
00105     this->basisDegree_ = order;
00106     this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
00107     this -> basisType_         = BASIS_FEM_FIAT;
00108     this -> basisCoordinates_  = COORDINATES_CARTESIAN;
00109     this -> basisTagsAreSet_   = false;
00110 
00111     // fill up the pt arrays with calls to the lattice
00112     EPointType pt = (pointType==POINTTYPE_EQUISPACED)?pointType:POINTTYPE_WARPBLEND;
00113     PointTools::getLattice<Scalar,FieldContainer<Scalar> >( ptsx_ ,
00114                                                             shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >()) ,
00115                                                             order ,
00116                                                             0 ,
00117                                                             pt );
00118           
00119     for (int i=0;i<order+1;i++)
00120       {
00121         ptsy_(i,0) = ptsx_(i,0);
00122       }
00123   }
00124 
00125   template<class Scalar, class ArrayScalar>
00126   void Basis_HGRAD_QUAD_Cn_FEM<Scalar,ArrayScalar>::initializeTags()
00127   {
00128     // Basis-dependent initializations
00129     int tagSize  = 4;        // size of DoF tag, i.e., number of fields in the tag
00130     int posScDim = 0;        // position in the tag, counting from 0, of the subcell dim 
00131     int posScOrd = 1;        // position in the tag, counting from 0, of the subcell ordinal
00132     int posDfOrd = 2;        // position in the tag, counting from 0, of DoF ordinal relative to the subcell
00133   
00134     // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration 
00135 
00136     std::vector<int> tags( tagSize * this->getCardinality() );
00137 
00138     // temporarily just put everything on the cell itself
00139     for (int i=0;i<this->getCardinality();i++) {
00140        tags[tagSize*i] = 2;
00141        tags[tagSize*i+1] = 0;
00142        tags[tagSize*i+2] = i;
00143        tags[tagSize*i+3] = this->getCardinality();
00144      }
00145 
00146     Basis<Scalar,ArrayScalar> &xBasis_ = *this->getBases()[0][0];
00147     Basis<Scalar,ArrayScalar> &yBasis_ = *this->getBases()[0][1];
00148 
00149     // now let's try to do it "right"
00150     // let's get the x and y bases and their dof
00151     const std::vector<std::vector<int> >& xdoftags = xBasis_.getAllDofTags();
00152     const std::vector<std::vector<int> >& ydoftags = yBasis_.getAllDofTags();
00153 
00154     std::map<int,std::map<int,int> > total_dof_per_entity;
00155     std::map<int,std::map<int,int> > current_dof_per_entity;
00156 
00157     for (int i=0;i<4;i++) {
00158       total_dof_per_entity[0][i] = 0;
00159       current_dof_per_entity[0][i] = 0;
00160     }
00161     for (int i=0;i<4;i++) {
00162       total_dof_per_entity[1][i] = 0;
00163       current_dof_per_entity[1][i] = 0;
00164     }
00165     total_dof_per_entity[2][0] = 0;
00166     current_dof_per_entity[2][0] = 0;
00167     
00168     // let's tally the total degrees of freedom on each entity
00169     for (int j=0;j<yBasis_.getCardinality();j++) {
00170       const int ydim = ydoftags[j][0];
00171       const int yent = ydoftags[j][1];
00172       for (int i=0;i<xBasis_.getCardinality();i++) {
00173         const int xdim = xdoftags[i][0];
00174         const int xent = xdoftags[i][1];
00175         int dofdim;
00176         int dofent;
00177         ProductTopology::lineProduct2d( xdim , xent , ydim , yent , dofdim , dofent );
00178         total_dof_per_entity[dofdim][dofent] += 1;
00179       }
00180     }
00181 
00182     int tagcur = 0;
00183     for (int j=0;j<yBasis_.getCardinality();j++) {
00184       const int ydim = ydoftags[j][0];
00185       const int yent = ydoftags[j][1];
00186       for (int i=0;i<xBasis_.getCardinality();i++) {
00187         const int xdim = xdoftags[i][0];
00188         const int xent = xdoftags[i][1];
00189         int dofdim;
00190         int dofent;
00191         ProductTopology::lineProduct2d( xdim , xent , ydim , yent , dofdim , dofent );
00192         tags[4*tagcur] = dofdim;
00193         tags[4*tagcur+1] = dofent;
00194         tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent];
00195         current_dof_per_entity[dofdim][dofent]++;
00196         tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent];
00197         tagcur++;
00198       }
00199     }
00200 
00201     setOrdinalTagData(this -> tagToOrdinal_,
00202                       this -> ordinalToTag_,
00203                       &(tags[0]),
00204                       this -> basisCardinality_,
00205                       tagSize,
00206                       posScDim,
00207                       posScOrd,
00208                       posDfOrd);
00209 
00210   }
00211 
00212   template<class Scalar, class ArrayScalar>
00213   void Basis_HGRAD_QUAD_Cn_FEM<Scalar,ArrayScalar>::getValues( ArrayScalar &outputValues ,
00214                                                                const ArrayScalar &inputPoints ,
00215                                                                const EOperator operatorType ) const 
00216   {
00217 #ifdef HAVE_INTREPID_DEBUG
00218     getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
00219                                               inputPoints,
00220                                               operatorType,
00221                                               this -> getBaseCellTopology(),
00222                                               this -> getCardinality() );
00223 #endif
00224 
00225     FieldContainer<Scalar> xInputPoints(inputPoints.dimension(0),1);
00226     FieldContainer<Scalar> yInputPoints(inputPoints.dimension(0),1);
00227 
00228     const Basis<Scalar,ArrayScalar> &xBasis_ = *this->bases_[0][0];
00229     const Basis<Scalar,ArrayScalar> &yBasis_ = *this->bases_[0][1];
00230 
00231     for (int i=0;i<inputPoints.dimension(0);i++) {
00232       xInputPoints(i,0) = inputPoints(i,0);
00233       yInputPoints(i,0) = inputPoints(i,1);
00234     }
00235 
00236     switch (operatorType) {
00237     case OPERATOR_VALUE:
00238       {
00239         FieldContainer<Scalar> xBasisValues(xBasis_.getCardinality(),xInputPoints.dimension(0));
00240         FieldContainer<Scalar> yBasisValues(yBasis_.getCardinality(),yInputPoints.dimension(0));
00241 
00242         xBasis_.getValues(xBasisValues,xInputPoints,OPERATOR_VALUE);
00243         yBasis_.getValues(yBasisValues,yInputPoints,OPERATOR_VALUE);
00244 
00245         int bfcur = 0;
00246         for (int j=0;j<yBasis_.getCardinality();j++) {
00247           for (int i=0;i<xBasis_.getCardinality();i++) {
00248             for (int k=0;k<inputPoints.dimension(0);k++) {
00249               outputValues(bfcur,k) = xBasisValues(i,k) * yBasisValues(j,k);
00250             }
00251             bfcur++;
00252           }
00253         }
00254       }
00255       break;
00256     case OPERATOR_GRAD:
00257     case OPERATOR_D1:
00258       {
00259         FieldContainer<Scalar> xBasisValues(xBasis_.getCardinality(),xInputPoints.dimension(0));
00260         FieldContainer<Scalar> yBasisValues(yBasis_.getCardinality(),yInputPoints.dimension(0));
00261         FieldContainer<Scalar> xBasisDerivs(xBasis_.getCardinality(),xInputPoints.dimension(0),1);
00262         FieldContainer<Scalar> yBasisDerivs(yBasis_.getCardinality(),yInputPoints.dimension(0),1);
00263 
00264         xBasis_.getValues(xBasisValues,xInputPoints,OPERATOR_VALUE);
00265         yBasis_.getValues(yBasisValues,yInputPoints,OPERATOR_VALUE);
00266         xBasis_.getValues(xBasisDerivs,xInputPoints,OPERATOR_D1);
00267         yBasis_.getValues(yBasisDerivs,yInputPoints,OPERATOR_D1);       
00268 
00269         // there are two multiindices: I need the (1,0) and (0,1) derivatives
00270         int bfcur = 0;
00271 
00272         for (int j=0;j<yBasis_.getCardinality();j++) {
00273           for (int i=0;i<xBasis_.getCardinality();i++) {
00274             for (int k=0;k<inputPoints.dimension(0);k++) {
00275               outputValues(bfcur,k,0) = xBasisDerivs(i,k,0) * yBasisValues(j,k);
00276               outputValues(bfcur,k,1) = xBasisValues(i,k) * yBasisDerivs(j,k,0);
00277             }
00278             bfcur++;
00279           }
00280         }
00281       }
00282       break;
00283     case OPERATOR_D2:
00284     case OPERATOR_D3:
00285     case OPERATOR_D4:
00286     case OPERATOR_D5: 
00287     case OPERATOR_D6:
00288     case OPERATOR_D7:
00289     case OPERATOR_D8:
00290     case OPERATOR_D9:
00291     case OPERATOR_D10:
00292       {
00293         FieldContainer<Scalar> xBasisValues(xBasis_.getCardinality(),xInputPoints.dimension(0));
00294         FieldContainer<Scalar> yBasisValues(yBasis_.getCardinality(),yInputPoints.dimension(0));
00295 
00296         Teuchos::Array<int> partialMult;
00297 
00298         for (int d=0;d<getDkCardinality(operatorType,2);d++) {
00299           getDkMultiplicities( partialMult , d , operatorType , 2 );
00300           if (partialMult[0] == 0) {
00301             xBasisValues.resize(xBasis_.getCardinality(),xInputPoints.dimension(0));
00302             xBasis_.getValues( xBasisValues , xInputPoints, OPERATOR_VALUE );
00303           }
00304           else {
00305             xBasisValues.resize(xBasis_.getCardinality(),xInputPoints.dimension(0),1);
00306             EOperator xop = (EOperator) ( (int) OPERATOR_D1 + partialMult[0] - 1 );
00307             xBasis_.getValues( xBasisValues , xInputPoints, xop );
00308             xBasisValues.resize(xBasis_.getCardinality(),xInputPoints.dimension(0));
00309           }
00310           if (partialMult[1] == 0) {
00311             yBasisValues.resize(yBasis_.getCardinality(),yInputPoints.dimension(0));
00312             yBasis_.getValues( yBasisValues , yInputPoints, OPERATOR_VALUE );
00313           }
00314           else {
00315             yBasisValues.resize(yBasis_.getCardinality(),yInputPoints.dimension(0),1);
00316             EOperator yop = (EOperator) ( (int) OPERATOR_D1 + partialMult[1] - 1 );
00317             yBasis_.getValues( yBasisValues , yInputPoints, yop );
00318             yBasisValues.resize(yBasis_.getCardinality(),yInputPoints.dimension(0));
00319           }
00320 
00321 
00322           int bfcur = 0;
00323           for (int j=0;j<yBasis_.getCardinality();j++) {
00324             for (int i=0;i<xBasis_.getCardinality();i++) {
00325               for (int k=0;k<inputPoints.dimension(0);k++) {
00326                 outputValues(bfcur,k,d) = xBasisValues(i,k) * yBasisValues(j,k);
00327               }
00328               bfcur++;
00329             }
00330           }
00331         }
00332       }
00333       break;
00334     case OPERATOR_CURL:
00335       {
00336         FieldContainer<Scalar> xBasisValues(xBasis_.getCardinality(),xInputPoints.dimension(0));
00337         FieldContainer<Scalar> yBasisValues(yBasis_.getCardinality(),yInputPoints.dimension(0));
00338         FieldContainer<Scalar> xBasisDerivs(xBasis_.getCardinality(),xInputPoints.dimension(0),1);
00339         FieldContainer<Scalar> yBasisDerivs(yBasis_.getCardinality(),yInputPoints.dimension(0),1);
00340 
00341         xBasis_.getValues(xBasisValues,xInputPoints,OPERATOR_VALUE);
00342         yBasis_.getValues(yBasisValues,yInputPoints,OPERATOR_VALUE);
00343         xBasis_.getValues(xBasisDerivs,xInputPoints,OPERATOR_D1);
00344         yBasis_.getValues(yBasisDerivs,yInputPoints,OPERATOR_D1);       
00345 
00346         // there are two multiindices: I need the (1,0) and (0,1) derivatives
00347         int bfcur = 0;
00348 
00349         for (int j=0;j<yBasis_.getCardinality();j++) {
00350           for (int i=0;i<xBasis_.getCardinality();i++) {
00351             for (int k=0;k<inputPoints.dimension(0);k++) {
00352               outputValues(bfcur,k,0) = xBasisValues(i,k) * yBasisDerivs(j,k,0);
00353               outputValues(bfcur,k,1) = -xBasisDerivs(i,k,0) * yBasisValues(j,k);
00354             }
00355             bfcur++;
00356           }
00357         }
00358       }
00359       break;      
00360     default:
00361         TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument,
00362                             ">>> ERROR (Basis_HGRAD_QUAD_Cn_FEM): Operator type not implemented");
00363         break;
00364     }
00365   }
00366 
00367   template<class Scalar,class ArrayScalar>
00368   void Basis_HGRAD_QUAD_Cn_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar&           outputValues,
00369                                                                const ArrayScalar &    inputPoints,
00370                                                                const ArrayScalar &    cellVertices,
00371                                                                const EOperator        operatorType) const {
00372     TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
00373                         ">>> ERROR (Basis_HGRAD_QUAD_Cn_FEM): FEM Basis calling an FVD member function");
00374   }
00375 
00376   template<class Scalar,class ArrayScalar>
00377   void Basis_HGRAD_QUAD_Cn_FEM<Scalar, ArrayScalar>::getDofCoords( ArrayScalar & dofCoords ) const
00378   {
00379     int cur = 0;
00380     for (int j=0;j<ptsy_.dimension(0);j++)
00381       {
00382         for (int i=0;i<ptsx_.dimension(0);i++)
00383           {
00384             dofCoords(cur,0) = ptsx_(i);
00385             dofCoords(cur,1) = ptsy_(j);
00386             cur++;
00387           }
00388       }
00389   }
00390 
00391   
00392 }// namespace Intrepid
00393 
00394 #endif