Intrepid
http://trilinos.sandia.gov/packages/docs/r10.6/packages/intrepid/src/Discretization/Basis/Intrepid_HCURL_TET_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_HCURL_TET_In_FEM<Scalar,ArrayScalar>::Basis_HCURL_TET_In_FEM( const int n ,
00041                                                                       const EPointType pointType ):
00042     Phis_( n ),
00043     coeffs_( (n+1)*(n+2)*(n+3)/2 , n*(n+2)*(n+3)/2  )
00044   {
00045      const int N = n*(n+2)*(n+3)/2;
00046      this -> basisCardinality_  = N;
00047      this -> basisDegree_       = n;
00048      this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<4> >() );
00049      this -> basisType_         = BASIS_FEM_FIAT;
00050      this -> basisCoordinates_  = COORDINATES_CARTESIAN;
00051      this -> basisTagsAreSet_   = false;
00052 
00053      const int littleN = n*(n+1)*(n+2)/2;    // dim of (P_{n-1})^3 -- smaller space
00054      const int bigN = (n+1)*(n+2)*(n+3)/2;   // dim of (P_{n})^3 -- larger space
00055      const int start_PkH = (n-1)*n*(n+1)/6;  // dim of P({n-2}), offset into 
00056      const int dim_PkH = n*(n+1)*(n+2)/6 - start_PkH;
00057      const int scalarLittleN = littleN/3;
00058      const int scalarBigN = bigN/3;
00059 
00060      // first, need to project the basis for Nedelec space onto the
00061      // orthogonal basis of degree n
00062      // get coefficients of PkHx
00063 
00064      Teuchos::SerialDenseMatrix<int,Scalar> V1(bigN, littleN + 3 * dim_PkH);
00065 
00066      // these two loops get the first three sets of basis functions
00067      for (int i=0;i<scalarLittleN;i++) {
00068        for (int k=0;k<3;k++) {
00069          V1(i+k*scalarBigN,i+k*scalarLittleN) = 1.0;
00070        }
00071      }
00072 
00073      // first 3*scalarLittleN columns are (P_{n-1})^3 space
00074 
00075 
00076      // now I need to integrate { (x,y,z) \times } against the big basis
00077      // first, get a cubature rule.
00078      CubatureDirectTetDefault<Scalar,FieldContainer<Scalar> > myCub( 2 * n );
00079      FieldContainer<Scalar> cubPoints( myCub.getNumPoints() , 3 );
00080      FieldContainer<Scalar> cubWeights( myCub.getNumPoints() );
00081      myCub.getCubature( cubPoints , cubWeights );
00082 
00083      // tabulate the scalar orthonormal basis at cubature points
00084      FieldContainer<Scalar> phisAtCubPoints( scalarBigN , myCub.getNumPoints() );
00085      Phis_.getValues( phisAtCubPoints , cubPoints , OPERATOR_VALUE );
00086 
00087 
00088 
00089      // first set of these functions will write into the first dimPkH columns of remainder
00090 
00091      for (int j=0;j<dim_PkH;j++) { // loop over homogeneous polynomials
00092        // write into second spatial component, where
00093        // I integrate z phi_j phi_i
00094        for (int i=0;i<scalarBigN;i++) {
00095          V1(scalarBigN+i,littleN+j) = 0.0;
00096          for (int k=0;k<myCub.getNumPoints();k++) {
00097            V1(scalarBigN+i,littleN+j) -= cubWeights(k) * cubPoints(k,2) 
00098              * phisAtCubPoints(start_PkH+j,k)
00099              * phisAtCubPoints(i,k);
00100          }
00101        }  
00102        // write into third spatial component (-y phi_j, phi_i)
00103        for (int i=0;i<scalarBigN;i++) {
00104          V1(2*scalarBigN+i,littleN+j) = 0.0;
00105          for (int k=0;k<myCub.getNumPoints();k++) {
00106            V1(2*scalarBigN+i,littleN+j) += cubWeights(k) * cubPoints(k,1) 
00107              * phisAtCubPoints(start_PkH+j,k)
00108              * phisAtCubPoints(i,k);
00109          }
00110        }  
00111      }
00112 
00113      // second set of basis functions, write into second set of dimPkH columns
00114      for (int j=0;j<dim_PkH;j++) { // loop over homogeneous polynomials
00115        // write into first spatial component, where
00116        // I integrate -z phi_j phi_i
00117        for (int i=0;i<scalarBigN;i++) {
00118          V1(i,littleN+dim_PkH+j) = 0.0;
00119          for (int k=0;k<myCub.getNumPoints();k++) {
00120            V1(i,littleN+dim_PkH+j) += cubWeights(k) * cubPoints(k,2) 
00121              * phisAtCubPoints(start_PkH+j,k)
00122              * phisAtCubPoints(i,k);
00123          }
00124        } 
00125      
00126        // third spatial component, x phi_j phi_i
00127        for (int i=0;i<scalarBigN;i++) {
00128          V1(2*scalarBigN+i,littleN+dim_PkH+j) = 0.0;
00129          for (int k=0;k<myCub.getNumPoints();k++) {
00130            V1(2*scalarBigN+i,littleN+dim_PkH+j) -= cubWeights(k) * cubPoints(k,0) 
00131              * phisAtCubPoints(start_PkH+j,k)
00132              * phisAtCubPoints(i,k);
00133          }
00134        }  
00135      }
00136     
00137      // third clump of dimPkH columns
00138      for (int j=0;j<dim_PkH;j++) { // loop over homogeneous polynomials
00139        // write into first spatial component, where
00140        // I integrate y phi_j phi_i
00141        for (int i=0;i<scalarBigN;i++) {
00142          V1(i,littleN+2*dim_PkH+j) = 0.0;
00143          for (int k=0;k<myCub.getNumPoints();k++) {
00144            V1(i,littleN+2*dim_PkH+j) -= cubWeights(k) * cubPoints(k,1) 
00145              * phisAtCubPoints(start_PkH+j,k)
00146              * phisAtCubPoints(i,k);
00147          }
00148        }  
00149        // second spatial component, -x phi_j phi_i
00150        for (int i=0;i<scalarBigN;i++) {
00151          V1(scalarBigN+i,littleN+2*dim_PkH+j) = 0.0;
00152          for (int k=0;k<myCub.getNumPoints();k++) {
00153            V1(scalarBigN+i,littleN+2*dim_PkH+j) += cubWeights(k) * cubPoints(k,0) 
00154              * phisAtCubPoints(start_PkH+j,k)
00155              * phisAtCubPoints(i,k);
00156          }
00157        }  
00158      }
00159 
00160      // now I need to set up an SVD to get a basis for the space
00161      Teuchos::SerialDenseMatrix<int,Scalar> S(bigN,1);
00162      Teuchos::SerialDenseMatrix<int,Scalar> U(bigN, bigN);
00163      Teuchos::SerialDenseMatrix<int,Scalar> Vt(bigN,bigN);
00164      Teuchos::SerialDenseMatrix<int,Scalar> work(5*bigN,1);
00165      Teuchos::SerialDenseMatrix<int,Scalar> rWork(1,1);
00166      int info;
00167 
00168      Teuchos::LAPACK<int,Scalar> lala;
00169 
00170      lala.GESVD( 'A',  
00171                  'N',
00172                  V1.numRows() ,
00173                  V1.numCols() ,
00174                  V1.values() ,
00175                  V1.stride() ,
00176                  S.values() ,
00177                  U.values() ,
00178                  U.stride() ,
00179                  Vt.values() ,
00180                  Vt.stride() ,
00181                  work.values() ,
00182                  5*bigN ,
00183                  rWork.values() ,
00184                  &info );
00185                         
00186      int num_nonzero_sv = 0;
00187      for (int i=0;i<bigN;i++) {
00188        if (S(i,0) > INTREPID_TOL) {
00189          num_nonzero_sv++;
00190        }
00191      }
00192      
00193      Teuchos::SerialDenseMatrix<int,Scalar> Uslender(bigN, num_nonzero_sv);
00194      for (int j=0;j<num_nonzero_sv;j++) {
00195        for (int i=0;i<bigN;i++) {
00196          Uslender(i,j) = U(i,j);
00197        }
00198      }
00199 
00200      // apply nodes to big space
00201      Teuchos::SerialDenseMatrix<int,Scalar> V2(N, bigN);
00202 
00203      shards::CellTopology edgeTop(shards::getCellTopologyData<shards::Line<2> >() );
00204      shards::CellTopology faceTop(shards::getCellTopologyData<shards::Triangle<3> >() );
00205 
00206 
00207      const int numPtsPerEdge = PointTools::getLatticeSize( edgeTop ,
00208                                                            n+1 ,
00209                                                            1 );
00210 
00211      const int numPtsPerFace = PointTools::getLatticeSize( faceTop ,
00212                                                            n+1 ,
00213                                                            1 );
00214 
00215      const int numPtsPerCell = PointTools::getLatticeSize( this->basisCellTopology_ ,
00216                                                            n+1 ,
00217                                                            1 );
00218     
00219      // these hold the reference domain points that will be mapped to each edge or face
00220      FieldContainer<Scalar> oneDPts( numPtsPerEdge , 1 );
00221      FieldContainer<Scalar> twoDPts( numPtsPerFace , 2 );
00222 
00223      if (pointType == POINTTYPE_WARPBLEND) {
00224        CubatureDirectLineGauss<Scalar> edgeRule( numPtsPerEdge );
00225        FieldContainer<Scalar> edgeCubWts( numPtsPerEdge );
00226        edgeRule.getCubature( oneDPts , edgeCubWts );
00227      }
00228      else if (pointType == POINTTYPE_EQUISPACED ) {
00229        PointTools::getLattice<Scalar,FieldContainer<Scalar> >( oneDPts , 
00230                                                                edgeTop ,
00231                                                                n+1 , 
00232                                                                1 ,
00233                                                                pointType );
00234      }
00235 
00236      PointTools::getLattice<Scalar,FieldContainer<Scalar> >( twoDPts ,
00237                                                              faceTop ,
00238                                                              n+1 ,
00239                                                              1 ,
00240                                                              pointType );
00241 
00242      FieldContainer<Scalar> edgePts( numPtsPerEdge , 3 );
00243      FieldContainer<Scalar> phisAtEdgePoints( scalarBigN , numPtsPerEdge );
00244 
00245      FieldContainer<Scalar> facePts( numPtsPerFace , 3 );
00246      FieldContainer<Scalar> phisAtFacePoints( scalarBigN , 
00247                                              numPtsPerFace );   
00248 
00249      FieldContainer<Scalar> edgeTan( 3 );
00250     
00251      // loop over the edges
00252      for (int edge=0;edge<6;edge++) {
00253        CellTools<Scalar>::getReferenceEdgeTangent( edgeTan ,
00254                                                    edge ,
00255                                                    this->basisCellTopology_ );
00256        /* multiply by 2.0 to account for a scaling in Pavel's definition */
00257        for (int j=0;j<3;j++) {
00258          edgeTan(j) *= 2.0;
00259        }
00260 
00261        CellTools<Scalar>::mapToReferenceSubcell( edgePts ,
00262                                                  oneDPts ,
00263                                                  1 ,
00264                                                  edge ,
00265                                                  this->basisCellTopology_ );
00266       
00267        Phis_.getValues( phisAtEdgePoints , edgePts , OPERATOR_VALUE );
00268    
00269        // loop over points (rows of V2)
00270        for (int j=0;j<numPtsPerEdge;j++) {
00271          // loop over orthonormal basis functions (columns of V2)
00272          for (int k=0;k<scalarBigN;k++) {
00273            for (int d=0;d<3;d++) {
00274              V2(edge*numPtsPerEdge+j,k+scalarBigN*d) = edgeTan(d) * phisAtEdgePoints(k,j);
00275            }
00276          }
00277        }
00278      }
00279 
00280    // handle the faces, if needed
00281     if (n > 1) {
00282       FieldContainer<Scalar> refFaceTanU(3);
00283       FieldContainer<Scalar> refFaceTanV(3);
00284       for (int face=0;face<4;face++) {
00285         CellTools<Scalar>::getReferenceFaceTangents( refFaceTanU ,
00286                                                     refFaceTanV ,
00287                                                     face ,
00288                                                     this->basisCellTopology_ );
00289         CellTools<Scalar>::mapToReferenceSubcell( facePts ,
00290                                                   twoDPts ,
00291                                                   2 ,
00292                                                   face ,
00293                                                   this->basisCellTopology_ );
00294         Phis_.getValues( phisAtFacePoints , facePts , OPERATOR_VALUE );
00295         for (int j=0;j<numPtsPerFace;j++) {
00296           for (int k=0;k<scalarBigN;k++) {
00297             for (int d=0;d<3;d++) {
00298               V2(6*numPtsPerEdge+2*face*numPtsPerFace+2*j,k+scalarBigN*d) =
00299                 refFaceTanU(d) * phisAtFacePoints(k,j);
00300               V2(6*numPtsPerEdge+2*face*numPtsPerFace+2*j+1,k+scalarBigN*d) =
00301                 refFaceTanV(d) * phisAtFacePoints(k,j);
00302             }
00303           }
00304         }
00305       }
00306    }
00307 
00308      // internal dof, if needed
00309      if (n > 2) {
00310        FieldContainer<Scalar> cellPoints( numPtsPerCell , 3 );
00311        PointTools::getLattice<Scalar,FieldContainer<Scalar> >( cellPoints ,
00312                                                                this->getBaseCellTopology() , 
00313                                                                n + 1 ,
00314                                                                1 ,
00315                                                                pointType );
00316        FieldContainer<Scalar> phisAtCellPoints( scalarBigN , numPtsPerCell );
00317        Phis_.getValues( phisAtCellPoints , cellPoints , OPERATOR_VALUE );
00318        for (int i=0;i<numPtsPerCell;i++) {
00319          for (int j=0;j<scalarBigN;j++) {
00320            for (int k=0;k<3;k++) {
00321              V2(6*numPtsPerEdge+8*numPtsPerFace+k*numPtsPerCell+i,k*scalarBigN+j) = phisAtCellPoints(j,i);
00322            }
00323          }
00324        }
00325      }
00326 
00327     Teuchos::SerialDenseMatrix<int,Scalar> Vsdm( N , N );
00328     
00329     // multiply V2 * U --> V
00330     Vsdm.multiply( Teuchos::NO_TRANS , Teuchos::NO_TRANS , 1.0 , V2 , Uslender , 0.0 );
00331 
00332     Teuchos::SerialDenseSolver<int,Scalar> solver;
00333     solver.setMatrix( rcp( &Vsdm , false ) );
00334 
00335     solver.invert( );
00336 
00337 
00338     Teuchos::SerialDenseMatrix<int,Scalar> Csdm( bigN , N );
00339     Csdm.multiply( Teuchos::NO_TRANS , Teuchos::NO_TRANS , 1.0 , Uslender , Vsdm , 0.0 );
00340 
00341     //std::cout << Csdm << "\n";
00342 
00343     for (int i=0;i<bigN;i++) {
00344       for (int j=0;j<N;j++) {
00345         coeffs_(i,j) = Csdm(i,j);
00346       }
00347     }
00348 
00349     //std::cout << coeffs_ << std::endl;
00350 
00351   }
00352     
00353   template<class Scalar, class ArrayScalar>
00354   void Basis_HCURL_TET_In_FEM<Scalar, ArrayScalar>::initializeTags() {
00355     // Basis-dependent initializations
00356     int tagSize  = 4;        // size of DoF tag, i.e., number of fields in the tag
00357     int posScDim = 0;        // position in the tag, counting from 0, of the subcell dim 
00358     int posScOrd = 1;        // position in the tag, counting from 0, of the subcell ordinal
00359     int posDfOrd = 2;        // position in the tag, counting from 0, of DoF ordinal relative to the subcell
00360     
00361     // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration 
00362     
00363     int *tags = new int[ tagSize * this->getCardinality() ];
00364     int *tag_cur = tags;
00365     const int deg = this->getDegree();
00366 
00367     shards::CellTopology edgeTop(shards::getCellTopologyData<shards::Line<2> >() );
00368     shards::CellTopology faceTop(shards::getCellTopologyData<shards::Triangle<3> >() );
00369 
00370 
00371     const int numPtsPerEdge = PointTools::getLatticeSize( edgeTop ,
00372                                                           deg+1 ,
00373                                                           1 );
00374 
00375     const int numPtsPerFace = PointTools::getLatticeSize( faceTop ,
00376                                                           deg+1 ,
00377                                                           1 );
00378 
00379     const int numPtsPerCell = PointTools::getLatticeSize( this->basisCellTopology_ ,
00380                                                           deg+1 ,
00381                                                           1 );
00382 
00383     // edge dof first
00384     for (int e=0;e<6;e++) {
00385       for (int i=0;i<numPtsPerEdge;i++) {
00386         tag_cur[0] = 1;  tag_cur[1] = e;  tag_cur[2] = i;  tag_cur[3] = numPtsPerEdge;
00387         tag_cur += tagSize;
00388       }
00389     }
00390 
00391     // face dof, 2 * numPtsPerFace dof per face
00392     for (int f=0;f<4;f++) {
00393       for (int i=0;i<2*numPtsPerFace;i++) {
00394         tag_cur[0] = 2;  tag_cur[1] = f;  tag_cur[2] = i;  tag_cur[3] = 2*numPtsPerFace;
00395         tag_cur+= tagSize;
00396       }
00397     }
00398     
00399     // internal dof, 3 * numPtsPerCell
00400     for (int i=0;i<3*numPtsPerCell;i++) {
00401       tag_cur[0] = 3;  tag_cur[1] = 0;  tag_cur[2] = i;  tag_cur[3] = 3*numPtsPerCell;
00402       tag_cur += tagSize;
00403     }
00404     Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
00405                                 this -> ordinalToTag_,
00406                                 tags,
00407                                 this -> basisCardinality_,
00408                                 tagSize,
00409                                 posScDim,
00410                                 posScOrd,
00411                                 posDfOrd);
00412 
00413     delete []tags;   
00414     
00415   }  
00416 
00417 
00418 
00419   template<class Scalar, class ArrayScalar> 
00420   void Basis_HCURL_TET_In_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar &        outputValues,
00421                                                               const ArrayScalar &  inputPoints,
00422                                                               const EOperator      operatorType) const {
00423   
00424     // Verify arguments
00425 #ifdef HAVE_INTREPID_DEBUG
00426     Intrepid::getValues_HCURL_Args<Scalar, ArrayScalar>(outputValues,
00427                                                         inputPoints,
00428                                                         operatorType,
00429                                                         this -> getBaseCellTopology(),
00430                                                         this -> getCardinality() );
00431 #endif
00432     const int numPts = inputPoints.dimension(0);
00433     const int deg = this -> getDegree();
00434     const int scalarBigN = (deg+1)*(deg+2)*(deg+3)/6;
00435 
00436     try {
00437       switch (operatorType) {
00438       case OPERATOR_VALUE:
00439         {
00440           FieldContainer<Scalar> phisCur( scalarBigN , numPts );
00441           Phis_.getValues( phisCur , inputPoints , OPERATOR_VALUE );
00442 
00443           for (int i=0;i<outputValues.dimension(0);i++) { // RT bf
00444             for (int j=0;j<outputValues.dimension(1);j++) {  // point
00445               for (int d=0;d<3;d++) {
00446                 outputValues(i,j,d) = 0.0;
00447               }
00448               for (int k=0;k<scalarBigN;k++) { // Dubiner bf
00449                 for (int d=0;d<3;d++) {
00450                   outputValues(i,j,d) += coeffs_(k+d*scalarBigN,i) * phisCur(k,j);
00451                 }
00452               }
00453             }
00454           }
00455         }
00456         break;
00457       case OPERATOR_CURL:
00458         {
00459           FieldContainer<Scalar> phisCur( scalarBigN , numPts , 3 );
00460           Phis_.getValues( phisCur , inputPoints , OPERATOR_GRAD );
00461           for (int i=0;i<outputValues.dimension(0);i++) { // bf loop
00462             for (int j=0;j<outputValues.dimension(1);j++) { // point loop
00463               outputValues(i,j,0) = 0.0;
00464               for (int k=0;k<scalarBigN;k++) {
00465                 outputValues(i,j,0) += coeffs_(k+2*scalarBigN,i) * phisCur(k,j,1);
00466                 outputValues(i,j,0) -= coeffs_(k+scalarBigN,i) * phisCur(k,j,2);
00467               }
00468               
00469               outputValues(i,j,1) = 0.0;
00470               for (int k=0;k<scalarBigN;k++) {
00471                 outputValues(i,j,1) += coeffs_(k,i) * phisCur(k,j,2);
00472                 outputValues(i,j,1) -= coeffs_(k+2*scalarBigN,i) * phisCur(k,j,0);
00473               }
00474 
00475               outputValues(i,j,2) = 0.0;
00476               for (int k=0;k<scalarBigN;k++) {
00477                 outputValues(i,j,2) += coeffs_(k+scalarBigN,i) * phisCur(k,j,0);
00478                 outputValues(i,j,2) -= coeffs_(k,i) * phisCur(k,j,1);
00479               }
00480             }
00481           }
00482         }
00483         break;
00484       default:
00485         TEST_FOR_EXCEPTION( true , std::invalid_argument,
00486                             ">>> ERROR (Basis_HCURL_TET_In_FEM): Operator type not implemented");
00487         break;
00488       }
00489     }
00490     catch (std::invalid_argument &exception){
00491       TEST_FOR_EXCEPTION( true , std::invalid_argument,
00492                           ">>> ERROR (Basis_HCURL_TET_In_FEM): Operator type not implemented");    
00493     }
00494 
00495   }
00496   
00497 
00498   
00499   template<class Scalar, class ArrayScalar>
00500   void Basis_HCURL_TET_In_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar&           outputValues,
00501                                                               const ArrayScalar &    inputPoints,
00502                                                               const ArrayScalar &    cellVertices,
00503                                                               const EOperator        operatorType) const {
00504     TEST_FOR_EXCEPTION( (true), std::logic_error,
00505                         ">>> ERROR (Basis_HCURL_TET_In_FEM): FEM Basis calling an FVD member function");
00506   }
00507 
00508 
00509 }// namespace Intrepid