Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/test/Discretization/Integration/test_20.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 template<class Scalar>
00062 class StdVector {
00063 private:
00064   Teuchos::RefCountPtr<std::vector<Scalar> >  std_vec_;
00065 
00066 public:
00067 
00068   StdVector( const Teuchos::RefCountPtr<std::vector<Scalar> > & std_vec ) 
00069   : std_vec_(std_vec) {}
00070 
00071   Teuchos::RefCountPtr<StdVector<Scalar> > Create() const {
00072     return Teuchos::rcp( new StdVector<Scalar>(
00073            Teuchos::rcp(new std::vector<Scalar>(std_vec_->size(),0))));
00074   }
00075 
00076   void Update( StdVector<Scalar> & s ) {
00077     int dimension  = (int)(std_vec_->size());
00078     for (int i=0; i<dimension; i++)
00079       (*std_vec_)[i] += s[i]; 
00080   }
00081 
00082   void Update( Scalar alpha, StdVector<Scalar> & s )  {
00083     int dimension  = (int)(std_vec_->size());
00084     for (int i=0; i<dimension; i++)
00085       (*std_vec_)[i] += alpha*s[i];
00086   }
00087 
00088   Scalar operator[](int i) {
00089     return (*std_vec_)[i];
00090   }
00091 
00092   void clear() {
00093     std_vec_->clear();
00094   }
00095 
00096   void resize(int n, Scalar p) {
00097     std_vec_->resize(n,p);
00098   }
00099 
00100   int size() {
00101     return (int)std_vec_->size();
00102   }
00103 
00104   void Set( Scalar alpha )  {
00105     int dimension  = (int)(std_vec_->size());
00106     for (int i=0; i<dimension; i++)
00107       (*std_vec_)[i] = alpha;
00108   }
00109 };
00110 
00111 template<class Scalar, class UserVector>
00112 class ASGdata : 
00113   public Intrepid::AdaptiveSparseGridInterface<Scalar,UserVector> {  
00114 public:  
00115   ~ASGdata() {}
00116 
00117   ASGdata(int dimension,std::vector<EIntrepidBurkardt> rule1D,
00118           std::vector<EIntrepidGrowth> growth1D, int maxLevel,
00119           bool isNormalized) : AdaptiveSparseGridInterface<Scalar,UserVector>(
00120           dimension,rule1D,growth1D,maxLevel,isNormalized) {}
00121 
00122   void eval_integrand(UserVector & output, std::vector<Scalar> & input) {
00123     output.clear(); output.resize(1,std::exp(-input[0]*input[0])
00124                                   +10.0*std::exp(-input[1]*input[1]));
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           AdaptiveSparseGridInterface<long double,StdVector<long double> > & 
00141           problem_data,long double TOL) {
00142 
00143   // Construct a Container for the adapted rule
00144   int dimension = problem_data.getDimension();
00145   std::vector<int> index(dimension,1);
00146   
00147   // Initialize global error indicator
00148   long double eta = 1.0;
00149   
00150   // Initialize the Active index set
00151   std::multimap<long double,std::vector<int> > activeIndex;  
00152   activeIndex.insert(std::pair<long double,std::vector<int> >(eta,index));
00153 
00154   // Initialize the old index set
00155   std::set<std::vector<int> > oldIndex;
00156 
00157   // Perform Adaptation
00158   while (eta > TOL) {
00159     eta = AdaptiveSparseGrid<long double,StdVector<long double> >::refine_grid(
00160         activeIndex,oldIndex,iv,eta,problem_data);
00161   }
00162   return eta;
00163 }
00164 
00165 int main(int argc, char *argv[]) {
00166   
00167   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00168 
00169   // This little trick lets us print to std::cout only if
00170   // a (dummy) command-line argument is provided.
00171   int iprint     = argc - 1;
00172   Teuchos::RCP<std::ostream> outStream;
00173   Teuchos::oblackholestream bhs; // outputs nothing
00174   if (iprint > 0)
00175     outStream = Teuchos::rcp(&std::cout, false);
00176   else
00177     outStream = Teuchos::rcp(&bhs, false);
00178 
00179   // Save the format state of the original std::cout.
00180   Teuchos::oblackholestream oldFormatState;
00181   oldFormatState.copyfmt(std::cout);
00182  
00183   *outStream \
00184   << "===============================================================================\n" \
00185   << "|                                                                             |\n" \
00186   << "|                         Unit Test (AdaptiveSparseGrid)                      |\n" \
00187   << "|                                                                             |\n" \
00188   << "|     1) Integrate a sum of Gaussians in 2D (Gerstner and Griebel).           |\n" \
00189   << "|                                                                             |\n" \
00190   << "|  Questions? Contact  Drew Kouri (dpkouri@sandia.gov) or                     |\n" \
00191   << "|                      Denis Ridzal (dridzal@sandia.gov).                     |\n" \
00192   << "|                                                                             |\n" \
00193   << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00194   << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00195   << "|                                                                             |\n" \
00196   << "===============================================================================\n"\
00197   << "| TEST 20: Integrate an anisotropic sum of Gaussians in 2D                    |\n"\
00198   << "===============================================================================\n";
00199 
00200 
00201   // internal variables:
00202   int         errorFlag    = 0;
00203   long double TOL          = INTREPID_TOL;
00204   int         dimension    = 2;
00205   int         maxLevel     = 25;
00206   bool        isNormalized = true;                  
00207 
00208   std::vector<EIntrepidBurkardt> rule1D(dimension,BURK_CLENSHAWCURTIS);
00209   std::vector<EIntrepidGrowth>   growth1D(dimension,GROWTH_FULLEXP);
00210  
00211   ASGdata<long double,StdVector<long double> > problem_data(
00212     dimension,rule1D,growth1D,maxLevel,isNormalized);  
00213   Teuchos::RCP<std::vector<long double> > integralValue = 
00214     Teuchos::rcp(new std::vector<long double>(1,0.0));
00215   StdVector<long double> sol(integralValue); sol.Set(0.0);
00216   problem_data.init(sol);
00217 
00218   long double eta = adaptSG(sol,problem_data,TOL); 
00219 
00220   long double analyticInt = (1.0+10.0)*std::sqrt(M_PI)/2.0*erff(1.0);
00221   long double abstol      = 1.0e1*std::sqrt(INTREPID_TOL);
00222   long double absdiff     = fabs(analyticInt-sol[0]);
00223   try { 
00224     *outStream << "Adaptive Sparse Grid exited with global error " 
00225                << std::scientific << std::setprecision(16) << eta << "\n"
00226                << "Approx = " << std::scientific << std::setprecision(16) << sol[0] 
00227                << ",  Exact = " << std::scientific << std::setprecision(16) << analyticInt << "\n"
00228                << "Error = " << std::scientific << std::setprecision(16) << absdiff << "   " 
00229                << "<?" << "   " << abstol << "\n"; 
00230     if (absdiff > abstol) {
00231       errorFlag++;
00232       *outStream << std::right << std::setw(104) << "^^^^---FAILURE!\n";
00233     }
00234   }
00235   catch (std::logic_error err) {    
00236     *outStream << err.what() << "\n";
00237     errorFlag = -1;
00238   };  
00239 
00240   if (errorFlag != 0)
00241     std::cout << "End Result: TEST FAILED\n";
00242   else
00243     std::cout << "End Result: TEST PASSED\n";
00244 
00245   // reset format state of std::cout
00246   std::cout.copyfmt(oldFormatState);
00247 
00248   return errorFlag;
00249 }