Intrepid
test_21.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
61template<class Scalar>
62class StdVector {
63private:
64 Teuchos::RefCountPtr<std::vector<Scalar> > std_vec_;
65
66public:
67
68 StdVector( const Teuchos::RefCountPtr<std::vector<Scalar> > & std_vec )
69 : std_vec_(std_vec) {}
70
71 Teuchos::RefCountPtr<StdVector<Scalar> > Create() const {
72 return Teuchos::rcp( new StdVector<Scalar>(
73 Teuchos::rcp(new std::vector<Scalar>(std_vec_->size(),0))));
74 }
75
76 void Update( StdVector<Scalar> & s ) {
77 int dimension = (int)(std_vec_->size());
78 for (int i=0; i<dimension; i++)
79 (*std_vec_)[i] += s[i];
80 }
81
82 void Update( Scalar alpha, StdVector<Scalar> & s ) {
83 int dimension = (int)(std_vec_->size());
84 for (int i=0; i<dimension; i++)
85 (*std_vec_)[i] += alpha*s[i];
86 }
87
88 Scalar operator[](int i) {
89 return (*std_vec_)[i];
90 }
91
92 void clear() {
93 std_vec_->clear();
94 }
95
96 void resize(int n, Scalar p) {
97 std_vec_->resize(n,p);
98 }
99
100 int size() {
101 return (int)std_vec_->size();
102 }
103
104 void Set( Scalar alpha ) {
105 int dimension = (int)(std_vec_->size());
106 for (int i=0; i<dimension; i++)
107 (*std_vec_)[i] = alpha;
108 }
109};
110
111template<class Scalar, class UserVector>
112class ASGdata :
113 public Intrepid::AdaptiveSparseGridInterface<Scalar,UserVector> {
114public:
115 ~ASGdata() {}
116
117 ASGdata(int dimension,std::vector<EIntrepidBurkardt> rule1D,
118 std::vector<EIntrepidGrowth> growth1D, int maxLevel,
119 bool isNormalized) : AdaptiveSparseGridInterface<Scalar,UserVector>(
120 dimension,rule1D,growth1D,maxLevel,isNormalized) {}
121
122 void eval_integrand(UserVector & output, std::vector<Scalar> & input) {
123 output.clear(); output.resize(1,std::exp(-input[0]*input[0])
124 +10.0*std::exp(-input[1]*input[1]));
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
139void 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,
143 long double TOL, int flag) {
144
145 // Construct a Container for the adapted rule
146 int dimension = problem_data.getDimension();
147 std::vector<int> index(dimension,1);
148
149 // Initialize global error indicator
150 long double eta = 1.0;
151
152 // Initialize the Active index set
153 activeIndex.insert(std::pair<long double,std::vector<int> >(eta,index));
154
155 std::vector<EIntrepidBurkardt> rule1D; problem_data.getRule(rule1D);
156 std::vector<EIntrepidGrowth> growth1D; problem_data.getGrowth(growth1D);
157 bool isNormalized = problem_data.isNormalized();
159 dimension,index,rule1D,growth1D,isNormalized);
160 std::cout << "BEFORE\n";
161 // Perform Adaptation
162 while (eta > TOL) {
163 if (flag==1) {
164 eta = AdaptiveSparseGrid<long double,StdVector<long double> >::refine_grid(activeIndex,oldIndex,iv,eta,problem_data);
165 }
166 else if (flag==2) {
167 eta = AdaptiveSparseGrid<long double,StdVector<long double> >::refine_grid(activeIndex,oldIndex,iv,cubRule,eta,problem_data);
168 }
169 }
170 std::cout << "AFTER\n";
171}
172
173int main(int argc, char *argv[]) {
174
175 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
176
177 // This little trick lets us print to std::cout only if
178 // a (dummy) command-line argument is provided.
179 int iprint = argc - 1;
180 Teuchos::RCP<std::ostream> outStream;
181 Teuchos::oblackholestream bhs; // outputs nothing
182 if (iprint > 0)
183 outStream = Teuchos::rcp(&std::cout, false);
184 else
185 outStream = Teuchos::rcp(&bhs, false);
186
187 // Save the format state of the original std::cout.
188 Teuchos::oblackholestream oldFormatState;
189 oldFormatState.copyfmt(std::cout);
190
191 *outStream \
192 << "===============================================================================\n" \
193 << "| |\n" \
194 << "| Unit Test (AdaptiveSparseGrid) |\n" \
195 << "| |\n" \
196 << "| 1) Integrate a sum of Gaussians in 2D and compare index sets. |\n" \
197 << "| |\n" \
198 << "| Questions? Contact Drew Kouri (dpkouri@sandia.gov) or |\n" \
199 << "| Denis Ridzal (dridzal@sandia.gov). |\n" \
200 << "| |\n" \
201 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
202 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
203 << "| |\n" \
204 << "===============================================================================\n"\
205 << "| TEST 21: Compare index sets for different instances of refine grid |\n"\
206 << "===============================================================================\n";
207
208
209 // internal variables:
210 int errorFlag = 0;
211 long double TOL = INTREPID_TOL;
212 int dimension = 2;
213 int maxLevel = 25;
214 bool isNormalized = true;
215
216 std::vector<EIntrepidBurkardt> rule1D(dimension,BURK_CLENSHAWCURTIS);
217 std::vector<EIntrepidGrowth> growth1D(dimension,GROWTH_FULLEXP);
218
219 ASGdata<long double,StdVector<long double> > problem_data(dimension,rule1D,growth1D,maxLevel,isNormalized);
220 Teuchos::RCP<std::vector<long double> > integralValue =
221 Teuchos::rcp(new std::vector<long double>(1,0.0));
222 StdVector<long double> sol(integralValue); sol.Set(0.0);
223 problem_data.init(sol);
224
225 try {
226
227 // Initialize the index sets
228 std::multimap<long double,std::vector<int> > activeIndex1;
229 std::set<std::vector<int> > oldIndex1;
230
231 adaptSG(sol,activeIndex1,oldIndex1,problem_data,TOL,1);
232
233 // Initialize the index sets
234 std::multimap<long double,std::vector<int> > activeIndex2;
235 std::set<std::vector<int> > oldIndex2;
236 adaptSG(sol,activeIndex2,oldIndex2,problem_data,TOL,2);
237 if (activeIndex1.size()!=activeIndex2.size()||oldIndex1.size()!=oldIndex2.size()) {
238 errorFlag++;
239 *outStream << std::right << std::setw(104) << "^^^^---FAILURE!\n";
240 }
241 else {
242 std::multimap<long double,std::vector<int> > inter1;
243 std::set_intersection(activeIndex1.begin(),activeIndex1.end(),
244 activeIndex2.begin(),activeIndex2.end(),
245 inserter(inter1,inter1.begin()),inter1.value_comp());
246 std::set<std::vector<int> > inter2;
247 std::set_intersection(oldIndex1.begin(),oldIndex1.end(),
248 oldIndex2.begin(),oldIndex2.end(),
249 inserter(inter2,inter2.begin()),inter2.value_comp());
250
251 if(inter1.size()!=activeIndex1.size()||inter1.size()!=activeIndex2.size()||
252 inter2.size()!=oldIndex1.size()||inter2.size()!=oldIndex2.size()) {
253 errorFlag++;
254 *outStream << std::right << std::setw(104) << "^^^^---FAILURE!\n";
255 }
256 else {
257 *outStream << "Index sets 1 and " << 2 << " are equal!\n";
258 }
259 }
260 }
261 catch (const std::logic_error & err) {
262 *outStream << err.what() << "\n";
263 errorFlag = -1;
264 };
265
266 if (errorFlag != 0)
267 std::cout << "End Result: TEST FAILED\n";
268 else
269 std::cout << "End Result: TEST PASSED\n";
270
271 // reset format state of std::cout
272 std::cout.copyfmt(oldFormatState);
273
274 return errorFlag;
275}
Header file for the Intrepid::AdaptiveSparseGrid class.
Intrepid utilities.
void eval_integrand(UserVector &output, std::vector< Scalar > &input)
Evaluate the integrand function.
Definition: test_21.cpp:122
Scalar error_indicator(UserVector &input)
User defined error indicator function.
Definition: test_21.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,...