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