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