Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/test/Discretization/Integration/test_21.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 void 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              long double TOL, int flag) {
00144 
00145   // Construct a Container for the adapted rule
00146   int dimension = problem_data.getDimension();
00147   std::vector<int> index(dimension,1);
00148   
00149   // Initialize global error indicator
00150   long double eta = 1.0;
00151   
00152   // Initialize the Active index set
00153   activeIndex.insert(std::pair<long double,std::vector<int> >(eta,index));
00154 
00155   std::vector<EIntrepidBurkardt> rule1D; problem_data.getRule(rule1D);
00156   std::vector<EIntrepidGrowth> growth1D; problem_data.getGrowth(growth1D);
00157   bool isNormalized = problem_data.isNormalized();
00158   CubatureTensorSorted<long double> cubRule(
00159             dimension,index,rule1D,growth1D,isNormalized);
00160   std::cout << "BEFORE\n";
00161   // Perform Adaptation
00162   while (eta > TOL) {
00163     if (flag==1) {
00164       eta = AdaptiveSparseGrid<long double,StdVector<long double> >::refine_grid(activeIndex,oldIndex,iv,eta,problem_data);
00165      }
00166     else if (flag==2) {
00167       eta = AdaptiveSparseGrid<long double,StdVector<long double> >::refine_grid(activeIndex,oldIndex,iv,cubRule,eta,problem_data);
00168      }
00169   }
00170   std::cout << "AFTER\n";
00171 }
00172 
00173 int main(int argc, char *argv[]) {
00174   
00175   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00176 
00177   // This little trick lets us print to std::cout only if
00178   // a (dummy) command-line argument is provided.
00179   int iprint     = argc - 1;
00180   Teuchos::RCP<std::ostream> outStream;
00181   Teuchos::oblackholestream bhs; // outputs nothing
00182   if (iprint > 0)
00183     outStream = Teuchos::rcp(&std::cout, false);
00184   else
00185     outStream = Teuchos::rcp(&bhs, false);
00186 
00187   // Save the format state of the original std::cout.
00188   Teuchos::oblackholestream oldFormatState;
00189   oldFormatState.copyfmt(std::cout);
00190  
00191   *outStream \
00192   << "===============================================================================\n" \
00193   << "|                                                                             |\n" \
00194   << "|                         Unit Test (AdaptiveSparseGrid)                      |\n" \
00195   << "|                                                                             |\n" \
00196   << "|     1) Integrate a sum of Gaussians in 2D and compare index sets.           |\n" \
00197   << "|                                                                             |\n" \
00198   << "|  Questions? Contact  Drew Kouri (dpkouri@sandia.gov) or                     |\n" \
00199   << "|                      Denis Ridzal (dridzal@sandia.gov).                     |\n" \
00200   << "|                                                                             |\n" \
00201   << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00202   << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00203   << "|                                                                             |\n" \
00204   << "===============================================================================\n"\
00205   << "| TEST 21: Compare index sets for different instances of refine grid          |\n"\
00206   << "===============================================================================\n";
00207 
00208 
00209   // internal variables:
00210   int         errorFlag    = 0;
00211   long double TOL          = INTREPID_TOL;
00212   int         dimension    = 2;
00213   int         maxLevel     = 25;
00214   bool        isNormalized = true;                  
00215 
00216   std::vector<EIntrepidBurkardt> rule1D(dimension,BURK_CLENSHAWCURTIS);
00217   std::vector<EIntrepidGrowth>   growth1D(dimension,GROWTH_FULLEXP);
00218 
00219   ASGdata<long double,StdVector<long double> > problem_data(dimension,rule1D,growth1D,maxLevel,isNormalized);  
00220   Teuchos::RCP<std::vector<long double> > integralValue = 
00221     Teuchos::rcp(new std::vector<long double>(1,0.0));
00222   StdVector<long double> sol(integralValue); sol.Set(0.0);
00223   problem_data.init(sol);
00224 
00225   try { 
00226     
00227     // Initialize the index sets
00228     std::multimap<long double,std::vector<int> > activeIndex1;
00229     std::set<std::vector<int> > oldIndex1;  
00230     
00231     adaptSG(sol,activeIndex1,oldIndex1,problem_data,TOL,1);
00232  
00233     // Initialize the index sets
00234     std::multimap<long double,std::vector<int> > activeIndex2;
00235     std::set<std::vector<int> > oldIndex2;  
00236     adaptSG(sol,activeIndex2,oldIndex2,problem_data,TOL,2);
00237     if (activeIndex1.size()!=activeIndex2.size()||oldIndex1.size()!=oldIndex2.size()) {
00238       errorFlag++;
00239       *outStream << std::right << std::setw(104) << "^^^^---FAILURE!\n";
00240     }
00241     else {
00242       std::multimap<long double,std::vector<int> > inter1;
00243       std::set_intersection(activeIndex1.begin(),activeIndex1.end(), 
00244                             activeIndex2.begin(),activeIndex2.end(),
00245                             inserter(inter1,inter1.begin()),inter1.value_comp());
00246       std::set<std::vector<int> > inter2;
00247       std::set_intersection(oldIndex1.begin(),oldIndex1.end(), 
00248                             oldIndex2.begin(),oldIndex2.end(),
00249                             inserter(inter2,inter2.begin()),inter2.value_comp());
00250       
00251       if(inter1.size()!=activeIndex1.size()||inter1.size()!=activeIndex2.size()||
00252          inter2.size()!=oldIndex1.size()||inter2.size()!=oldIndex2.size()) {
00253         errorFlag++;
00254         *outStream << std::right << std::setw(104) << "^^^^---FAILURE!\n";
00255       }
00256       else {
00257         *outStream << "Index sets 1 and " << 2 << " are equal!\n";
00258       } 
00259     }
00260   }
00261   catch (std::logic_error err) {    
00262     *outStream << err.what() << "\n";
00263     errorFlag = -1;
00264   };  
00265 
00266   if (errorFlag != 0)
00267     std::cout << "End Result: TEST FAILED\n";
00268   else
00269     std::cout << "End Result: TEST PASSED\n";
00270 
00271   // reset format state of std::cout
00272   std::cout.copyfmt(oldFormatState);
00273 
00274   return errorFlag;
00275 }