Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/test/Discretization/TensorProductSpaceTools/test_01.cpp
00001 #include "Teuchos_Array.hpp"
00002 #include "Teuchos_RCP.hpp"
00003 #include "Teuchos_oblackholestream.hpp"
00004 #include "Teuchos_GlobalMPISession.hpp"
00005 #include "Intrepid_TensorProductSpaceTools.hpp"
00006 #include "Intrepid_HGRAD_QUAD_Cn_FEM.hpp"
00007 #include "Intrepid_HGRAD_HEX_Cn_FEM.hpp"
00008 #include "Intrepid_CubaturePolylib.hpp"
00009 #include "Intrepid_Utils.hpp"
00010 #include "Intrepid_Types.hpp"
00011 
00012 
00013 using Teuchos::Array;
00014 using Intrepid::FieldContainer;
00015 using Intrepid::Basis;
00016 using Intrepid::TensorBasis;
00017 
00018 #define INTREPID_TEST_COMMAND( S )                                                                                  \
00019 {                                                                                                                   \
00020   try {                                                                                                             \
00021     S ;                                                                                                             \
00022   }                                                                                                                 \
00023   catch (std::logic_error err) {                                                                                    \
00024       *outStream << "Expected Error ----------------------------------------------------------------\n";            \
00025       *outStream << err.what() << '\n';                                                                             \
00026       *outStream << "-------------------------------------------------------------------------------" << "\n\n";    \
00027   };                                                                                                                \
00028 }
00029 
00030 
00031 int main( int argc , char **argv )
00032 {  
00033   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00034   
00035   // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
00036   int iprint     = argc - 1;
00037   
00038   Teuchos::RCP<std::ostream> outStream;
00039   Teuchos::oblackholestream bhs; // outputs nothing
00040   
00041   if (iprint > 0)
00042     outStream = Teuchos::rcp(&std::cout, false);
00043   else
00044     outStream = Teuchos::rcp(&bhs, false);
00045   
00046   // Save the format state of the original std::cout.
00047   Teuchos::oblackholestream oldFormatState;
00048   oldFormatState.copyfmt(std::cout);
00049 
00050 
00051   *outStream \
00052     << "===============================================================================\n" \
00053     << "|                                                                             |\n" \
00054     << "|                      Unit Test TensorProductSpace Tools                     |\n" \
00055     << "|                                                                             |\n" \
00056     << "|     Tests sum-factored polynomial evaluation and integration                |\n" \
00057     << "|                                                                             |\n" \
00058     << "|  Questions? Contact  Pavel Bochev (pbboche@sandia.gov) or                   |\n" \
00059     << "|                      Denis Ridzal (dridzal@sandia.gov) or                   |\n" \
00060     << "|                      Robert Kirby (robert.c.kirby@ttu.edu)                  |\n" \
00061     << "|                                                                             |\n" \
00062     << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00063     << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00064     << "|                                                                             |\n" \
00065     << "===============================================================================\n";
00066 
00067 
00068 
00069   int errorFlag = 0;
00070 
00071   Array<RCP<TensorBasis<double,FieldContainer<double> > > > basesByDim(4);
00072   Array<RCP<FieldContainer<double> > > cubPtsByDim(4);
00073 
00074   Intrepid::CubaturePolylib<double> cpl(2,Intrepid::PL_GAUSS_LOBATTO);
00075 
00076   FieldContainer<double> cubPoints( cpl.getNumPoints() ,1 );
00077   FieldContainer<double> cubWeights( cpl.getNumPoints() );
00078 
00079   cpl.getCubature( cubPoints, cubWeights );
00080 
00081   basesByDim[2] = Teuchos::rcp( new Intrepid::Basis_HGRAD_QUAD_Cn_FEM<double,FieldContainer<double> >( 2 , Intrepid::POINTTYPE_SPECTRAL ) );
00082   basesByDim[3] = Teuchos::rcp( new Intrepid::Basis_HGRAD_HEX_Cn_FEM<double,FieldContainer<double> >( 2 , Intrepid::POINTTYPE_SPECTRAL ) );
00083 
00084 
00085   // get points
00086   FieldContainer<double> quadPts( cpl.getNumPoints() * cpl.getNumPoints() , 2 );
00087   for (int j=0;j<cpl.getNumPoints();j++)
00088     {
00089       for (int i=0;i<cpl.getNumPoints();i++)
00090         {
00091           int index = j*cpl.getNumPoints() + i;
00092           quadPts(index,0) = cubPoints(i,0);
00093           quadPts(index,1) = cubPoints(j,0);
00094         }
00095     }
00096 
00097   FieldContainer<double> cubPts( cpl.getNumPoints() * cpl.getNumPoints() * cpl.getNumPoints() , 3 );
00098   for (int k=0;k<cpl.getNumPoints();k++)
00099     {
00100       for (int j=0;j<cpl.getNumPoints();j++)
00101         {
00102           for (int i=0;i<cpl.getNumPoints();i++)
00103             {
00104               int index = k* cpl.getNumPoints() * cpl.getNumPoints() + j*cpl.getNumPoints() + i;
00105               cubPts(index,0) = cubPoints(i,0);
00106               cubPts(index,1) = cubPoints(j,0);
00107               cubPts(index,2) = cubPoints(k,0);
00108             }
00109         }
00110     }
00111 
00112   cubPtsByDim[2] = Teuchos::rcp( &quadPts , false );
00113   cubPtsByDim[3] = Teuchos::rcp( &cubPts , false );
00114 
00115   int space_dim = 2;
00116 
00117   Array<Array<RCP<Basis<double,FieldContainer<double> > > > > &bases = basesByDim[space_dim]->getBases();
00118 
00119   FieldContainer<double> coeff(1,1,basesByDim[space_dim]->getCardinality());
00120 
00121 
00122   
00123   Array<RCP<FieldContainer<double> > > pts( space_dim );
00124   pts[0] = Teuchos::rcp( &cubPoints, false );
00125   for (int i=1;i<space_dim;i++)
00126     {
00127       pts[i] = pts[0];
00128     }
00129 
00130   Array<RCP<FieldContainer<double> > > wts(space_dim);
00131   wts[0] = Teuchos::rcp( &cubWeights , false );
00132   for (int i=1;i<space_dim;i++)
00133     {
00134       wts[i] = wts[0];
00135     }
00136 
00137   FieldContainer<double> Phix(bases[0][0]->getCardinality(),
00138                               cpl.getNumPoints() );
00139   FieldContainer<double> Phiy(bases[0][1]->getCardinality(),
00140                               cpl.getNumPoints() );
00141   FieldContainer<double> DPhix(bases[0][0]->getCardinality(),
00142                                cpl.getNumPoints(), 1 );
00143   FieldContainer<double> DPhiy(bases[0][1]->getCardinality(),
00144                                cpl.getNumPoints(), 1 );
00145 
00146   bases[0][0]->getValues( Phix , cubPoints, Intrepid::OPERATOR_VALUE );
00147   bases[0][1]->getValues( Phiy , cubPoints, Intrepid::OPERATOR_VALUE );
00148   bases[0][0]->getValues( DPhix , cubPoints, Intrepid::OPERATOR_D1 );
00149   bases[0][1]->getValues( DPhiy , cubPoints, Intrepid::OPERATOR_D1 );
00150 
00151   Array<RCP<FieldContainer<double> > > basisVals(2);
00152   basisVals[0] = Teuchos::rcp( &Phix , false );
00153   basisVals[1] = Teuchos::rcp( &Phiy , false );
00154 
00155   Array<RCP<FieldContainer<double> > > basisDVals(2);
00156   basisDVals[0] = Teuchos::rcp( &DPhix , false );
00157   basisDVals[1] = Teuchos::rcp( &DPhiy , false );
00158 
00159   FieldContainer<double> vals(1,1,pts[0]->size() * pts[1]->size() );
00160 
00161   // first basis function is the polynomial.
00162   coeff(0,0,0) = 1.0;
00163 
00164   Intrepid::TensorProductSpaceTools::evaluate<double,FieldContainer<double>,FieldContainer<double>,FieldContainer<double> >( vals , coeff , basisVals );
00165 
00166   FieldContainer<double> grads( 1 , 1, pts[0]->size() * pts[1]->size() , 2 );
00167 
00168   Intrepid::TensorProductSpaceTools::evaluateGradient<double,FieldContainer<double>,FieldContainer<double>,FieldContainer<double> >( grads , 
00169                                                                                                                                      coeff , 
00170                                                                                                                                      basisVals ,
00171                                                                                                                                      basisDVals );
00172 
00173   // confirm by comparing to actual gradients
00174   FieldContainer<double> fullVals(basesByDim[space_dim]->getCardinality(),
00175                                   basesByDim[space_dim]->getCardinality());
00176   FieldContainer<double> fullGrads(basesByDim[space_dim]->getCardinality(),
00177                                    basesByDim[space_dim]->getCardinality(),
00178                                    space_dim );
00179 
00180 
00181   basesByDim[space_dim]->getValues( fullVals ,
00182                     quadPts ,
00183                     Intrepid::OPERATOR_VALUE );
00184   basesByDim[space_dim]->getValues( fullGrads ,
00185                     quadPts ,
00186                     Intrepid::OPERATOR_GRAD );
00187 
00188   for (int i=0;i<fullVals.dimension(1);i++)
00189     {
00190       if (std::abs( fullVals(0,i) - vals(0,0,i) ) > Intrepid::INTREPID_TOL ) 
00191         {
00192           errorFlag++;
00193           *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00194           
00195           // Output the multi-index of the value where the error is:
00196           *outStream << " Evaluating first bf At multi-index { ";
00197           *outStream << i;
00198           *outStream << "}  brute force value: " << fullVals(0,i)
00199                      << " but tensor-product  value: " << vals(0,i) << "\n";
00200           *outStream << "Difference: " << std::abs( fullVals(0,i) - vals(0,i) ) << "\n";
00201           }
00202     }
00203 
00204   for (int i=0;i<fullGrads.dimension(1);i++)
00205     {
00206       for (int j=0;j<fullGrads.dimension(2);j++)
00207         {
00208           if (std::abs( fullGrads(0,i,j) - grads(0,0,i,j) ) > Intrepid::INTREPID_TOL ) 
00209             {
00210               errorFlag++;
00211               *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00212               
00213               // Output the multi-index of the value where the error is:
00214               *outStream << " Evaluating first bf At multi-index { ";
00215               *outStream << i << " " << j;
00216               *outStream << "}  brute force value: " << fullGrads(0,i,j)
00217                          << " but tensor-product  value: " << grads(0,0,i,j) << "\n";
00218               *outStream << "Difference: " << std::abs( fullGrads(0,i,j) - grads(0,0,i,j) ) << "\n";
00219             }
00220         }
00221     }
00222 
00223   
00224   // now test moments.  
00225   // I've already evaluated the first basis function at the quadrature points.
00226   // why not use it?
00227 
00228   FieldContainer<double> momentsNaive(1,basesByDim[2]->getCardinality());
00229   for (int i=0;i<basesByDim[2]->getCardinality();i++)
00230     {
00231       momentsNaive(0,i) = 0.0;
00232       for (int qpty=0;qpty<cubPoints.dimension(0);qpty++)
00233         {
00234           for (int qptx=0;qptx<cubPoints.dimension(0);qptx++)
00235             {
00236               momentsNaive(0,i) += cubWeights(qpty) * cubWeights(qptx) *
00237                 vals( 0, 0, qpty*cubPoints.dimension(0)+qptx )
00238                 * fullVals(i,qpty*cubPoints.dimension(0)+qptx);
00239             }
00240         }
00241     }
00242 
00243   FieldContainer<double> momentsClever(1,1,basesByDim[space_dim]->getCardinality());
00244   Intrepid::TensorProductSpaceTools::moments<double,FieldContainer<double>,FieldContainer<double>,FieldContainer<double>,FieldContainer<double> >( momentsClever , 
00245                                                                                                                                                    vals ,
00246                                                                                                                                                    basisVals ,
00247                                                                                                                                                    wts );
00248   for (int j=0;j<momentsClever.dimension(0);j++)
00249     {
00250       if (std::abs( momentsClever(0,0,j) - momentsNaive(0,j) ) > Intrepid::INTREPID_TOL ) 
00251         {
00252           errorFlag++;
00253           *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00254           
00255           // Output the multi-index of the value where the error is:
00256           *outStream << " At multi-index { ";
00257           *outStream << " " << j;
00258           *outStream << "}  brute force value: " << momentsNaive(0,j)
00259                      << " but sum-factored value: " << momentsClever(0,0,j) << "\n";
00260           *outStream << "Difference: " << std::abs( momentsNaive(0,j) - momentsClever(0,0,j) ) << "\n";
00261         }
00262     }
00263 
00264   if (errorFlag != 0)
00265     {
00266       std::cout << "End Result: TEST FAILED\n";
00267     }
00268   else
00269     {
00270       std::cout << "End Result: TEST PASSED\n";
00271     }
00272   
00273   // reset format state of std::cout
00274   std::cout.copyfmt(oldFormatState);
00275   
00276   return errorFlag;
00277 
00278 }