Intrepid
test_23.cpp
Go to the documentation of this file.
1 // @HEADER
2 // ************************************************************************
3 //
4 // Intrepid Package
5 // Copyright (2007) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Pavel Bochev (pbboche@sandia.gov)
38 // Denis Ridzal (dridzal@sandia.gov), or
39 // Kara Peterson (kjpeter@sandia.gov)
40 //
41 // ************************************************************************
42 // @HEADER
43 
44 
52 //#include "Intrepid_CubatureLineSorted.hpp"
53 #include "Intrepid_Utils.hpp"
54 #include "Teuchos_oblackholestream.hpp"
55 #include "Teuchos_RCP.hpp"
56 #include "Teuchos_RefCountPtr.hpp"
57 #include "Teuchos_GlobalMPISession.hpp"
58 
59 using namespace Intrepid;
60 
61 
62 template<class Scalar>
63 class StdVector {
64 private:
65  Teuchos::RefCountPtr<std::vector<Scalar> > std_vec_;
66 
67 public:
68 
69  StdVector( const Teuchos::RefCountPtr<std::vector<Scalar> > & std_vec )
70  : std_vec_(std_vec) {}
71 
72  Teuchos::RefCountPtr<StdVector<Scalar> > Create() const {
73  return Teuchos::rcp( new StdVector<Scalar>(
74  Teuchos::rcp(new std::vector<Scalar>(std_vec_->size(),0))));
75  }
76 
77  void Update( StdVector<Scalar> & s ) {
78  int dimension = (int)(std_vec_->size());
79  for (int i=0; i<dimension; i++)
80  (*std_vec_)[i] += s[i];
81  }
82 
83  void Update( Scalar alpha, StdVector<Scalar> & s ) {
84  int dimension = (int)(std_vec_->size());
85  for (int i=0; i<dimension; i++)
86  (*std_vec_)[i] += alpha*s[i];
87  }
88 
89  Scalar operator[](int i) {
90  return (*std_vec_)[i];
91  }
92 
93  void clear() {
94  std_vec_->clear();
95  }
96 
97  void resize(int n, Scalar p) {
98  std_vec_->resize(n,p);
99  }
100 
101  int size() {
102  return (int)std_vec_->size();
103  }
104 
105  void Set( Scalar alpha ) {
106  int dimension = (int)(std_vec_->size());
107  for (int i=0; i<dimension; i++)
108  (*std_vec_)[i] = alpha;
109  }
110 };
111 
112 template<class Scalar, class UserVector>
113 class ASGdata :
114  public Intrepid::AdaptiveSparseGridInterface<Scalar,UserVector> {
115 public:
116  ~ASGdata() {}
117 
118  ASGdata(int dimension,std::vector<EIntrepidBurkardt> rule1D,
119  std::vector<EIntrepidGrowth> growth1D, int maxLevel,
120  bool isNormalized) : AdaptiveSparseGridInterface<Scalar,UserVector>(
121  dimension,rule1D,growth1D,maxLevel,isNormalized) {}
122 
123  void eval_integrand(UserVector & output, std::vector<Scalar> & input) {
124  output.clear(); output.resize(1,powl(input[0]+input[1],(long double)6.0));
125  }
126  Scalar error_indicator(UserVector & input) {
127  int dimension = (int)input.size();
128  Scalar norm2 = 0.0;
129  for (int i=0; i<dimension; i++)
130  norm2 += input[i]*input[i];
131 
134  norm2 = std::sqrt(norm2)/ID;
135  return norm2;
136  }
137 };
138 
139 long double adaptSG(StdVector<long double> & iv,
140  std::multimap<long double,std::vector<int> > & activeIndex,
141  std::set<std::vector<int> > & oldIndex,
142  AdaptiveSparseGridInterface<long double,StdVector<long double> > & problem_data,
144  long double TOL) {
145 
146  // Construct a Container for the adapted rule
147  int dimension = problem_data.getDimension();
148  std::vector<int> index(dimension,1);
149 
150  // Initialize global error indicator
151  long double eta = 1.0;
152 
153  // Initialize the Active index set
154  activeIndex.insert(std::pair<long double,std::vector<int> >(eta,index));
155 
156  // Perform Adaptation
157  while (eta > TOL) {
159  activeIndex,oldIndex,
160  iv,cubRule,
161  eta,problem_data);
162  }
163  cubRule.normalize();
164  return eta;
165 }
166 
167 long double evalQuad(CubatureTensorSorted<long double> & lineCub) {
168 
169  int size = lineCub.getNumPoints();
170  int dimension = lineCub.getDimension();
171  FieldContainer<long double> cubPoints(size,dimension);
172  FieldContainer<long double> cubWeights(size);
173  lineCub.getCubature(cubPoints,cubWeights);
174 
175  long double Q = 0.0;
176  for (int k=0; k<size; k++)
177  Q += cubWeights(k)*powl(cubPoints(k,0)+cubPoints(k,1),(long double)6.0);
178 
179  return Q;
180 }
181 
182 int main(int argc, char *argv[]) {
183 
184  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
185 
186  // This little trick lets us print to std::cout only if
187  // a (dummy) command-line argument is provided.
188  int iprint = argc - 1;
189  Teuchos::RCP<std::ostream> outStream;
190  Teuchos::oblackholestream bhs; // outputs nothing
191  if (iprint > 0)
192  outStream = Teuchos::rcp(&std::cout, false);
193  else
194  outStream = Teuchos::rcp(&bhs, false);
195 
196  // Save the format state of the original std::cout.
197  Teuchos::oblackholestream oldFormatState;
198  oldFormatState.copyfmt(std::cout);
199 
200  *outStream \
201  << "===============================================================================\n" \
202  << "| |\n" \
203  << "| Unit Test (AdaptiveSparseGrid) |\n" \
204  << "| |\n" \
205  << "| 1) Integrate a sum of Gaussians in 2D and compare index sets. |\n" \
206  << "| |\n" \
207  << "| Questions? Contact Drew Kouri (dpkouri@sandia.gov) or |\n" \
208  << "| Denis Ridzal (dridzal@sandia.gov). |\n" \
209  << "| |\n" \
210  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
211  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
212  << "| |\n" \
213  << "===============================================================================\n"\
214  << "| TEST 23: Compare index sets for different instances of refine grid |\n"\
215  << "===============================================================================\n";
216 
217 
218  // internal variables:
219  int errorFlag = 0;
220  long double TOL = INTREPID_TOL;
221  int dimension = 2;
222  int maxLevel = 4;
223  bool isNormalized = false;
224 
225  std::vector<EIntrepidBurkardt> rule1D(dimension,BURK_CLENSHAWCURTIS);
226  std::vector<EIntrepidGrowth> growth1D(dimension,GROWTH_FULLEXP);
227 
228  ASGdata<long double,StdVector<long double> > problem_data(dimension,rule1D,
229  growth1D,maxLevel,
230  isNormalized);
231  Teuchos::RCP<std::vector<long double> > integralValue =
232  Teuchos::rcp(new std::vector<long double>(1,0.0));
233  StdVector<long double> sol(integralValue); sol.Set(0.0);
234  problem_data.init(sol);
235 
236  try {
237 
238  // Initialize the index sets
239  std::multimap<long double,std::vector<int> > activeIndex1;
240  std::set<std::vector<int> > oldIndex1;
241  std::vector<int> index(dimension,1);
242  CubatureTensorSorted<long double> adaptedRule(dimension,index,rule1D,
243  growth1D,isNormalized);
244  adaptSG(sol,activeIndex1,oldIndex1,problem_data,adaptedRule,TOL);
245  long double Q1 = sol[0];
246 
247  CubatureTensorSorted<long double> fullRule(0,dimension);
249  fullRule,dimension,
250  maxLevel,rule1D,
251  growth1D,isNormalized);
252  long double Q2 = evalQuad(fullRule);
253  fullRule.normalize();
254 
255  long double diff = std::abs(Q1-Q2);
256 
257  *outStream << "Q1 = " << Q1 << " Q2 = " << Q2
258  << " |Q1-Q2| = " << diff << "\n";
259 
260  int size1 = adaptedRule.getNumPoints();
261  FieldContainer<long double> aPoints(size1,dimension);
262  FieldContainer<long double> aWeights(size1);
263  adaptedRule.getCubature(aPoints,aWeights);
264 
265  *outStream << "\n\nAdapted Rule Nodes and Weights\n";
266  for (int i=0; i<size1; i++)
267  *outStream << aPoints(i,0) << "\t" << aPoints(i,1)
268  << "\t" << aWeights(i) << "\n";
269 
270  int size2 = fullRule.getNumPoints();
271  FieldContainer<long double> fPoints(size2,dimension);
272  FieldContainer<long double> fWeights(size2);
273  fullRule.getCubature(fPoints,fWeights);
274 
275  *outStream << "\n\nFull Rule Nodes and Weights\n";
276  for (int i=0; i<size2; i++)
277  *outStream << fPoints(i,0) << "\t" << fPoints(i,1)
278  << "\t" << fWeights(i) << "\n";
279 
280  *outStream << "\n\nSize of adapted rule = " << size1
281  << " Size of full rule = " << size2 << "\n";
282  if (diff > TOL*std::abs(Q2)||size1!=size2) {
283  errorFlag++;
284  *outStream << std::right << std::setw(104) << "^^^^---FAILURE!\n";
285  }
286  else {
287  long double sum1 = 0.0, sum2 = 0.0;
288  for (int i=0; i<size1; i++) {
289  //diff = std::abs(fWeights(i)-aWeights(i));
290  sum1 += fWeights(i);
291  sum2 += aWeights(i);
292  }
293  *outStream << "Check if weights are normalized:"
294  << " Adapted Rule Sum = " << sum2
295  << " Full Rule Sum = " << sum1 << "\n";
296  if (std::abs(sum1-1.0) > TOL || std::abs(sum2-1.0) > TOL) {
297  errorFlag++;
298  *outStream << std::right << std::setw(104) << "^^^^---FAILURE!\n";
299  }
300  }
301  }
302  catch (std::logic_error err) {
303  *outStream << err.what() << "\n";
304  errorFlag = -1;
305  };
306 
307  if (errorFlag != 0)
308  std::cout << "End Result: TEST FAILED\n";
309  else
310  std::cout << "End Result: TEST PASSED\n";
311 
312  // reset format state of std::cout
313  std::cout.copyfmt(oldFormatState);
314 
315  return errorFlag;
316 }
int getNumPoints() const
Returns the number of cubature points.
void normalize()
Normalize CubatureLineSorted weights.
Builds general adaptive sparse grid rules (Gerstner and Griebel) using the 1D cubature rules in the I...
Scalar getInitialDiff()
Return initial error indicator.
Header file for the Intrepid::AdaptiveSparseGrid class.
Intrepid utilities.
Scalar error_indicator(UserVector &input)
User defined error indicator function.
Definition: test_23.cpp:126
Implementation of a templated lexicographical container for a multi-indexed scalar quantity...
int getDimension() const
Returns dimension of domain of integration.
Utilizes 1D cubature (integration) rules contained in the library sandia_rules (John Burkardt...
void getCubature(ArrayPoint &cubPoints, ArrayWeight &cubWeights) const
Returns cubature points and weights (return arrays must be pre-sized/pre-allocated).
void eval_integrand(UserVector &output, std::vector< Scalar > &input)
Evaluate the integrand function.
Definition: test_23.cpp:123