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
59using namespace Intrepid;
60
61
62template<class Scalar>
63class StdVector {
64private:
65 Teuchos::RefCountPtr<std::vector<Scalar> > std_vec_;
66
67public:
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
112template<class Scalar, class UserVector>
113class ASGdata :
114 public Intrepid::AdaptiveSparseGridInterface<Scalar,UserVector> {
115public:
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
139long 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
167long 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
182int 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 (const 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}
Header file for the Intrepid::AdaptiveSparseGrid class.
Intrepid utilities.
void eval_integrand(UserVector &output, std::vector< Scalar > &input)
Evaluate the integrand function.
Definition: test_23.cpp:123
Scalar error_indicator(UserVector &input)
User defined error indicator function.
Definition: test_23.cpp:126
bool isNormalized()
Return whether or not cubature weights are normalized.
Builds general adaptive sparse grid rules (Gerstner and Griebel) using the 1D cubature rules in the I...
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 normalize()
Normalize CubatureLineSorted weights.
int getNumPoints() const
Returns the number of cubature points.
int getDimension() const
Returns dimension of domain of integration.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....