Intrepid
http://trilinos.sandia.gov/packages/docs/r11.2/packages/intrepid/src/Discretization/Basis/Intrepid_HGRAD_LINE_Cn_FEMDef.hpp
00001 #ifndef INTREPID_HGRAD_LINE_CN_FEMDEF_HPP
00002 #define INTREPID_HGRAD_LINE_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 namespace Intrepid {
00052 
00053   template<class Scalar, class ArrayScalar>
00054   Basis_HGRAD_LINE_Cn_FEM<Scalar,ArrayScalar>::Basis_HGRAD_LINE_Cn_FEM( const int n ,
00055                                                                         const ArrayScalar &pts ):
00056     latticePts_( n+1 , 1 ),
00057     Phis_( n ),
00058     V_(n+1,n+1),
00059     Vinv_(n+1,n+1)
00060   {
00061     const int N = n+1;
00062     this -> basisCardinality_  = N;
00063     this -> basisDegree_       = n;
00064     this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >() );
00065     this -> basisType_         = BASIS_FEM_FIAT;
00066     this -> basisCoordinates_  = COORDINATES_CARTESIAN;
00067     this -> basisTagsAreSet_   = false;
00068 
00069 
00070     // check validity of points
00071     for (int i=0;i<n;i++) {
00072       TEUCHOS_TEST_FOR_EXCEPTION( pts(i,0) >= pts(i+1,0) ,
00073                           std::runtime_error ,
00074                           "Intrepid::Basis_HGRAD_LINE_Cn_FEM Illegal points given to constructor" );
00075     }
00076 
00077     // copy points int latticePts, correcting endpoints if needed
00078     if (std::abs(pts(0,0)+1.0) < INTREPID_TOL) {
00079       latticePts_(0,0) = -1.0;
00080     }
00081     else {
00082       latticePts_(0,0) = pts(0,0);
00083     }
00084     for (int i=1;i<n;i++) {
00085       latticePts_(i,0) = pts(i,0);
00086     }
00087     if (std::abs(pts(n,0)-1.0) < INTREPID_TOL) {
00088       latticePts_(n,0) = 1.0;
00089     }
00090     else {
00091       latticePts_(n,0) = pts(n,0);
00092     }
00093     
00094     // form Vandermonde matrix.  Actually, this is the transpose of the VDM,
00095     // so we transpose on copy below.
00096   
00097     Phis_.getValues( V_ , latticePts_ , OPERATOR_VALUE );
00098 
00099     // now I need to copy V into a Teuchos array to do the inversion
00100     Teuchos::SerialDenseMatrix<int,Scalar> Vsdm(N,N);
00101     for (int i=0;i<N;i++) {
00102       for (int j=0;j<N;j++) {
00103         Vsdm(i,j) = V_(i,j);
00104       }
00105     }
00106 
00107     // invert the matrix
00108     Teuchos::SerialDenseSolver<int,Scalar> solver;
00109     solver.setMatrix( rcp( &Vsdm , false ) );
00110     solver.invert( );
00111 
00112     // now I need to copy the inverse into Vinv
00113     for (int i=0;i<N;i++) {
00114       for (int j=0;j<N;j++) {
00115         Vinv_(i,j) = Vsdm(j,i);
00116       }
00117     }
00118 
00119   }  
00120 
00121   template<class Scalar, class ArrayScalar>
00122   Basis_HGRAD_LINE_Cn_FEM<Scalar,ArrayScalar>::Basis_HGRAD_LINE_Cn_FEM( const int n ,
00123                                                                         const EPointType &pointType ):
00124     latticePts_( n+1 , 1 ),
00125     Phis_( n ),
00126     V_(n+1,n+1),
00127     Vinv_(n+1,n+1)
00128   {
00129     const int N = n+1;
00130     this -> basisCardinality_  = N;
00131     this -> basisDegree_       = n;
00132     this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >() );
00133     this -> basisType_         = BASIS_FEM_FIAT;
00134     this -> basisCoordinates_  = COORDINATES_CARTESIAN;
00135     this -> basisTagsAreSet_   = false;
00136 
00137     switch(pointType) {
00138     case POINTTYPE_EQUISPACED:
00139       PointTools::getLattice<Scalar,FieldContainer<Scalar> >( latticePts_ ,  this->basisCellTopology_ , n , 0 , POINTTYPE_EQUISPACED );
00140       break;
00141     case POINTTYPE_SPECTRAL: 
00142       PointTools::getLattice<Scalar,FieldContainer<Scalar> >( latticePts_ ,  this->basisCellTopology_ , n , 0 , POINTTYPE_WARPBLEND );
00143       break;
00144     case POINTTYPE_SPECTRAL_OPEN: 
00145       PointTools::getGaussPoints<Scalar,FieldContainer<Scalar> >( latticePts_ , n );
00146       break;
00147     default:
00148       TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument , "Basis_HGRAD_LINE_Cn_FEM:: invalid point type" );
00149       break;
00150     }
00151 
00152     // form Vandermonde matrix.  Actually, this is the transpose of the VDM,
00153     // so we transpose on copy below.
00154   
00155     Phis_.getValues( V_ , latticePts_ , OPERATOR_VALUE );
00156 
00157     // now I need to copy V into a Teuchos array to do the inversion
00158     Teuchos::SerialDenseMatrix<int,Scalar> Vsdm(N,N);
00159     for (int i=0;i<N;i++) {
00160       for (int j=0;j<N;j++) {
00161         Vsdm(i,j) = V_(i,j);
00162       }
00163     }
00164 
00165     // invert the matrix
00166     Teuchos::SerialDenseSolver<int,Scalar> solver;
00167     solver.setMatrix( rcp( &Vsdm , false ) );
00168     solver.invert( );
00169 
00170     // now I need to copy the inverse into Vinv
00171     for (int i=0;i<N;i++) {
00172       for (int j=0;j<N;j++) {
00173         Vinv_(i,j) = Vsdm(j,i);
00174       }
00175     }
00176   }  
00177   
00178   
00179   template<class Scalar, class ArrayScalar>
00180   void Basis_HGRAD_LINE_Cn_FEM<Scalar, ArrayScalar>::initializeTags() {
00181   
00182     // Basis-dependent initializations
00183     int tagSize  = 4;        // size of DoF tag, i.e., number of fields in the tag
00184     int posScDim = 0;        // position in the tag, counting from 0, of the subcell dim 
00185     int posScOrd = 1;        // position in the tag, counting from 0, of the subcell ordinal
00186     int posDfOrd = 2;        // position in the tag, counting from 0, of DoF ordinal relative to the subcell
00187   
00188     // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration 
00189 
00190     int *tags = new int[ tagSize * this->getCardinality() ];
00191 
00192     int internal_dof;
00193     int edge_dof;
00194 
00195     const int n = this->getDegree();
00196 
00197     // now we check the points for association 
00198     if (latticePts_(0,0) == -1.0) {
00199       tags[0] = 0;
00200       tags[1] = 0;
00201       tags[2] = 0;
00202       tags[3] = 1;
00203       edge_dof = 1;
00204       internal_dof = n-1;
00205     }
00206     else {
00207       tags[0] = 1;
00208       tags[1] = 0;
00209       tags[2] = 0;
00210       tags[3] = n+1;
00211       edge_dof = 0;
00212       internal_dof = n+1;
00213     }
00214     for (int i=1;i<n;i++) {
00215       tags[4*i] = 1;
00216       tags[4*i+1] = 0;
00217       tags[4*i+2] = -edge_dof + i;
00218       tags[4*i+3] = internal_dof;
00219     }
00220     if (latticePts_(n,0) == 1.0) {
00221       tags[4*n] = 0;
00222       tags[4*n+1] = 1;
00223       tags[4*n+2] = 0;
00224       tags[4*n+3] = 1;
00225     }
00226     else {
00227       tags[4*n] = 1;
00228       tags[4*n+1] = 0;
00229       tags[4*n+2] = n;
00230       tags[4*n+3] = n;
00231     }    
00232     
00233     Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
00234                                 this -> ordinalToTag_,
00235                                 tags,
00236                                 this -> basisCardinality_,
00237                                 tagSize,
00238                                 posScDim,
00239                                 posScOrd,
00240                                 posDfOrd);
00241 
00242     delete []tags;
00243   
00244   }  
00245 
00246 
00247 
00248   template<class Scalar, class ArrayScalar> 
00249   void Basis_HGRAD_LINE_Cn_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar &        outputValues,
00250                                                               const ArrayScalar &  inputPoints,
00251                                                               const EOperator      operatorType) const {
00252   
00253     // Verify arguments
00254 #ifdef HAVE_INTREPID_DEBUG
00255     Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
00256                                                         inputPoints,
00257                                                         operatorType,
00258                                                         this -> getBaseCellTopology(),
00259                                                         this -> getCardinality() );
00260 #endif
00261     const int numPts = inputPoints.dimension(0);
00262     const int numBf = this->getCardinality();
00263 
00264     try {
00265       switch (operatorType) {
00266       case OPERATOR_VALUE:
00267         {
00268           FieldContainer<Scalar> phisCur( numBf , numPts );
00269           Phis_.getValues( phisCur , inputPoints , operatorType );
00270           for (int i=0;i<outputValues.dimension(0);i++) {
00271             for (int j=0;j<outputValues.dimension(1);j++) {
00272               outputValues(i,j) = 0.0;
00273               for (int k=0;k<this->getCardinality();k++) {
00274                 outputValues(i,j) += this->Vinv_(k,i) * phisCur(k,j);
00275               }
00276             }
00277           }
00278         }
00279         break;
00280       case OPERATOR_GRAD:
00281       case OPERATOR_D1:
00282       case OPERATOR_D2:
00283       case OPERATOR_D3:
00284       case OPERATOR_D4:
00285       case OPERATOR_D5:
00286       case OPERATOR_D6:
00287       case OPERATOR_D7:
00288       case OPERATOR_D8:
00289       case OPERATOR_D9:
00290       case OPERATOR_D10:
00291         {
00292           const int dkcard = 
00293             (operatorType == OPERATOR_GRAD)? getDkCardinality(OPERATOR_D1,1): getDkCardinality(operatorType,1);
00294           
00295           FieldContainer<Scalar> phisCur( numBf , numPts , dkcard );
00296           Phis_.getValues( phisCur , inputPoints , operatorType );
00297 
00298           for (int i=0;i<outputValues.dimension(0);i++) {
00299             for (int j=0;j<outputValues.dimension(1);j++) {
00300               for (int k=0;k<outputValues.dimension(2);k++) {
00301                 outputValues(i,j,k) = 0.0;
00302                 for (int l=0;l<this->getCardinality();l++) {
00303                   outputValues(i,j,k) += this->Vinv_(l,i) * phisCur(l,j,k);
00304                 }
00305               }
00306             }
00307           }
00308         }
00309         break;
00310       default:
00311         TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument,
00312                             ">>> ERROR (Basis_HGRAD_LINE_Cn_FEM): Operator type not implemented" );
00313         break;
00314       }
00315     }
00316     catch (std::invalid_argument &exception){
00317       TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument,
00318                           ">>> ERROR (Basis_HGRAD_LINE_Cn_FEM): Operator failed");    
00319     }
00320 
00321   }
00322   
00323 
00324   
00325   template<class Scalar, class ArrayScalar>
00326   void Basis_HGRAD_LINE_Cn_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar&           outputValues,
00327                                                               const ArrayScalar &    inputPoints,
00328                                                               const ArrayScalar &    cellVertices,
00329                                                               const EOperator        operatorType) const {
00330     TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
00331                         ">>> ERROR (Basis_HGRAD_LINE_Cn_FEM): FEM Basis calling an FVD member function");
00332   }
00333 
00334 
00335   template<class Scalar, class ArrayScalar>
00336   void Basis_HGRAD_LINE_Cn_FEM<Scalar,ArrayScalar>::getDofCoords( ArrayScalar & dofCoords ) const
00337   {
00338     for (int i=0;i<latticePts_.dimension(0);i++)
00339       {
00340         for (int j=0;j<latticePts_.dimension(1);j++)
00341           {
00342             dofCoords(i,j) = latticePts_(i,j);
00343           }
00344       }
00345     return;
00346   }
00347 
00348 }// namespace Intrepid
00349 #endif