Intrepid
http://trilinos.sandia.gov/packages/docs/r11.2/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     Array<Array<RCP<Basis<Scalar,ArrayScalar > > > > bases(2);
00068     bases[0].resize(2); bases[1].resize(2);
00069     bases[0][0] = rcp( &closedBasis_ , false );
00070     bases[0][1] = rcp( &openBasis_ , false );
00071     bases[1][0] = rcp( &openBasis_ , false );
00072     bases[1][1] = rcp( &closedBasis_ , false );
00073     this->setBases( bases );
00074 
00075   }
00076 
00077   template<class Scalar, class ArrayScalar>
00078   Basis_HDIV_QUAD_In_FEM<Scalar,ArrayScalar>::Basis_HDIV_QUAD_In_FEM( int order , const EPointType &pointType ):
00079     closedBasis_( order , pointType==POINTTYPE_SPECTRAL?POINTTYPE_SPECTRAL:POINTTYPE_EQUISPACED ),
00080     openBasis_( order-1 , pointType==POINTTYPE_SPECTRAL?POINTTYPE_SPECTRAL_OPEN:POINTTYPE_EQUISPACED ),
00081     closedPts_( order+1 , 1 ),
00082     openPts_( order , 1 )
00083   {
00084     this -> basisDegree_       = order;
00085     this -> basisCardinality_  = 2 * closedBasis_.getCardinality() * openBasis_.getCardinality(); 
00086     this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
00087     this -> basisType_         = BASIS_FEM_FIAT;
00088     this -> basisCoordinates_  = COORDINATES_CARTESIAN;
00089     this -> basisTagsAreSet_   = false;
00090 
00091     PointTools::getLattice<Scalar,FieldContainer<Scalar> >( closedPts_ ,
00092                                                             shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >()) ,
00093                                                             order ,
00094                                                             0 ,
00095                                                             pointType==POINTTYPE_SPECTRAL?POINTTYPE_WARPBLEND:POINTTYPE_EQUISPACED );
00096 
00097     if (pointType == POINTTYPE_SPECTRAL)
00098       {
00099         PointTools::getGaussPoints<Scalar,FieldContainer<Scalar> >( openPts_ ,
00100                                                                     order - 1 );
00101       }
00102     else
00103       {
00104         PointTools::getLattice<Scalar,FieldContainer<Scalar> >( openPts_ ,
00105                                                                 shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >()) ,
00106                                                                 order - 1,
00107                                                                 0 ,
00108                                                                 POINTTYPE_EQUISPACED );
00109 
00110       }
00111 
00112     Array<Array<RCP<Basis<Scalar,ArrayScalar > > > > bases(2);
00113     bases[0].resize(2); bases[1].resize(2);
00114     bases[0][0] = rcp( &closedBasis_ , false );
00115     bases[0][1] = rcp( &openBasis_ , false );
00116     bases[1][0] = rcp( &openBasis_ , false );
00117     bases[1][1] = rcp( &closedBasis_ , false );
00118     this->setBases( bases );
00119 
00120   }
00121   
00122   template<class Scalar, class ArrayScalar>
00123   void Basis_HDIV_QUAD_In_FEM<Scalar, ArrayScalar>::initializeTags() {
00124     
00125     // Basis-dependent intializations
00126     int tagSize  = 4;        // size of DoF tag
00127     int posScDim = 0;        // position in the tag, counting from 0, of the subcell dim 
00128     int posScOrd = 1;        // position in the tag, counting from 0, of the subcell ordinal
00129     int posDfOrd = 2;        // position in the tag, counting from 0, of DoF ordinal relative to the subcell
00130     
00131     std::vector<int> tags( tagSize * this->getCardinality() );
00132     
00133     const std::vector<std::vector<int> >& closedDofTags = closedBasis_.getAllDofTags();
00134     const std::vector<std::vector<int> >& openDofTags = openBasis_.getAllDofTags();
00135 
00136     std::map<int,std::map<int,int> > total_dof_per_entity;
00137     std::map<int,std::map<int,int> > current_dof_per_entity;
00138 
00139     for (int i=0;i<4;i++) {
00140       total_dof_per_entity[0][i] = 0;
00141       current_dof_per_entity[0][i] = 0;
00142     }
00143     for (int i=0;i<4;i++) {
00144       total_dof_per_entity[1][i] = 0;
00145       current_dof_per_entity[1][i] = 0;
00146     }
00147     total_dof_per_entity[2][0] = 0;
00148     current_dof_per_entity[2][0] = 0;
00149 
00150     // tally dof on each facet.  none on vertex
00151     for (int i=0;i<4;i++) {
00152       total_dof_per_entity[1][i] = openBasis_.getCardinality();
00153     }
00154 
00155     total_dof_per_entity[2][0] = this->getCardinality() - 4 * openBasis_.getCardinality();
00156 
00157     int tagcur = 0;
00158     // loop over the x-component basis functions, which are (psi(x)phi(y),0)
00159     // for psi in the closed basis and phi in the open
00160     for (int j=0;j<openBasis_.getCardinality();j++) {
00161       const int odim = openDofTags[j][0];
00162       const int oent = openDofTags[j][1];
00163       for (int i=0;i<closedBasis_.getCardinality();i++) {
00164         const int cdim = closedDofTags[i][0];
00165         const int cent = closedDofTags[i][1];
00166         int dofdim;
00167         int dofent;
00168         ProductTopology::lineProduct2d(cdim,cent,odim,oent,dofdim,dofent);
00169         tags[4*tagcur] = dofdim;
00170         tags[4*tagcur+1] = dofent;
00171         tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent];
00172         current_dof_per_entity[dofdim][dofent]++;
00173         tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent];
00174         tagcur++;
00175       }
00176     }
00177     // now we have to do it for the y-component basis functions, which are
00178     // (0,phi(x)psi(y)) for psi in the closed basis and phi in the open
00179     for (int j=0;j<closedBasis_.getCardinality();j++) {
00180       const int cdim = closedDofTags[j][0];
00181       const int cent = closedDofTags[j][1];
00182       for (int i=0;i<openBasis_.getCardinality();i++) {
00183         const int odim = openDofTags[i][0];
00184         const int oent = openDofTags[i][1];
00185         int dofdim;
00186         int dofent;
00187         ProductTopology::lineProduct2d(odim,oent,cdim,cent,dofdim,dofent);
00188         tags[4*tagcur] = dofdim;
00189         tags[4*tagcur+1] = dofent;
00190         tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent];
00191         current_dof_per_entity[dofdim][dofent]++;
00192         tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent];
00193         tagcur++;
00194       }
00195     }
00196 
00197 //     for (int i=0;i<this->getCardinality();i++) {
00198 //       for (int j=0;j<4;j++) {
00199 //      std::cout << tags[4*i+j] << " ";
00200 //       }
00201 //       std::cout << std::endl;
00202 //     }
00203   
00204     // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
00205     Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
00206                                 this -> ordinalToTag_,
00207                                 &(tags[0]),
00208                                 this -> basisCardinality_,
00209                                 tagSize,
00210                                 posScDim,
00211                                 posScOrd,
00212                                 posDfOrd);
00213   }
00214 
00215 
00216   template<class Scalar, class ArrayScalar>
00217   void Basis_HDIV_QUAD_In_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar &        outputValues,
00218                                                               const ArrayScalar &  inputPoints,
00219                                                               const EOperator      operatorType) const {
00220     
00221     // Verify arguments
00222 #ifdef HAVE_INTREPID_DEBUG
00223     Intrepid::getValues_HDIV_Args<Scalar, ArrayScalar>(outputValues,
00224                                                        inputPoints,
00225                                                        operatorType,
00226                                                        this -> getBaseCellTopology(),
00227                                                        this -> getCardinality() );
00228 #endif
00229     
00230     // Number of evaluation points = dim 0 of inputPoints
00231     int dim0 = inputPoints.dimension(0);
00232     
00233     // separate out points
00234     FieldContainer<Scalar> xPoints(dim0,1);
00235     FieldContainer<Scalar> yPoints(dim0,1);
00236     
00237     for (int i=0;i<dim0;i++) {
00238       xPoints(i,0) = inputPoints(i,0);
00239       yPoints(i,0) = inputPoints(i,1);
00240     }
00241     
00242     switch (operatorType) {
00243     case OPERATOR_VALUE:
00244       {
00245         FieldContainer<Scalar> closedBasisValsXPts( closedBasis_.getCardinality() , dim0 );
00246         FieldContainer<Scalar> closedBasisValsYPts( closedBasis_.getCardinality() , dim0 );
00247         FieldContainer<Scalar> openBasisValsXPts( openBasis_.getCardinality() , dim0 );
00248         FieldContainer<Scalar> openBasisValsYPts( openBasis_.getCardinality() , dim0 );
00249         
00250         closedBasis_.getValues( closedBasisValsXPts , xPoints , OPERATOR_VALUE );
00251         closedBasis_.getValues( closedBasisValsYPts , yPoints , OPERATOR_VALUE );
00252         openBasis_.getValues( openBasisValsXPts , xPoints , OPERATOR_VALUE );
00253         openBasis_.getValues( openBasisValsYPts , yPoints , OPERATOR_VALUE );
00254         
00255         int bfcur = 0;
00256         // x component bfs are (closed(x) open(y),0)
00257         for (int j=0;j<openBasis_.getCardinality();j++) {
00258           for (int i=0;i<closedBasis_.getCardinality();i++) {
00259             for (int l=0;l<dim0;l++) {
00260               outputValues(bfcur,l,0) = closedBasisValsXPts(i,l) * openBasisValsYPts(j,l);
00261               outputValues(bfcur,l,1) = 0.0;
00262             }
00263             bfcur++;
00264           }
00265         }
00266         
00267         // y component bfs are (0,open(x) closed(y))
00268         for (int j=0;j<closedBasis_.getCardinality();j++) {
00269           for (int i=0;i<openBasis_.getCardinality();i++) {
00270             for (int l=0;l<dim0;l++) {
00271               outputValues(bfcur,l,0) = 0.0;
00272               outputValues(bfcur,l,1) = openBasisValsXPts(i,l) * closedBasisValsYPts(j,l);
00273             }
00274             bfcur++;
00275           }
00276         }
00277       }
00278       break;
00279     case OPERATOR_DIV:
00280       {
00281         FieldContainer<Scalar> closedBasisDerivsXPts( closedBasis_.getCardinality() , dim0 , 1 );
00282         FieldContainer<Scalar> closedBasisDerivsYPts( closedBasis_.getCardinality() , dim0 , 1 );
00283         FieldContainer<Scalar> openBasisValsXPts( openBasis_.getCardinality() , dim0 );
00284         FieldContainer<Scalar> openBasisValsYPts( openBasis_.getCardinality() , dim0 );
00285         
00286         closedBasis_.getValues( closedBasisDerivsXPts , xPoints , OPERATOR_D1 );
00287         closedBasis_.getValues( closedBasisDerivsYPts , yPoints , OPERATOR_D1 );
00288         openBasis_.getValues( openBasisValsXPts , xPoints , OPERATOR_VALUE );
00289         openBasis_.getValues( openBasisValsYPts , yPoints , OPERATOR_VALUE );
00290         
00291         int bfcur = 0;
00292         
00293         // x component basis functions first
00294         for (int j=0;j<openBasis_.getCardinality();j++) {
00295           for (int i=0;i<closedBasis_.getCardinality();i++) {
00296             for (int l=0;l<dim0;l++) {
00297               outputValues(bfcur,l) = closedBasisDerivsXPts(i,l,0) * openBasisValsYPts(j,l);
00298             }
00299             bfcur++;
00300           }
00301         }
00302         
00303         // now y component basis functions
00304         for (int j=0;j<closedBasis_.getCardinality();j++) {
00305           for (int i=0;i<openBasis_.getCardinality();i++) {
00306             for (int l=0;l<dim0;l++) {
00307               outputValues(bfcur,l) = openBasisValsXPts(i,l) * closedBasisDerivsYPts(j,l,0);
00308             }
00309             bfcur++;
00310           }
00311         }
00312       }
00313       break;
00314     case OPERATOR_CURL:
00315       TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_CURL), std::invalid_argument,
00316                           ">>> ERROR (Basis_HDIV_QUAD_In_FEM): CURL is invalid operator for HDIV Basis Functions");
00317       break;
00318       
00319     case OPERATOR_GRAD:
00320       TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_GRAD), std::invalid_argument,
00321                           ">>> ERROR (Basis_HDIV_QUAD_In_FEM): GRAD is invalid operator for HDIV Basis Functions");
00322       break;
00323       
00324     case OPERATOR_D1:
00325     case OPERATOR_D2:
00326     case OPERATOR_D3:
00327     case OPERATOR_D4:
00328     case OPERATOR_D5:
00329     case OPERATOR_D6:
00330     case OPERATOR_D7:
00331     case OPERATOR_D8:
00332     case OPERATOR_D9:
00333     case OPERATOR_D10:
00334       TEUCHOS_TEST_FOR_EXCEPTION( ( (operatorType == OPERATOR_D1)    ||
00335                             (operatorType == OPERATOR_D2)    ||
00336                             (operatorType == OPERATOR_D3)    ||
00337                             (operatorType == OPERATOR_D4)    ||
00338                             (operatorType == OPERATOR_D5)    ||
00339                             (operatorType == OPERATOR_D6)    ||
00340                             (operatorType == OPERATOR_D7)    ||
00341                             (operatorType == OPERATOR_D8)    ||
00342                             (operatorType == OPERATOR_D9)    ||
00343                             (operatorType == OPERATOR_D10) ),
00344                           std::invalid_argument,
00345                           ">>> ERROR (Basis_HDIV_QUAD_In_FEM): Invalid operator type");
00346       break;
00347       
00348     default:
00349       TEUCHOS_TEST_FOR_EXCEPTION( ( (operatorType != OPERATOR_VALUE) &&
00350                             (operatorType != OPERATOR_GRAD)  &&
00351                             (operatorType != OPERATOR_CURL)  &&
00352                             (operatorType != OPERATOR_DIV)   &&
00353                             (operatorType != OPERATOR_D1)    &&
00354                             (operatorType != OPERATOR_D2)    &&
00355                             (operatorType != OPERATOR_D3)    &&
00356                             (operatorType != OPERATOR_D4)    &&
00357                             (operatorType != OPERATOR_D5)    &&
00358                             (operatorType != OPERATOR_D6)    &&
00359                             (operatorType != OPERATOR_D7)    &&
00360                             (operatorType != OPERATOR_D8)    &&
00361                             (operatorType != OPERATOR_D9)    &&
00362                             (operatorType != OPERATOR_D10) ),
00363                           std::invalid_argument,
00364                           ">>> ERROR (Basis_HDIV_QUAD_In_FEM): Invalid operator type");
00365     }
00366   }
00367   
00368   
00369   
00370   template<class Scalar, class ArrayScalar>
00371   void Basis_HDIV_QUAD_In_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar&           outputValues,
00372                                                               const ArrayScalar &    inputPoints,
00373                                                               const ArrayScalar &    cellVertices,
00374                                                               const EOperator        operatorType) const {
00375     TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
00376                         ">>> ERROR (Basis_HDIV_QUAD_In_FEM): FEM Basis calling an FVD member function");
00377   }
00378   
00379   template<class Scalar, class ArrayScalar>
00380   void Basis_HDIV_QUAD_In_FEM<Scalar,ArrayScalar>::getDofCoords(ArrayScalar & DofCoords) const
00381   {
00382     // x-component basis functions
00383     int cur = 0;
00384 
00385     for (int j=0;j<openPts_.dimension(0);j++)
00386       {
00387         for (int i=0;i<closedPts_.dimension(0);i++)
00388           {
00389             DofCoords(cur,0) = closedPts_(i,0);
00390             DofCoords(cur,1) = openPts_(j,0);
00391             cur++;
00392           }
00393       }
00394 
00395     // y-component basis functions
00396     for (int j=0;j<closedPts_.dimension(0);j++)
00397       {
00398         for (int i=0;i<openPts_.dimension(0);i++)
00399           {
00400             DofCoords(cur,0) = openPts_(i,0);
00401             DofCoords(cur,1) = closedPts_(j,0);
00402             cur++;
00403           }
00404       }
00405 
00406     return;
00407   }
00408 
00409   
00410 }// namespace Intrepid