Intrepid
http://trilinos.sandia.gov/packages/docs/r10.10/packages/intrepid/test/Discretization/Basis/HDIV_TET_In_FEM/test_02.cpp
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 
00049 #include "Intrepid_FieldContainer.hpp"
00050 #include "Intrepid_HGRAD_TET_Cn_FEM_ORTH.hpp"
00051 #include "Intrepid_HDIV_TET_In_FEM.hpp"
00052 #include "Intrepid_DefaultCubatureFactory.hpp"
00053 #include "Intrepid_RealSpaceTools.hpp"
00054 #include "Intrepid_ArrayTools.hpp"
00055 #include "Intrepid_FunctionSpaceTools.hpp"
00056 #include "Intrepid_CellTools.hpp"
00057 #include "Teuchos_oblackholestream.hpp"
00058 #include "Teuchos_RCP.hpp"
00059 #include "Teuchos_GlobalMPISession.hpp"
00060 #include "Teuchos_SerialDenseMatrix.hpp"
00061 #include "Teuchos_SerialDenseVector.hpp"
00062 #include "Teuchos_LAPACK.hpp"
00063 
00064 using namespace std;
00065 using namespace Intrepid;
00066 
00067 void rhsFunc( FieldContainer<double> &, const FieldContainer<double> &, int, int, int );
00068 void u_exact( FieldContainer<double> &, const FieldContainer<double> &, int, int, int );
00069 
00070 // This is the rhs for (div tau,w) = (f,w),
00071 // which makes f the negative Laplacian of scalar solution
00072 void rhsFunc( FieldContainer<double> &result, 
00073               const FieldContainer<double> &points,
00074               int xd,
00075               int yd ,
00076               int zd)
00077 {
00078   for (int cell=0;cell<result.dimension(0);cell++) {
00079     for (int pt=0;pt<result.dimension(1);pt++) {
00080       result(cell,pt) = 0.0;
00081       if (xd >=2) {
00082         result(cell,pt) += xd*(xd-1)*pow(points(cell,pt,0),xd-2)*pow(points(cell,pt,1),yd)
00083           *pow(points(cell,pt,2),zd);
00084       }
00085       if (yd >=2) {
00086         result(cell,pt) += yd*(yd-1)*pow(points(cell,pt,0),xd)*pow(points(cell,pt,1),yd-2)
00087           *pow(points(cell,pt,2),zd);
00088       }
00089       if (zd>=2) {
00090         result(cell,pt) += zd*(zd-1)*pow(points(cell,pt,0),xd)*pow(points(cell,pt,1),yd)
00091           *pow(points(cell,pt,2),zd-2);
00092         
00093       }
00094     }
00095   }
00096 }
00097 
00098 void u_exact( FieldContainer<double> &result, 
00099               const FieldContainer<double> &points,
00100               int xd,
00101               int yd,
00102               int zd)
00103 {
00104   for (int cell=0;cell<result.dimension(0);cell++){
00105     for (int pt=0;pt<result.dimension(1);pt++) {
00106       result(cell,pt) = std::pow(points(cell,pt,0),xd)*std::pow(points(cell,pt,1),yd)
00107         *std::pow(points(cell,pt,2),zd);
00108     }
00109   }
00110   return;
00111 }
00112 
00113 int main(int argc, char *argv[]) {
00114   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00115   
00116   // This little trick lets us print to std::cout only if
00117   // a (dummy) command-line argument is provided.
00118   int iprint = argc - 1;
00119   Teuchos::RCP<std::ostream> outStream;
00120   Teuchos::oblackholestream bhs; // outputs nothing
00121   if (iprint > 0)
00122     outStream = Teuchos::rcp(&std::cout, false);
00123   else
00124     outStream = Teuchos::rcp(&bhs, false);
00125   
00126   // Save the format state of the original std::cout.
00127   Teuchos::oblackholestream oldFormatState;
00128   oldFormatState.copyfmt(std::cout);
00129   
00130   *outStream \
00131     << "===============================================================================\n" \
00132     << "|                                                                             |\n" \
00133     << "|                  Unit Test (Basis_HGRAD_TET_In_FEM)                         |\n" \
00134     << "|                                                                             |\n" \
00135     << "| 1) Patch test involving H(div) matrices                                     |\n" \
00136     << "|    for the Dirichlet problem on a tetrahedral patch                         |\n" \
00137     << "|    Omega with boundary Gamma.                                               |\n" \
00138     << "|                                                                             |\n" \
00139     << "|   Questions? Contact Pavel Bochev (pbboche@sandia.gov),                     |\n" \
00140     << "|                      Robert Kirby (robert.c.kirby@ttu.edu),                 |\n" \
00141     << "|                      Denis Ridzal (dridzal@sandia.gov),                     |\n" \
00142     << "|                      Kara Peterson (kjpeter@sandia.gov).                    |\n" \
00143     << "|                                                                             |\n" \
00144     << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00145     << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00146     << "|                                                                             |\n" \
00147     << "===============================================================================\n" \
00148     << "| TEST 1: Patch test                                                          |\n" \
00149     << "===============================================================================\n";
00150   
00151   
00152   int errorFlag = 0;
00153   
00154   outStream -> precision(16);
00155   
00156   try {
00157     DefaultCubatureFactory<double> cubFactory;                                          // create cubature factory
00158     shards::CellTopology cell(shards::getCellTopologyData< shards::Tetrahedron<> >());  // create parent cell topology
00159     shards::CellTopology side(shards::getCellTopologyData< shards::Triangle<> >());     // create relevant subcell (side) topology
00160     
00161     int cellDim = cell.getDimension();
00162     int sideDim = side.getDimension();
00163     
00164     int min_order = 0;
00165     int max_order = 5;
00166     
00167     int numIntervals = max_order;
00168     int numInterpPoints = ((numIntervals + 1)*(numIntervals + 2)*(numIntervals+3))/6;
00169     FieldContainer<double> interp_points_ref(numInterpPoints, cellDim);
00170     int counter = 0;
00171     for (int j=0; j<=numIntervals; j++) {
00172       for (int i=0; i<=numIntervals-j; i++) {
00173         for (int k=0;k<numIntervals-j-i;k++) {
00174           interp_points_ref(counter,0) = i*(1.0/numIntervals);
00175           interp_points_ref(counter,1) = j*(1.0/numIntervals);
00176           interp_points_ref(counter,2) = k*(1.0/numIntervals);
00177           counter++;
00178         }
00179       }
00180     }
00181     
00182     //interp_points_ref.resize(numInterpPoints,cellDim);
00183     
00184     for (int basis_order=min_order;basis_order<=max_order;basis_order++) {
00185       // create bases
00186       Teuchos::RCP<Basis<double,FieldContainer<double> > > vectorBasis =
00187         Teuchos::rcp(new Basis_HDIV_TET_In_FEM<double,FieldContainer<double> >(basis_order+1,POINTTYPE_EQUISPACED) );
00188       Teuchos::RCP<Basis<double,FieldContainer<double> > > scalarBasis =
00189         Teuchos::rcp(new Basis_HGRAD_TET_Cn_FEM_ORTH<double,FieldContainer<double> >(basis_order) );
00190       
00191       int numVectorFields = vectorBasis->getCardinality();
00192       int numScalarFields = scalarBasis->getCardinality();
00193       int numTotalFields = numVectorFields + numScalarFields;
00194       
00195       // create cubatures
00196       Teuchos::RCP<Cubature<double> > cellCub = cubFactory.create(cell, 2*(basis_order+1));
00197       Teuchos::RCP<Cubature<double> > sideCub = cubFactory.create(side, 2*(basis_order+1));
00198       
00199       int numCubPointsCell = cellCub->getNumPoints();
00200       int numCubPointsSide = sideCub->getNumPoints();
00201       
00202       // hold cubature information
00203       FieldContainer<double> cub_points_cell(numCubPointsCell, cellDim);
00204       FieldContainer<double> cub_weights_cell(numCubPointsCell);
00205       FieldContainer<double> cub_points_side( numCubPointsSide, sideDim );
00206       FieldContainer<double> cub_weights_side( numCubPointsSide );
00207       FieldContainer<double> cub_points_side_refcell( numCubPointsSide , cellDim );
00208       
00209       // hold basis function information on refcell
00210       FieldContainer<double> value_of_v_basis_at_cub_points_cell(numVectorFields, numCubPointsCell, cellDim );
00211       FieldContainer<double> w_value_of_v_basis_at_cub_points_cell(1, numVectorFields, numCubPointsCell, cellDim);
00212       FieldContainer<double> div_of_v_basis_at_cub_points_cell( numVectorFields, numCubPointsCell );
00213       FieldContainer<double> w_div_of_v_basis_at_cub_points_cell( 1, numVectorFields , numCubPointsCell );
00214       FieldContainer<double> value_of_s_basis_at_cub_points_cell(numScalarFields,numCubPointsCell);
00215       FieldContainer<double> w_value_of_s_basis_at_cub_points_cell(1,numScalarFields,numCubPointsCell);
00216       
00217       // containers for side integration:
00218       // I just need the normal component of the vector basis
00219       // and the exact solution at the cub points
00220       FieldContainer<double> value_of_v_basis_at_cub_points_side(numVectorFields,numCubPointsSide,cellDim);
00221       FieldContainer<double> n_of_v_basis_at_cub_points_side(numVectorFields,numCubPointsSide);
00222       FieldContainer<double> w_n_of_v_basis_at_cub_points_side(1,numVectorFields,numCubPointsSide);
00223       FieldContainer<double> diri_data_at_cub_points_side(1,numCubPointsSide);
00224       FieldContainer<double> side_normal(cellDim);
00225       
00226       // holds rhs data
00227       FieldContainer<double> rhs_at_cub_points_cell(1,numCubPointsCell);
00228       
00229       // FEM matrices and vectors
00230       FieldContainer<double> fe_matrix_M(1,numVectorFields,numVectorFields);
00231       FieldContainer<double> fe_matrix_B(1,numVectorFields,numScalarFields);
00232       FieldContainer<double> fe_matrix(1,numTotalFields,numTotalFields);
00233       
00234       FieldContainer<double> rhs_vector_vec(1,numVectorFields);
00235       FieldContainer<double> rhs_vector_scal(1,numScalarFields);
00236       FieldContainer<double> rhs_and_soln_vec(1,numTotalFields);
00237       
00238       FieldContainer<int> ipiv(numTotalFields);
00239       FieldContainer<double> value_of_s_basis_at_interp_points( numScalarFields , numInterpPoints);
00240       FieldContainer<double> interpolant( 1 , numInterpPoints );
00241       
00242       // set test tolerance
00243       double zero = (basis_order+1)*(basis_order+1)*100*INTREPID_TOL;
00244       
00245       // build matrices outside the loop, and then just do the rhs
00246       // for each iteration
00247       
00248       cellCub->getCubature(cub_points_cell, cub_weights_cell);
00249       sideCub->getCubature(cub_points_side, cub_weights_side);
00250       
00251       // need the vector basis & its divergences
00252       vectorBasis->getValues(value_of_v_basis_at_cub_points_cell,
00253                              cub_points_cell,
00254                              OPERATOR_VALUE);
00255       vectorBasis->getValues(div_of_v_basis_at_cub_points_cell,
00256                              cub_points_cell,
00257                              OPERATOR_DIV);
00258       
00259       // need the scalar basis as well
00260       scalarBasis->getValues(value_of_s_basis_at_cub_points_cell,
00261                              cub_points_cell,
00262                              OPERATOR_VALUE);
00263       
00264       // construct mass matrix
00265       cub_weights_cell.resize(1,numCubPointsCell);
00266       FunctionSpaceTools::multiplyMeasure<double>(w_value_of_v_basis_at_cub_points_cell ,
00267                                                   cub_weights_cell ,
00268                                                   value_of_v_basis_at_cub_points_cell ); 
00269       cub_weights_cell.resize(numCubPointsCell);
00270       
00271       
00272       value_of_v_basis_at_cub_points_cell.resize( 1 , numVectorFields , numCubPointsCell , cellDim );
00273       FunctionSpaceTools::integrate<double>(fe_matrix_M,
00274                                             w_value_of_v_basis_at_cub_points_cell ,
00275                                             value_of_v_basis_at_cub_points_cell ,
00276                                             COMP_BLAS );
00277       value_of_v_basis_at_cub_points_cell.resize( numVectorFields , numCubPointsCell , cellDim );
00278       
00279       // div matrix
00280       cub_weights_cell.resize(1,numCubPointsCell);
00281       FunctionSpaceTools::multiplyMeasure<double>(w_div_of_v_basis_at_cub_points_cell,
00282                                                   cub_weights_cell,
00283                                                   div_of_v_basis_at_cub_points_cell);
00284       cub_weights_cell.resize(numCubPointsCell);
00285       
00286       value_of_s_basis_at_cub_points_cell.resize(1,numScalarFields,numCubPointsCell);
00287       FunctionSpaceTools::integrate<double>(fe_matrix_B,
00288                                             w_div_of_v_basis_at_cub_points_cell ,
00289                                             value_of_s_basis_at_cub_points_cell ,
00290                                             COMP_BLAS );
00291       value_of_s_basis_at_cub_points_cell.resize(numScalarFields,numCubPointsCell);
00292       
00293       
00294       // construct div matrix
00295       
00296       for (int x_order=0;x_order<=basis_order;x_order++) {
00297         for (int y_order=0;y_order<=basis_order-x_order;y_order++) {
00298           for (int z_order=0;z_order<=basis_order-x_order-y_order;z_order++) {
00299             
00300             
00301             // reset global matrix since I destroyed it in LU factorization.
00302             fe_matrix.initialize();
00303             // insert mass matrix into global matrix
00304             for (int i=0;i<numVectorFields;i++) {
00305               for (int j=0;j<numVectorFields;j++) {
00306                 fe_matrix(0,i,j) = fe_matrix_M(0,i,j);
00307               }
00308             }
00309             
00310             // insert div matrix into global matrix
00311             for (int i=0;i<numVectorFields;i++) {
00312               for (int j=0;j<numScalarFields;j++) {
00313                 fe_matrix(0,i,numVectorFields+j)=-fe_matrix_B(0,i,j);
00314                 fe_matrix(0,j+numVectorFields,i)=fe_matrix_B(0,i,j);
00315               }
00316             }
00317             
00318             // clear old vector data
00319             rhs_vector_vec.initialize();
00320             rhs_vector_scal.initialize();
00321             rhs_and_soln_vec.initialize();
00322             
00323             // now get rhs vector
00324             // rhs_vector_scal is just (rhs,w) for w in the scalar basis
00325             // I already have the scalar basis tabulated.
00326             cub_points_cell.resize(1,numCubPointsCell,cellDim);
00327             rhsFunc(rhs_at_cub_points_cell,
00328                     cub_points_cell,
00329                     x_order,
00330                     y_order,
00331                     z_order);
00332             
00333             cub_points_cell.resize(numCubPointsCell,cellDim);
00334             
00335             cub_weights_cell.resize(1,numCubPointsCell);
00336             FunctionSpaceTools::multiplyMeasure<double>(w_value_of_s_basis_at_cub_points_cell,
00337                                                         cub_weights_cell,
00338                                                         value_of_s_basis_at_cub_points_cell);
00339             cub_weights_cell.resize(numCubPointsCell);
00340             FunctionSpaceTools::integrate<double>(rhs_vector_scal,
00341                                                   rhs_at_cub_points_cell,
00342                                                   w_value_of_s_basis_at_cub_points_cell,
00343                                                   COMP_BLAS);
00344             
00345             for (int i=0;i<numScalarFields;i++) {
00346               rhs_and_soln_vec(0,numVectorFields+i) = rhs_vector_scal(0,i);
00347             }
00348             
00349             
00350             // now get <u,v.n> on boundary
00351             for (unsigned side_cur=0;side_cur<4;side_cur++) {
00352               // map side cubature to current side
00353               CellTools<double>::mapToReferenceSubcell( cub_points_side_refcell ,
00354                                                         cub_points_side ,
00355                                                         sideDim ,
00356                                                         (int)side_cur ,
00357                                                         cell );
00358               
00359               // Evaluate dirichlet data
00360               cub_points_side_refcell.resize(1,numCubPointsSide,cellDim);
00361               u_exact(diri_data_at_cub_points_side,
00362                       cub_points_side_refcell,x_order,y_order,z_order);
00363               cub_points_side_refcell.resize(numCubPointsSide,cellDim);
00364               
00365               // get normal direction, this has the edge weight factored into it already
00366               CellTools<double>::getReferenceSideNormal(side_normal , 
00367                                                         (int)side_cur,cell );
00368               
00369               // v.n at cub points on side
00370               vectorBasis->getValues(value_of_v_basis_at_cub_points_side ,
00371                                     cub_points_side_refcell ,
00372                                     OPERATOR_VALUE );
00373               
00374               
00375               for (int i=0;i<numVectorFields;i++) {
00376                 for (int j=0;j<numCubPointsSide;j++) {
00377                   n_of_v_basis_at_cub_points_side(i,j) = 0.0;
00378                   for (int k=0;k<cellDim;k++) {
00379                     n_of_v_basis_at_cub_points_side(i,j) += side_normal(k) * 
00380                       value_of_v_basis_at_cub_points_side(i,j,k);
00381                   }
00382                 } 
00383               }
00384               
00385               cub_weights_side.resize(1,numCubPointsSide);
00386               FunctionSpaceTools::multiplyMeasure<double>(w_n_of_v_basis_at_cub_points_side,
00387                                                           cub_weights_side,
00388                                                           n_of_v_basis_at_cub_points_side);
00389               cub_weights_side.resize(numCubPointsSide);
00390               
00391               FunctionSpaceTools::integrate<double>(rhs_vector_vec,
00392                                                     diri_data_at_cub_points_side,
00393                                                     w_n_of_v_basis_at_cub_points_side,
00394                                                     COMP_BLAS,
00395                                                     false);
00396               for (int i=0;i<numVectorFields;i++) {
00397                 rhs_and_soln_vec(0,i) -= rhs_vector_vec(0,i);
00398               }
00399               
00400             }
00401             
00402             // solve linear system
00403             int info = 0;
00404             Teuchos::LAPACK<int, double> solver;
00405             solver.GESV(numTotalFields, 1, &fe_matrix[0], numTotalFields, &ipiv(0), &rhs_and_soln_vec[0], 
00406                         numTotalFields, &info);
00407             
00408             // compute interpolant; the scalar entries are last
00409             scalarBasis->getValues(value_of_s_basis_at_interp_points,
00410                                   interp_points_ref,
00411                                   OPERATOR_VALUE);
00412             for (int pt=0;pt<numInterpPoints;pt++) {
00413               interpolant(0,pt)=0.0;
00414               for (int i=0;i<numScalarFields;i++) {
00415                 interpolant(0,pt) += rhs_and_soln_vec(0,numVectorFields+i)
00416                   * value_of_s_basis_at_interp_points(i,pt);
00417               }
00418             }
00419             
00420             interp_points_ref.resize(1,numInterpPoints,cellDim);
00421             // get exact solution for comparison
00422             FieldContainer<double> exact_solution(1,numInterpPoints);
00423             u_exact( exact_solution , interp_points_ref , x_order, y_order, z_order);
00424             interp_points_ref.resize(numInterpPoints,cellDim);
00425             
00426             RealSpaceTools<double>::add(interpolant,exact_solution);
00427             
00428             double nrm= RealSpaceTools<double>::vectorNorm(&interpolant[0],interpolant.dimension(1), NORM_TWO);
00429 
00430             *outStream << "\nNorm-2 error between scalar components of exact solution of order ("
00431                        << x_order << ", " << y_order << ", " << z_order
00432                        << ") and finite element interpolant of order " << basis_order << ": "
00433                        << nrm << "\n";
00434 
00435             if (nrm > zero) {
00436               *outStream << "\n\nPatch test failed for solution polynomial order ("
00437                          << x_order << ", " << y_order << ", " << z_order << ") and basis order (scalar, vector)  ("
00438                          << basis_order << ", " << basis_order+1 << ")\n\n";
00439               errorFlag++;
00440             }
00441             
00442           }
00443         }
00444       }
00445     }
00446     
00447   }
00448   
00449   catch (std::logic_error err) {
00450     *outStream << err.what() << "\n\n";
00451     errorFlag = -1000;
00452   };
00453   
00454   if (errorFlag != 0)
00455     std::cout << "End Result: TEST FAILED\n";
00456   else
00457     std::cout << "End Result: TEST PASSED\n";
00458   
00459   // reset format state of std::cout
00460   std::cout.copyfmt(oldFormatState);
00461   
00462   return errorFlag;
00463 }