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