Intrepid
http://trilinos.sandia.gov/packages/docs/r10.6/packages/intrepid/src/Discretization/Basis/Intrepid_HDIV_TRI_In_FEMDef.hpp
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 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Pavel Bochev (pbboche@sandia.gov) or
00025 //                    Denis Ridzal (dridzal@sandia.gov)
00026 //                    Kara Peterson (kjpeter@sandia.gov).
00027 //
00028 // ************************************************************************
00029 // @HEADER
00030 
00037 namespace Intrepid {
00038 
00039   template<class Scalar, class ArrayScalar>
00040   Basis_HDIV_TRI_In_FEM<Scalar,ArrayScalar>::Basis_HDIV_TRI_In_FEM( const int n ,
00041                                                                     const EPointType pointType ):
00042     Phis( n ),
00043     coeffs( (n+1)*(n+2) , n*(n+2) )
00044   {
00045     const int N = n*(n+2);
00046     this -> basisCardinality_  = N;
00047     this -> basisDegree_       = n;
00048     this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Triangle<3> >() );
00049     this -> basisType_         = BASIS_FEM_FIAT;
00050     this -> basisCoordinates_  = COORDINATES_CARTESIAN;
00051     this -> basisTagsAreSet_   = false;
00052 
00053 
00054     const int littleN = n*(n+1);   // dim of (P_{n-1})^2 -- smaller space
00055     const int bigN = (n+1)*(n+2);  // dim of (P_{n})^2 -- larger space
00056     const int scalarSmallestN = (n-1)*n / 2;
00057     const int scalarLittleN = littleN/2;
00058     const int scalarBigN = bigN/2;
00059 
00060     // first, need to project the basis for RT space onto the
00061     // orthogonal basis of degree n
00062     // get coefficients of PkHx
00063 
00064     Teuchos::SerialDenseMatrix<int,Scalar> V1(bigN, N);
00065 
00066     // basis for the space is 
00067     // { (phi_i,0) }_{i=0}^{scalarLittleN-1} ,
00068     // { (0,phi_i) }_{i=0}^{scalarLittleN-1} ,
00069     // { (x,y) . phi_i}_{i=scalarLittleN}^{scalarBigN-1}
00070     // columns of V1 are expansion of this basis in terms of the basis
00071     // for P_{n}^2
00072 
00073     // these two loops get the first two sets of basis functions
00074     for (int i=0;i<scalarLittleN;i++) {
00075       V1(i,i) = 1.0;
00076       V1(scalarBigN+i,scalarLittleN+i) = 1.0;
00077     }
00078 
00079     // now I need to integrate { (x,y) phi } against the big basis
00080     // first, get a cubature rule.
00081     CubatureDirectTriDefault<Scalar,FieldContainer<Scalar> > myCub( 2 * n );
00082     FieldContainer<Scalar> cubPoints( myCub.getNumPoints() , 2 );
00083     FieldContainer<Scalar> cubWeights( myCub.getNumPoints() );
00084     myCub.getCubature( cubPoints , cubWeights );
00085 
00086     // tabulate the scalar orthonormal basis at cubature points
00087     FieldContainer<Scalar> phisAtCubPoints( scalarBigN , myCub.getNumPoints() );
00088     Phis.getValues( phisAtCubPoints , cubPoints , OPERATOR_VALUE );
00089 
00090     // now do the integration
00091     for (int i=0;i<n;i++) {
00092       for (int j=0;j<scalarBigN;j++) { // int (x,y) phi_i \cdot (phi_j,0)
00093         V1(j,littleN+i) = 0.0;
00094         for (int k=0;k<myCub.getNumPoints();k++) {
00095           V1(j,littleN+i) += 
00096             cubWeights(k) * cubPoints(k,0) 
00097             * phisAtCubPoints(scalarSmallestN+i,k) 
00098             * phisAtCubPoints(j,k);
00099         }
00100       }
00101       for (int j=0;j<scalarBigN;j++) {  // int (x,y) phi_i \cdot (0,phi_j)
00102         V1(j+scalarBigN,littleN+i) = 0.0;
00103         for (int k=0;k<myCub.getNumPoints();k++) {
00104           V1(j+scalarBigN,littleN+i) += 
00105             cubWeights(k) * cubPoints(k,1) 
00106             * phisAtCubPoints(scalarSmallestN+i,k) 
00107             * phisAtCubPoints(j,k);
00108         }
00109       }
00110     }
00111 
00112     //std::cout << V1 << "\n";
00113 
00114     
00115     // next, apply the RT nodes (rows) to the basis for (P_n)^2 (columns)
00116     Teuchos::SerialDenseMatrix<int,Scalar> V2(N , bigN);
00117 
00118     // first 3 * degree nodes are normals at each edge
00119     // get the points on the line
00120     FieldContainer<Scalar> linePts( n , 1 );
00121     if (pointType == POINTTYPE_WARPBLEND) {
00122       CubatureDirectLineGauss<Scalar> edgeRule( n );
00123       FieldContainer<Scalar> edgeCubWts( n );
00124       edgeRule.getCubature( linePts , edgeCubWts );
00125     }
00126     else if (pointType == POINTTYPE_EQUISPACED ) {
00127       shards::CellTopology linetop(shards::getCellTopologyData<shards::Line<2> >() );
00128 
00129       PointTools::getLattice<Scalar,FieldContainer<Scalar> >( linePts , 
00130                                                               linetop ,
00131                                                               n+1 , 1 ,
00132                                                               POINTTYPE_EQUISPACED );
00133     }
00134     // holds the image of the line points 
00135     FieldContainer<Scalar> edgePts( n , 2 );
00136     FieldContainer<Scalar> phisAtEdgePoints( scalarBigN , n );
00137 
00138     // these are scaled by the appropriate edge lengths.
00139     const Scalar nx[] = {0.0,1.0,-1.0};
00140     const Scalar ny[] = {-1.0,1.0,0.0};
00141     
00142     for (int i=0;i<3;i++) {  // loop over edges
00143       CellTools<Scalar>::mapToReferenceSubcell( edgePts ,
00144                                                 linePts ,
00145                                                 1 ,
00146                                                 i ,
00147                                                 this->basisCellTopology_ );
00148 
00149       Phis.getValues( phisAtEdgePoints , edgePts , OPERATOR_VALUE );
00150 
00151       // loop over points (rows of V2)
00152       for (int j=0;j<n;j++) {
00153         // loop over orthonormal basis functions (columns of V2)
00154         for (int k=0;k<scalarBigN;k++) {
00155           V2(n*i+j,k) = nx[i] * phisAtEdgePoints(k,j);
00156           V2(n*i+j,k+scalarBigN) = ny[i] * phisAtEdgePoints(k,j);
00157         }
00158       }
00159     }
00160 
00161     // next map the points to each edge
00162 
00163 
00164     // remaining nodes are divided into two pieces:  point value of x
00165     // components and point values of y components.  These are
00166     // evaluated at the interior of a lattice of degree + 1, For then
00167     // the degree == 1 space corresponds classicaly to RT0 and so gets
00168     // no internal nodes, and degree == 2 corresponds to RT1 and needs
00169     // one internal node per vector component.
00170     const int numInternalPoints = PointTools::getLatticeSize( this->getBaseCellTopology() ,
00171                                                               n + 1 ,
00172                                                               1 );
00173 
00174     if (numInternalPoints > 0) {
00175       FieldContainer<Scalar> internalPoints( numInternalPoints , 2 );
00176       PointTools::getLattice<Scalar,FieldContainer<Scalar> >( internalPoints ,
00177                                                               this->getBaseCellTopology() , 
00178                                                               n + 1 ,
00179                                                               1 ,
00180                                                               pointType );
00181       
00182       FieldContainer<Scalar> phisAtInternalPoints( scalarBigN , numInternalPoints );
00183       Phis.getValues( phisAtInternalPoints , internalPoints , OPERATOR_VALUE );
00184       
00185       // copy values into right positions of V2
00186       for (int i=0;i<numInternalPoints;i++) {
00187         for (int j=0;j<scalarBigN;j++) {
00188           // x component
00189           V2(3*n+i,j) = phisAtInternalPoints(j,i);
00190           // y component
00191           V2(3*n+numInternalPoints+i,scalarBigN+j) = phisAtInternalPoints(j,i);
00192         }
00193       }
00194     }
00195 //     std::cout << "Nodes on big basis\n";
00196 //     std::cout << V2 << "\n";
00197 //     std::cout << "End nodes\n";
00198 
00199     Teuchos::SerialDenseMatrix<int,Scalar> Vsdm( N , N );
00200 
00201     // multiply V2 * V1 --> V
00202     Vsdm.multiply( Teuchos::NO_TRANS , Teuchos::NO_TRANS , 1.0 , V2 , V1 , 0.0 );
00203 
00204 //     std::cout << "Vandermonde:\n";
00205 //     std::cout << Vsdm << "\n";
00206 //     std::cout << "End Vandermonde\n";
00207     
00208     Teuchos::SerialDenseSolver<int,Scalar> solver;
00209     solver.setMatrix( rcp( &Vsdm , false ) );
00210     solver.invert( );
00211 
00212     Teuchos::SerialDenseMatrix<int,Scalar> Csdm( bigN , N );
00213     Csdm.multiply( Teuchos::NO_TRANS , Teuchos::NO_TRANS , 1.0 , V1 , Vsdm , 0.0 );
00214 
00215     //    std::cout << Csdm << "\n";
00216 
00217     for (int i=0;i<bigN;i++) {
00218       for (int j=0;j<N;j++) {
00219         coeffs(i,j) = Csdm(i,j);
00220       }
00221     }
00222   }  
00223     
00224   template<class Scalar, class ArrayScalar>
00225   void Basis_HDIV_TRI_In_FEM<Scalar, ArrayScalar>::initializeTags() {
00226   
00227     // Basis-dependent initializations
00228     int tagSize  = 4;        // size of DoF tag, i.e., number of fields in the tag
00229     int posScDim = 0;        // position in the tag, counting from 0, of the subcell dim 
00230     int posScOrd = 1;        // position in the tag, counting from 0, of the subcell ordinal
00231     int posDfOrd = 2;        // position in the tag, counting from 0, of DoF ordinal relative to the subcell
00232   
00233     // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration 
00234 
00235     int *tags = new int[ tagSize * this->getCardinality() ];
00236     int *tag_cur = tags;
00237     const int degree = this->getDegree();
00238 
00239     // there are degree internal dofs on each edge -- normals.  Let's do them
00240     for (int ed=0;ed<3;ed++) {
00241       for (int i=0;i<degree;i++) {
00242         tag_cur[0] = 1;  tag_cur[1] = ed;  tag_cur[2] = i;  tag_cur[3] = degree;
00243         tag_cur += tagSize;
00244       }
00245     }
00246                           
00247     // end edge dofs
00248 
00249     // the rest of the dofs are internal
00250     const int numFaceDof = (degree-1)*degree;
00251     int faceDofCur = 0;
00252     for (int i=3*degree;i<degree*(degree+1);i++) {
00253       tag_cur[0] = 2;  tag_cur[1] = 0;  tag_cur[2] = faceDofCur;  tag_cur[3] = numFaceDof;
00254       tag_cur += tagSize;
00255       faceDofCur++;
00256     }
00257     
00258     
00259     Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
00260                                 this -> ordinalToTag_,
00261                                 tags,
00262                                 this -> basisCardinality_,
00263                                 tagSize,
00264                                 posScDim,
00265                                 posScOrd,
00266                                 posDfOrd);
00267 
00268     delete []tags;
00269   
00270   }  
00271 
00272 
00273 
00274   template<class Scalar, class ArrayScalar> 
00275   void Basis_HDIV_TRI_In_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar &        outputValues,
00276                                                               const ArrayScalar &  inputPoints,
00277                                                               const EOperator      operatorType) const {
00278   
00279     // Verify arguments
00280 #ifdef HAVE_INTREPID_DEBUG
00281     Intrepid::getValues_HDIV_Args<Scalar, ArrayScalar>(outputValues,
00282                                                         inputPoints,
00283                                                         operatorType,
00284                                                         this -> getBaseCellTopology(),
00285                                                         this -> getCardinality() );
00286 #endif
00287     const int numPts = inputPoints.dimension(0);
00288     const int deg = this -> getDegree();
00289     const int scalarBigN = (deg+1)*(deg+2)/2;
00290 
00291     try {
00292       switch (operatorType) {
00293       case OPERATOR_VALUE:
00294         {
00295           FieldContainer<Scalar> phisCur( scalarBigN , numPts );
00296           Phis.getValues( phisCur , inputPoints , OPERATOR_VALUE );
00297 
00298           for (int i=0;i<outputValues.dimension(0);i++) { // RT bf
00299             for (int j=0;j<outputValues.dimension(1);j++) {  // point
00300               outputValues(i,j,0) = 0.0;
00301               outputValues(i,j,1) = 0.0;
00302               for (int k=0;k<scalarBigN;k++) { // Dubiner bf
00303                 outputValues(i,j,0) += coeffs(k,i) * phisCur(k,j);
00304                 outputValues(i,j,1) += coeffs(k+scalarBigN,i) * phisCur(k,j);
00305               }
00306             }
00307           }
00308         }
00309         break;
00310       case OPERATOR_DIV:
00311         {
00312           FieldContainer<Scalar> phisCur( scalarBigN , numPts , 2 );
00313           Phis.getValues( phisCur , inputPoints , OPERATOR_GRAD );
00314           for (int i=0;i<outputValues.dimension(0);i++) { // bf loop
00315             for (int j=0;j<outputValues.dimension(1);j++) { // point loop
00316               // dx of x component
00317               outputValues(i,j) = 0.0;
00318               for (int k=0;k<scalarBigN;k++) {
00319                 outputValues(i,j) += coeffs(k,i) * phisCur(k,j,0);
00320               }
00321               // dy of y component
00322               for (int k=0;k<scalarBigN;k++) {
00323                 outputValues(i,j) += coeffs(k+scalarBigN,i) * phisCur(k,j,1);
00324               }
00325             }
00326           }
00327         }
00328         break;
00329       default:
00330         TEST_FOR_EXCEPTION( true , std::invalid_argument,
00331                             ">>> ERROR (Basis_HDIV_TRI_In_FEM): Operator type not implemented");
00332         break;
00333       }
00334     }
00335     catch (std::invalid_argument &exception){
00336       TEST_FOR_EXCEPTION( true , std::invalid_argument,
00337                           ">>> ERROR (Basis_HDIV_TRI_In_FEM): Operator type not implemented");    
00338     }
00339 
00340   }
00341   
00342 
00343   
00344   template<class Scalar, class ArrayScalar>
00345   void Basis_HDIV_TRI_In_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar&           outputValues,
00346                                                               const ArrayScalar &    inputPoints,
00347                                                               const ArrayScalar &    cellVertices,
00348                                                               const EOperator        operatorType) const {
00349     TEST_FOR_EXCEPTION( (true), std::logic_error,
00350                         ">>> ERROR (Basis_HDIV_TRI_In_FEM): FEM Basis calling an FVD member function");
00351   }
00352 
00353 
00354 }// namespace Intrepid