Intrepid
http://trilinos.sandia.gov/packages/docs/r10.10/packages/intrepid/test/Discretization/Integration/test_23.cpp
Go to the documentation of this file.
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 
00044 
00051 #include "Intrepid_AdaptiveSparseGrid.hpp"
00052 //#include "Intrepid_CubatureLineSorted.hpp"
00053 #include "Intrepid_Utils.hpp"
00054 #include "Teuchos_oblackholestream.hpp"
00055 #include "Teuchos_RCP.hpp"
00056 #include "Teuchos_RefCountPtr.hpp"
00057 #include "Teuchos_GlobalMPISession.hpp"
00058 
00059 using namespace Intrepid;
00060 
00061 
00062 template<class Scalar>
00063 class StdVector {
00064 private:
00065   Teuchos::RefCountPtr<std::vector<Scalar> >  std_vec_;
00066 
00067 public:
00068 
00069   StdVector( const Teuchos::RefCountPtr<std::vector<Scalar> > & std_vec ) 
00070   : std_vec_(std_vec) {}
00071 
00072   Teuchos::RefCountPtr<StdVector<Scalar> > Create() const {
00073     return Teuchos::rcp( new StdVector<Scalar>(
00074            Teuchos::rcp(new std::vector<Scalar>(std_vec_->size(),0))));
00075   }
00076 
00077   void Update( StdVector<Scalar> & s ) {
00078     int dimension  = (int)(std_vec_->size());
00079     for (int i=0; i<dimension; i++)
00080       (*std_vec_)[i] += s[i]; 
00081   }
00082 
00083   void Update( Scalar alpha, StdVector<Scalar> & s )  {
00084     int dimension  = (int)(std_vec_->size());
00085     for (int i=0; i<dimension; i++)
00086       (*std_vec_)[i] += alpha*s[i];
00087   }
00088 
00089   Scalar operator[](int i) {
00090     return (*std_vec_)[i];
00091   }
00092 
00093   void clear() {
00094     std_vec_->clear();
00095   }
00096 
00097   void resize(int n, Scalar p) {
00098     std_vec_->resize(n,p);
00099   }
00100 
00101   int size() {
00102     return (int)std_vec_->size();
00103   }
00104 
00105   void Set( Scalar alpha )  {
00106     int dimension  = (int)(std_vec_->size());
00107     for (int i=0; i<dimension; i++)
00108       (*std_vec_)[i] = alpha;
00109   }
00110 };
00111 
00112 template<class Scalar, class UserVector>
00113 class ASGdata : 
00114   public Intrepid::AdaptiveSparseGridInterface<Scalar,UserVector> {  
00115 public:  
00116   ~ASGdata() {}
00117 
00118   ASGdata(int dimension,std::vector<EIntrepidBurkardt> rule1D,
00119           std::vector<EIntrepidGrowth> growth1D, int maxLevel,
00120           bool isNormalized) : AdaptiveSparseGridInterface<Scalar,UserVector>(
00121           dimension,rule1D,growth1D,maxLevel,isNormalized) {}
00122 
00123   void eval_integrand(UserVector & output, std::vector<Scalar> & input) {
00124     output.clear(); output.resize(1,powl(input[0]+input[1],(long double)6.0));
00125   }  
00126   Scalar error_indicator(UserVector & input) {
00127     int dimension = (int)input.size();
00128     Scalar norm2  = 0.0;
00129     for (int i=0; i<dimension; i++)
00130       norm2 += input[i]*input[i];
00131     
00132     Scalar ID = AdaptiveSparseGridInterface<Scalar,UserVector>::
00133       getInitialDiff();
00134     norm2 = std::sqrt(norm2)/ID;
00135     return norm2;
00136   }
00137 };
00138 
00139 long double adaptSG(StdVector<long double> & iv,
00140                     std::multimap<long double,std::vector<int> > & activeIndex,
00141                     std::set<std::vector<int> > & oldIndex,
00142                     AdaptiveSparseGridInterface<long double,StdVector<long double> > & problem_data,
00143                     CubatureTensorSorted<long double> & cubRule,
00144                     long double TOL) {
00145 
00146   // Construct a Container for the adapted rule
00147   int dimension = problem_data.getDimension();
00148   std::vector<int> index(dimension,1);
00149   
00150   // Initialize global error indicator
00151   long double eta = 1.0;
00152   
00153   // Initialize the Active index set
00154   activeIndex.insert(std::pair<long double,std::vector<int> >(eta,index));
00155 
00156   // Perform Adaptation
00157   while (eta > TOL) {
00158     eta = AdaptiveSparseGrid<long double,StdVector<long double> >::refine_grid(
00159                                                        activeIndex,oldIndex,
00160                                                        iv,cubRule,
00161                                                        eta,problem_data);
00162   }
00163   cubRule.normalize();
00164   return eta;
00165 }
00166 
00167 long double evalQuad(CubatureTensorSorted<long double> & lineCub) {
00168 
00169   int size = lineCub.getNumPoints();
00170   int dimension = lineCub.getDimension();
00171   FieldContainer<long double> cubPoints(size,dimension);
00172   FieldContainer<long double> cubWeights(size);
00173   lineCub.getCubature(cubPoints,cubWeights);
00174   
00175   long double Q = 0.0;
00176   for (int k=0; k<size; k++) 
00177     Q += cubWeights(k)*powl(cubPoints(k,0)+cubPoints(k,1),(long double)6.0);
00178 
00179   return Q;
00180 }
00181 
00182 int main(int argc, char *argv[]) {
00183   
00184   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00185 
00186   // This little trick lets us print to std::cout only if
00187   // a (dummy) command-line argument is provided.
00188   int iprint     = argc - 1;
00189   Teuchos::RCP<std::ostream> outStream;
00190   Teuchos::oblackholestream bhs; // outputs nothing
00191   if (iprint > 0)
00192     outStream = Teuchos::rcp(&std::cout, false);
00193   else
00194     outStream = Teuchos::rcp(&bhs, false);
00195 
00196   // Save the format state of the original std::cout.
00197   Teuchos::oblackholestream oldFormatState;
00198   oldFormatState.copyfmt(std::cout);
00199  
00200   *outStream \
00201   << "===============================================================================\n" \
00202   << "|                                                                             |\n" \
00203   << "|                         Unit Test (AdaptiveSparseGrid)                      |\n" \
00204   << "|                                                                             |\n" \
00205   << "|     1) Integrate a sum of Gaussians in 2D and compare index sets.           |\n" \
00206   << "|                                                                             |\n" \
00207   << "|  Questions? Contact  Drew Kouri (dpkouri@sandia.gov) or                     |\n" \
00208   << "|                      Denis Ridzal (dridzal@sandia.gov).                     |\n" \
00209   << "|                                                                             |\n" \
00210   << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00211   << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00212   << "|                                                                             |\n" \
00213   << "===============================================================================\n"\
00214   << "| TEST 23: Compare index sets for different instances of refine grid          |\n"\
00215   << "===============================================================================\n";
00216 
00217 
00218   // internal variables:
00219   int         errorFlag    = 0;
00220   long double TOL          = INTREPID_TOL;
00221   int         dimension    = 2;
00222   int         maxLevel     = 4;
00223   bool        isNormalized = false;                  
00224 
00225   std::vector<EIntrepidBurkardt> rule1D(dimension,BURK_CLENSHAWCURTIS);
00226   std::vector<EIntrepidGrowth>   growth1D(dimension,GROWTH_FULLEXP);
00227 
00228   ASGdata<long double,StdVector<long double> > problem_data(dimension,rule1D,
00229                                                             growth1D,maxLevel,
00230                                                             isNormalized);  
00231   Teuchos::RCP<std::vector<long double> > integralValue = 
00232     Teuchos::rcp(new std::vector<long double>(1,0.0));
00233   StdVector<long double> sol(integralValue); sol.Set(0.0);
00234   problem_data.init(sol);
00235 
00236   try { 
00237     
00238     // Initialize the index sets
00239     std::multimap<long double,std::vector<int> > activeIndex1;
00240     std::set<std::vector<int> > oldIndex1;  
00241     std::vector<int> index(dimension,1);
00242     CubatureTensorSorted<long double> adaptedRule(dimension,index,rule1D,
00243                                                   growth1D,isNormalized);
00244     adaptSG(sol,activeIndex1,oldIndex1,problem_data,adaptedRule,TOL);
00245     long double Q1  = sol[0];
00246     
00247     CubatureTensorSorted<long double> fullRule(0,dimension);
00248     AdaptiveSparseGrid<long double,StdVector<long double> >::buildSparseGrid(
00249                                                      fullRule,dimension,
00250                                                      maxLevel,rule1D,
00251                                                      growth1D,isNormalized);
00252     long double Q2 = evalQuad(fullRule);
00253     fullRule.normalize();
00254     
00255     long double diff = fabs(Q1-Q2);
00256 
00257     *outStream << "Q1 = " << Q1 << "   Q2 = " << Q2 
00258                << "   |Q1-Q2| = " << diff << "\n";
00259 
00260     int size1 = adaptedRule.getNumPoints(); 
00261     FieldContainer<long double> aPoints(size1,dimension);
00262     FieldContainer<long double> aWeights(size1);
00263     adaptedRule.getCubature(aPoints,aWeights);
00264 
00265     *outStream << "\n\nAdapted Rule Nodes and Weights\n";
00266     for (int i=0; i<size1; i++) 
00267       *outStream << aPoints(i,0) << "\t" << aPoints(i,1) 
00268                  << "\t" << aWeights(i) << "\n";
00269 
00270     int size2 = fullRule.getNumPoints();
00271     FieldContainer<long double> fPoints(size2,dimension);
00272     FieldContainer<long double> fWeights(size2);
00273     fullRule.getCubature(fPoints,fWeights);
00274 
00275     *outStream << "\n\nFull Rule Nodes and Weights\n";
00276     for (int i=0; i<size2; i++) 
00277       *outStream << fPoints(i,0) << "\t" << fPoints(i,1) 
00278                  << "\t" << fWeights(i) << "\n";  
00279 
00280     *outStream << "\n\nSize of adapted rule = " << size1 
00281                << "    Size of full rule = " << size2 << "\n";
00282     if (diff > TOL*fabs(Q2)||size1!=size2) {
00283       errorFlag++;
00284       *outStream << std::right << std::setw(104) << "^^^^---FAILURE!\n";
00285     }
00286     else {
00287       long double sum1 = 0.0, sum2 = 0.0;
00288       for (int i=0; i<size1; i++) {
00289         //diff = fabs(fWeights(i)-aWeights(i));
00290         sum1 += fWeights(i);
00291         sum2 += aWeights(i);
00292       }
00293       *outStream << "Check if weights are normalized:" 
00294                  << "  Adapted Rule Sum = " << sum2
00295                  << "  Full Rule Sum = " << sum1 << "\n";
00296       if (fabs(sum1-1.0) > TOL || fabs(sum2-1.0) > TOL) {
00297         errorFlag++;
00298         *outStream << std::right << std::setw(104) << "^^^^---FAILURE!\n";
00299       }
00300     }
00301   }
00302   catch (std::logic_error err) {    
00303     *outStream << err.what() << "\n";
00304     errorFlag = -1;
00305   };  
00306 
00307   if (errorFlag != 0)
00308     std::cout << "End Result: TEST FAILED\n";
00309   else
00310     std::cout << "End Result: TEST PASSED\n";
00311 
00312   // reset format state of std::cout
00313   std::cout.copyfmt(oldFormatState);
00314 
00315   return errorFlag;
00316 }