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