Intrepid
http://trilinos.sandia.gov/packages/docs/r10.10/packages/intrepid/test/Discretization/Basis/HCURL_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_HCURL_TET_In_FEM.hpp"
00051 #include "Intrepid_DefaultCubatureFactory.hpp"
00052 #include "Intrepid_RealSpaceTools.hpp"
00053 #include "Intrepid_ArrayTools.hpp"
00054 #include "Intrepid_FunctionSpaceTools.hpp"
00055 #include "Intrepid_CellTools.hpp"
00056 #include "Teuchos_oblackholestream.hpp"
00057 #include "Teuchos_RCP.hpp"
00058 #include "Teuchos_GlobalMPISession.hpp"
00059 #include "Teuchos_SerialDenseMatrix.hpp"
00060 #include "Teuchos_SerialDenseVector.hpp"
00061 #include "Teuchos_LAPACK.hpp"
00062 
00063 using namespace std;
00064 using namespace Intrepid;
00065 
00066 void rhsFunc( FieldContainer<double> &, const FieldContainer<double> &, int, int, int, int );
00067 void u_exact( FieldContainer<double> &, const FieldContainer<double> &, int , int, int, int );
00068 
00069 void u_exact( FieldContainer<double> &result, 
00070               const FieldContainer<double> &points,
00071               int comp, 
00072               int xd,
00073               int yd,
00074               int zd)
00075 {
00076   for (int cell=0;cell<result.dimension(0);cell++){
00077     for (int pt=0;pt<result.dimension(1);pt++) {
00078       result(cell,pt,comp) = std::pow(points(cell,pt,0),xd)*std::pow(points(cell,pt,1),yd)
00079         *std::pow(points(cell,pt,2),zd);
00080     }
00081   }
00082   return;
00083 }
00084 
00085 void rhsFunc( FieldContainer<double> & result , 
00086               const FieldContainer<double> &points ,
00087               int comp,
00088               int xd,
00089               int yd,
00090               int zd )
00091 {
00092   u_exact( result , points , comp , xd , yd , zd );
00093 }
00094 
00095 int main(int argc, char *argv[]) {
00096   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00097   
00098   // This little trick lets us print to std::cout only if
00099   // a (dummy) command-line argument is provided.
00100   int iprint = argc - 1;
00101   Teuchos::RCP<std::ostream> outStream;
00102   Teuchos::oblackholestream bhs; // outputs nothing
00103   if (iprint > 0)
00104     outStream = Teuchos::rcp(&std::cout, false);
00105   else
00106     outStream = Teuchos::rcp(&bhs, false);
00107   
00108   // Save the format state of the original std::cout.
00109   Teuchos::oblackholestream oldFormatState;
00110   oldFormatState.copyfmt(std::cout);
00111   
00112   *outStream \
00113     << "===============================================================================\n" \
00114     << "|                                                                             |\n" \
00115     << "|                  Unit Test (Basis_HCURL_TET_In_FEM)                         |\n" \
00116     << "|                                                                             |\n" \
00117     << "| 1) Patch test involving H(curl) matrices                                    |\n" \
00118     << "|                                                                             |\n" \
00119     << "|   Questions? Contact Pavel Bochev (pbboche@sandia.gov),                     |\n" \
00120     << "|                      Robert Kirby (robert.c.kirby@ttu.edu),                 |\n" \
00121     << "|                      Denis Ridzal (dridzal@sandia.gov),                     |\n" \
00122     << "|                      Kara Peterson (kjpeter@sandia.gov).                    |\n" \
00123     << "|                                                                             |\n" \
00124     << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00125     << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00126     << "|                                                                             |\n" \
00127     << "===============================================================================\n" \
00128     << "| TEST 2: Patch test for mass matrices                                        |\n" \
00129     << "===============================================================================\n";
00130   
00131   
00132   int errorFlag = 0;
00133   
00134   outStream -> precision(16);
00135   
00136   try {
00137     DefaultCubatureFactory<double> cubFactory;                                          // create cubature factory
00138     shards::CellTopology cell(shards::getCellTopologyData< shards::Tetrahedron<> >());  // create parent cell topology
00139     
00140     int cellDim = cell.getDimension();
00141     
00142     int min_order = 1;
00143     int max_order = 5;
00144     
00145     int numIntervals = max_order;
00146     int numInterpPoints = ((numIntervals + 1)*(numIntervals + 2)*(numIntervals+3))/6;
00147     FieldContainer<double> interp_points_ref(numInterpPoints, cellDim);
00148     int counter = 0;
00149     for (int j=0; j<=numIntervals; j++) {
00150       for (int i=0; i<=numIntervals-j; i++) {
00151         for (int k=0;k<numIntervals-j-i;k++) {
00152           interp_points_ref(counter,0) = i*(1.0/numIntervals);
00153           interp_points_ref(counter,1) = j*(1.0/numIntervals);
00154           interp_points_ref(counter,2) = k*(1.0/numIntervals);
00155           counter++;
00156         }
00157       }
00158     }
00159     
00160     for (int basis_order=min_order;basis_order<=max_order;basis_order++) {
00161       // create basis
00162       Teuchos::RCP<Basis<double,FieldContainer<double> > > basis =
00163         Teuchos::rcp(new Basis_HCURL_TET_In_FEM<double,FieldContainer<double> >(basis_order,POINTTYPE_EQUISPACED) );
00164       
00165       int numFields = basis->getCardinality();
00166       
00167       // create cubatures
00168       Teuchos::RCP<Cubature<double> > cellCub = cubFactory.create(cell, 2*(basis_order+1));
00169       
00170       int numCubPointsCell = cellCub->getNumPoints();
00171       
00172       // hold cubature information
00173       FieldContainer<double> cub_points_cell(numCubPointsCell, cellDim);
00174       FieldContainer<double> cub_weights_cell(numCubPointsCell);
00175       
00176       // hold basis function information on refcell
00177       FieldContainer<double> value_of_basis_at_cub_points_cell(numFields, numCubPointsCell, cellDim );
00178       FieldContainer<double> w_value_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell, cellDim);
00179 
00180       // holds rhs data
00181       FieldContainer<double> rhs_at_cub_points_cell(1,numCubPointsCell,cellDim);
00182       
00183       // FEM mass matrix
00184       FieldContainer<double> fe_matrix_bak(1,numFields,numFields);
00185       FieldContainer<double> fe_matrix(1,numFields,numFields);
00186       FieldContainer<double> rhs_and_soln_vec(1,numFields);
00187       
00188       FieldContainer<int> ipiv(numFields);
00189       FieldContainer<double> value_of_basis_at_interp_points( numFields , numInterpPoints , cellDim);
00190       FieldContainer<double> interpolant( 1, numInterpPoints , cellDim );
00191 
00192       int info = 0;
00193       Teuchos::LAPACK<int, double> solver;
00194       
00195       // set test tolerance
00196       double zero = (basis_order+1)*(basis_order+1)*1000*INTREPID_TOL;
00197       
00198       // build matrices outside the loop, and then just do the rhs
00199       // for each iteration
00200       cellCub->getCubature(cub_points_cell, cub_weights_cell);
00201       
00202       // need the vector basis
00203       basis->getValues(value_of_basis_at_cub_points_cell,
00204                        cub_points_cell,
00205                        OPERATOR_VALUE);
00206       basis->getValues( value_of_basis_at_interp_points ,
00207                         interp_points_ref ,
00208                         OPERATOR_VALUE );
00209                         
00210                         
00211 
00212 
00213       // construct mass matrix
00214       cub_weights_cell.resize(1,numCubPointsCell);
00215       FunctionSpaceTools::multiplyMeasure<double>(w_value_of_basis_at_cub_points_cell ,
00216                                                   cub_weights_cell ,
00217                                                   value_of_basis_at_cub_points_cell ); 
00218       cub_weights_cell.resize(numCubPointsCell);
00219       
00220       
00221       value_of_basis_at_cub_points_cell.resize( 1 , numFields , numCubPointsCell , cellDim );
00222       FunctionSpaceTools::integrate<double>(fe_matrix_bak,
00223                                             w_value_of_basis_at_cub_points_cell ,
00224                                             value_of_basis_at_cub_points_cell ,
00225                                             COMP_BLAS );
00226       value_of_basis_at_cub_points_cell.resize( numFields , numCubPointsCell , cellDim );
00227       
00228 
00229       //std::cout << fe_matrix_bak << std::endl;
00230 
00231       for (int x_order=0;x_order<basis_order;x_order++) {
00232         for (int y_order=0;y_order<basis_order-x_order;y_order++) {
00233           for (int z_order=0;z_order<basis_order-x_order-y_order;z_order++) {
00234             for (int comp=0;comp<cellDim;comp++) {
00235               fe_matrix.initialize();
00236               // copy mass matrix 
00237               for (int i=0;i<numFields;i++) {
00238                 for (int j=0;j<numFields;j++) {
00239                   fe_matrix(0,i,j) = fe_matrix_bak(0,i,j);
00240                 }
00241               }
00242               
00243               // clear old vector data
00244               rhs_and_soln_vec.initialize();
00245               
00246               // now get rhs vector
00247               
00248               cub_points_cell.resize(1,numCubPointsCell,cellDim);
00249              
00250               rhs_at_cub_points_cell.initialize();
00251               rhsFunc(rhs_at_cub_points_cell,
00252                       cub_points_cell,
00253                       comp, 
00254                       x_order,
00255                       y_order,
00256                       z_order);
00257             
00258               cub_points_cell.resize(numCubPointsCell,cellDim);
00259 
00260               cub_weights_cell.resize(numCubPointsCell);
00261 
00262               FunctionSpaceTools::integrate<double>(rhs_and_soln_vec,
00263                                                     rhs_at_cub_points_cell,
00264                                                     w_value_of_basis_at_cub_points_cell,
00265                                                     COMP_BLAS);
00266               
00267               // solve linear system
00268 
00269 //            solver.GESV(numFields, 1, &fe_matrix[0], numFields, &ipiv(0), &rhs_and_soln_vec[0], 
00270 //                        numFields, &info);
00271               solver.POTRF('L',numFields,&fe_matrix[0],numFields,&info);
00272               solver.POTRS('L',numFields,1,&fe_matrix[0],numFields,&rhs_and_soln_vec[0],numFields,&info);
00273               
00274               interp_points_ref.resize(1,numInterpPoints,cellDim);
00275               // get exact solution for comparison
00276               FieldContainer<double> exact_solution(1,numInterpPoints,cellDim);
00277               exact_solution.initialize();
00278               u_exact( exact_solution , interp_points_ref , comp , x_order, y_order, z_order);
00279               interp_points_ref.resize(numInterpPoints,cellDim);
00280 
00281               // compute interpolant
00282               // first evaluate basis at interpolation points
00283               value_of_basis_at_interp_points.resize(1,numFields,numInterpPoints,cellDim);
00284               FunctionSpaceTools::evaluate<double>( interpolant , 
00285                                                     rhs_and_soln_vec ,
00286                                                     value_of_basis_at_interp_points );
00287               value_of_basis_at_interp_points.resize(numFields,numInterpPoints,cellDim);
00288               
00289               RealSpaceTools<double>::subtract(interpolant,exact_solution);
00290               
00291               double nrm= RealSpaceTools<double>::vectorNorm(&interpolant[0],interpolant.dimension(1), NORM_TWO);
00292               
00293               *outStream << "\nNorm-2 error between scalar components of exact solution of order ("
00294                          << x_order << ", " << y_order << ", " << z_order
00295                          << ") in component " << comp 
00296                          << " and finite element interpolant of order " << basis_order << ": "
00297                          << nrm << "\n";
00298               
00299               if (nrm > zero) {
00300                 *outStream << "\n\nPatch test failed for solution polynomial order ("
00301                            << x_order << ", " << y_order << ", " << z_order << ") and basis order (scalar, vector)  ("
00302                            << basis_order << ", " << basis_order+1 << ")\n\n";
00303                 errorFlag++;
00304               }
00305             }
00306           }
00307         }
00308       }
00309     }
00310     
00311   }
00312   
00313   catch (std::logic_error err) {
00314     *outStream << err.what() << "\n\n";
00315     errorFlag = -1000;
00316   };
00317   
00318   if (errorFlag != 0)
00319     std::cout << "End Result: TEST FAILED\n";
00320   else
00321     std::cout << "End Result: TEST PASSED\n";
00322   
00323   // reset format state of std::cout
00324   std::cout.copyfmt(oldFormatState);
00325   
00326   return errorFlag;
00327 }