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