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