Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/src/Discretization/Basis/Intrepid_HDIV_QUAD_In_FEMDef.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ************************************************************************
00003 //
00004 //                           Intrepid Package
00005 //                 Copyright (2007) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Pavel Bochev  (pbboche@sandia.gov)
00038 //                    Denis Ridzal  (dridzal@sandia.gov), or
00039 //                    Kara Peterson (kjpeter@sandia.gov)
00040 //
00041 // ************************************************************************
00042 // @HEADER
00043 
00049 namespace Intrepid {
00050 
00051   template<class Scalar, class ArrayScalar>
00052   Basis_HDIV_QUAD_In_FEM<Scalar,ArrayScalar>::Basis_HDIV_QUAD_In_FEM( int order ,
00053                                                                       const ArrayScalar & ptsClosed ,
00054                                                                       const ArrayScalar & ptsOpen):
00055     closedBasis_( order , ptsClosed ),
00056     openBasis_( order-1 , ptsOpen ),
00057     closedPts_( ptsClosed ),
00058     openPts_( ptsOpen )
00059   {
00060     this -> basisDegree_       = order;
00061     this -> basisCardinality_  = 2 * closedBasis_.getCardinality() * openBasis_.getCardinality(); 
00062     this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
00063     this -> basisType_         = BASIS_FEM_FIAT;
00064     this -> basisCoordinates_  = COORDINATES_CARTESIAN;
00065     this -> basisTagsAreSet_   = false;
00066   }
00067 
00068   template<class Scalar, class ArrayScalar>
00069   Basis_HDIV_QUAD_In_FEM<Scalar,ArrayScalar>::Basis_HDIV_QUAD_In_FEM( int order , const EPointType &pointType ):
00070     closedBasis_( order , pointType==POINTTYPE_SPECTRAL?POINTTYPE_SPECTRAL:POINTTYPE_EQUISPACED ),
00071     openBasis_( order-1 , pointType==POINTTYPE_SPECTRAL?POINTTYPE_SPECTRAL_OPEN:POINTTYPE_EQUISPACED ),
00072     closedPts_( order+1 , 1 ),
00073     openPts_( order , 1 )
00074   {
00075     this -> basisDegree_       = order;
00076     this -> basisCardinality_  = 2 * closedBasis_.getCardinality() * openBasis_.getCardinality(); 
00077     this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
00078     this -> basisType_         = BASIS_FEM_FIAT;
00079     this -> basisCoordinates_  = COORDINATES_CARTESIAN;
00080     this -> basisTagsAreSet_   = false;
00081 
00082     PointTools::getLattice<Scalar,FieldContainer<Scalar> >( closedPts_ ,
00083                                                             shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >()) ,
00084                                                             order ,
00085                                                             0 ,
00086                                                             pointType==POINTTYPE_SPECTRAL?POINTTYPE_WARPBLEND:POINTTYPE_EQUISPACED );
00087 
00088     if (pointType == POINTTYPE_SPECTRAL)
00089       {
00090         PointTools::getGaussPoints<Scalar,FieldContainer<Scalar> >( openPts_ ,
00091                                                                     order - 1 );
00092       }
00093     else
00094       {
00095         PointTools::getLattice<Scalar,FieldContainer<Scalar> >( openPts_ ,
00096                                                                 shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >()) ,
00097                                                                 order - 1,
00098                                                                 0 ,
00099                                                                 POINTTYPE_EQUISPACED );
00100 
00101       }    
00102 
00103   }
00104   
00105   template<class Scalar, class ArrayScalar>
00106   void Basis_HDIV_QUAD_In_FEM<Scalar, ArrayScalar>::initializeTags() {
00107     
00108     // Basis-dependent intializations
00109     int tagSize  = 4;        // size of DoF tag
00110     int posScDim = 0;        // position in the tag, counting from 0, of the subcell dim 
00111     int posScOrd = 1;        // position in the tag, counting from 0, of the subcell ordinal
00112     int posDfOrd = 2;        // position in the tag, counting from 0, of DoF ordinal relative to the subcell
00113     
00114     std::vector<int> tags( tagSize * this->getCardinality() );
00115     
00116     const std::vector<std::vector<int> >& closedDofTags = closedBasis_.getAllDofTags();
00117     const std::vector<std::vector<int> >& openDofTags = openBasis_.getAllDofTags();
00118 
00119     std::map<int,std::map<int,int> > total_dof_per_entity;
00120     std::map<int,std::map<int,int> > current_dof_per_entity;
00121 
00122     for (int i=0;i<4;i++) {
00123       total_dof_per_entity[0][i] = 0;
00124       current_dof_per_entity[0][i] = 0;
00125     }
00126     for (int i=0;i<4;i++) {
00127       total_dof_per_entity[1][i] = 0;
00128       current_dof_per_entity[1][i] = 0;
00129     }
00130     total_dof_per_entity[2][0] = 0;
00131     current_dof_per_entity[2][0] = 0;
00132 
00133     // tally dof on each facet.  none on vertex
00134     for (int i=0;i<4;i++) {
00135       total_dof_per_entity[1][i] = openBasis_.getCardinality();
00136     }
00137 
00138     total_dof_per_entity[2][0] = this->getCardinality() - 4 * openBasis_.getCardinality();
00139 
00140     int tagcur = 0;
00141     // loop over the x-component basis functions, which are (psi(x)phi(y),0)
00142     // for psi in the closed basis and phi in the open
00143     for (int j=0;j<openBasis_.getCardinality();j++) {
00144       const int odim = openDofTags[j][0];
00145       const int oent = openDofTags[j][1];
00146       for (int i=0;i<closedBasis_.getCardinality();i++) {
00147         const int cdim = closedDofTags[i][0];
00148         const int cent = closedDofTags[i][1];
00149         int dofdim;
00150         int dofent;
00151         ProductTopology::lineProduct2d(cdim,cent,odim,oent,dofdim,dofent);
00152         tags[4*tagcur] = dofdim;
00153         tags[4*tagcur+1] = dofent;
00154         tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent];
00155         current_dof_per_entity[dofdim][dofent]++;
00156         tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent];
00157         tagcur++;
00158       }
00159     }
00160     // now we have to do it for the y-component basis functions, which are
00161     // (0,phi(x)psi(y)) for psi in the closed basis and phi in the open
00162     for (int j=0;j<closedBasis_.getCardinality();j++) {
00163       const int cdim = closedDofTags[j][0];
00164       const int cent = closedDofTags[j][1];
00165       for (int i=0;i<openBasis_.getCardinality();i++) {
00166         const int odim = openDofTags[i][0];
00167         const int oent = openDofTags[i][1];
00168         int dofdim;
00169         int dofent;
00170         ProductTopology::lineProduct2d(odim,oent,cdim,cent,dofdim,dofent);
00171         tags[4*tagcur] = dofdim;
00172         tags[4*tagcur+1] = dofent;
00173         tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent];
00174         current_dof_per_entity[dofdim][dofent]++;
00175         tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent];
00176         tagcur++;
00177       }
00178     }
00179 
00180 //     for (int i=0;i<this->getCardinality();i++) {
00181 //       for (int j=0;j<4;j++) {
00182 //      std::cout << tags[4*i+j] << " ";
00183 //       }
00184 //       std::cout << std::endl;
00185 //     }
00186   
00187     // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
00188     Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
00189                                 this -> ordinalToTag_,
00190                                 &(tags[0]),
00191                                 this -> basisCardinality_,
00192                                 tagSize,
00193                                 posScDim,
00194                                 posScOrd,
00195                                 posDfOrd);
00196   }
00197 
00198 
00199   template<class Scalar, class ArrayScalar>
00200   void Basis_HDIV_QUAD_In_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar &        outputValues,
00201                                                               const ArrayScalar &  inputPoints,
00202                                                               const EOperator      operatorType) const {
00203     
00204     // Verify arguments
00205 #ifdef HAVE_INTREPID_DEBUG
00206     Intrepid::getValues_HDIV_Args<Scalar, ArrayScalar>(outputValues,
00207                                                        inputPoints,
00208                                                        operatorType,
00209                                                        this -> getBaseCellTopology(),
00210                                                        this -> getCardinality() );
00211 #endif
00212     
00213     // Number of evaluation points = dim 0 of inputPoints
00214     int dim0 = inputPoints.dimension(0);
00215     
00216     // separate out points
00217     FieldContainer<Scalar> xPoints(dim0,1);
00218     FieldContainer<Scalar> yPoints(dim0,1);
00219     
00220     for (int i=0;i<dim0;i++) {
00221       xPoints(i,0) = inputPoints(i,0);
00222       yPoints(i,0) = inputPoints(i,1);
00223     }
00224     
00225     switch (operatorType) {
00226     case OPERATOR_VALUE:
00227       {
00228         FieldContainer<Scalar> closedBasisValsXPts( closedBasis_.getCardinality() , dim0 );
00229         FieldContainer<Scalar> closedBasisValsYPts( closedBasis_.getCardinality() , dim0 );
00230         FieldContainer<Scalar> openBasisValsXPts( openBasis_.getCardinality() , dim0 );
00231         FieldContainer<Scalar> openBasisValsYPts( openBasis_.getCardinality() , dim0 );
00232         
00233         closedBasis_.getValues( closedBasisValsXPts , xPoints , OPERATOR_VALUE );
00234         closedBasis_.getValues( closedBasisValsYPts , yPoints , OPERATOR_VALUE );
00235         openBasis_.getValues( openBasisValsXPts , xPoints , OPERATOR_VALUE );
00236         openBasis_.getValues( openBasisValsYPts , yPoints , OPERATOR_VALUE );
00237         
00238         int bfcur = 0;
00239         // x component bfs are (closed(x) open(y),0)
00240         for (int j=0;j<openBasis_.getCardinality();j++) {
00241           for (int i=0;i<closedBasis_.getCardinality();i++) {
00242             for (int l=0;l<dim0;l++) {
00243               outputValues(bfcur,l,0) = closedBasisValsXPts(i,l) * openBasisValsYPts(j,l);
00244               outputValues(bfcur,l,1) = 0.0;
00245             }
00246             bfcur++;
00247           }
00248         }
00249         
00250         // y component bfs are (0,open(x) closed(y))
00251         for (int j=0;j<closedBasis_.getCardinality();j++) {
00252           for (int i=0;i<openBasis_.getCardinality();i++) {
00253             for (int l=0;l<dim0;l++) {
00254               outputValues(bfcur,l,0) = 0.0;
00255               outputValues(bfcur,l,1) = openBasisValsXPts(i,l) * closedBasisValsYPts(j,l);
00256             }
00257             bfcur++;
00258           }
00259         }
00260       }
00261       break;
00262     case OPERATOR_DIV:
00263       {
00264         FieldContainer<Scalar> closedBasisDerivsXPts( closedBasis_.getCardinality() , dim0 , 1 );
00265         FieldContainer<Scalar> closedBasisDerivsYPts( closedBasis_.getCardinality() , dim0 , 1 );
00266         FieldContainer<Scalar> openBasisValsXPts( openBasis_.getCardinality() , dim0 );
00267         FieldContainer<Scalar> openBasisValsYPts( openBasis_.getCardinality() , dim0 );
00268         
00269         closedBasis_.getValues( closedBasisDerivsXPts , xPoints , OPERATOR_D1 );
00270         closedBasis_.getValues( closedBasisDerivsYPts , yPoints , OPERATOR_D1 );
00271         openBasis_.getValues( openBasisValsXPts , xPoints , OPERATOR_VALUE );
00272         openBasis_.getValues( openBasisValsYPts , yPoints , OPERATOR_VALUE );
00273         
00274         int bfcur = 0;
00275         
00276         // x component basis functions first
00277         for (int j=0;j<openBasis_.getCardinality();j++) {
00278           for (int i=0;i<closedBasis_.getCardinality();i++) {
00279             for (int l=0;l<dim0;l++) {
00280               outputValues(bfcur,l) = closedBasisDerivsXPts(i,l,0) * openBasisValsYPts(j,l);
00281             }
00282             bfcur++;
00283           }
00284         }
00285         
00286         // now y component basis functions
00287         for (int j=0;j<closedBasis_.getCardinality();j++) {
00288           for (int i=0;i<openBasis_.getCardinality();i++) {
00289             for (int l=0;l<dim0;l++) {
00290               outputValues(bfcur,l) = openBasisValsXPts(i,l) * closedBasisDerivsYPts(j,l,0);
00291             }
00292             bfcur++;
00293           }
00294         }
00295       }
00296       break;
00297     case OPERATOR_CURL:
00298       TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_CURL), std::invalid_argument,
00299                           ">>> ERROR (Basis_HDIV_QUAD_In_FEM): CURL is invalid operator for HDIV Basis Functions");
00300       break;
00301       
00302     case OPERATOR_GRAD:
00303       TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_GRAD), std::invalid_argument,
00304                           ">>> ERROR (Basis_HDIV_QUAD_In_FEM): GRAD is invalid operator for HDIV Basis Functions");
00305       break;
00306       
00307     case OPERATOR_D1:
00308     case OPERATOR_D2:
00309     case OPERATOR_D3:
00310     case OPERATOR_D4:
00311     case OPERATOR_D5:
00312     case OPERATOR_D6:
00313     case OPERATOR_D7:
00314     case OPERATOR_D8:
00315     case OPERATOR_D9:
00316     case OPERATOR_D10:
00317       TEUCHOS_TEST_FOR_EXCEPTION( ( (operatorType == OPERATOR_D1)    ||
00318                             (operatorType == OPERATOR_D2)    ||
00319                             (operatorType == OPERATOR_D3)    ||
00320                             (operatorType == OPERATOR_D4)    ||
00321                             (operatorType == OPERATOR_D5)    ||
00322                             (operatorType == OPERATOR_D6)    ||
00323                             (operatorType == OPERATOR_D7)    ||
00324                             (operatorType == OPERATOR_D8)    ||
00325                             (operatorType == OPERATOR_D9)    ||
00326                             (operatorType == OPERATOR_D10) ),
00327                           std::invalid_argument,
00328                           ">>> ERROR (Basis_HDIV_QUAD_In_FEM): Invalid operator type");
00329       break;
00330       
00331     default:
00332       TEUCHOS_TEST_FOR_EXCEPTION( ( (operatorType != OPERATOR_VALUE) &&
00333                             (operatorType != OPERATOR_GRAD)  &&
00334                             (operatorType != OPERATOR_CURL)  &&
00335                             (operatorType != OPERATOR_DIV)   &&
00336                             (operatorType != OPERATOR_D1)    &&
00337                             (operatorType != OPERATOR_D2)    &&
00338                             (operatorType != OPERATOR_D3)    &&
00339                             (operatorType != OPERATOR_D4)    &&
00340                             (operatorType != OPERATOR_D5)    &&
00341                             (operatorType != OPERATOR_D6)    &&
00342                             (operatorType != OPERATOR_D7)    &&
00343                             (operatorType != OPERATOR_D8)    &&
00344                             (operatorType != OPERATOR_D9)    &&
00345                             (operatorType != OPERATOR_D10) ),
00346                           std::invalid_argument,
00347                           ">>> ERROR (Basis_HDIV_QUAD_In_FEM): Invalid operator type");
00348     }
00349   }
00350   
00351   
00352   
00353   template<class Scalar, class ArrayScalar>
00354   void Basis_HDIV_QUAD_In_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar&           outputValues,
00355                                                               const ArrayScalar &    inputPoints,
00356                                                               const ArrayScalar &    cellVertices,
00357                                                               const EOperator        operatorType) const {
00358     TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
00359                         ">>> ERROR (Basis_HDIV_QUAD_In_FEM): FEM Basis calling an FVD member function");
00360   }
00361   
00362   template<class Scalar, class ArrayScalar>
00363   void Basis_HDIV_QUAD_In_FEM<Scalar,ArrayScalar>::getDofCoords(ArrayScalar & DofCoords) const
00364   {
00365     // x-component basis functions
00366     int cur = 0;
00367 
00368     for (int j=0;j<openPts_.dimension(0);j++)
00369       {
00370         for (int i=0;i<closedPts_.dimension(0);i++)
00371           {
00372             DofCoords(cur,0) = closedPts_(i,0);
00373             DofCoords(cur,1) = openPts_(j,0);
00374             cur++;
00375           }
00376       }
00377 
00378     // y-component basis functions
00379     for (int j=0;j<closedPts_.dimension(0);j++)
00380       {
00381         for (int i=0;i<openPts_.dimension(0);i++)
00382           {
00383             DofCoords(cur,0) = openPts_(i,0);
00384             DofCoords(cur,1) = closedPts_(j,0);
00385             cur++;
00386           }
00387       }
00388 
00389     return;
00390   }
00391 
00392   
00393 }// namespace Intrepid