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