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