Intrepid
Intrepid_AdaptiveSparseGridDef.hpp
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
49namespace Intrepid {
50
51template<class Scalar, class UserVector>
53 std::vector<int> index,
54 int direction,
55 std::set<std::vector<int> > inOldIndex,
56 AdaptiveSparseGridInterface<Scalar,UserVector> & problem_data) {
57 /*
58 Check if inOldIndex remains admissible if index is added, i.e.
59 index-ek in inOldIndex for all k=1,...,dim.
60 */
61 int dimension = problem_data.getDimension();
62 for (int i=0; i<dimension; i++) {
63 if (index[i]>1 && i!=direction) {
64 index[i]--;
65 if (!inOldIndex.count(index)) {
66 return false;
67 }
68 index[i]++;
69 }
70 }
71 return true;
72}
73
74template<class Scalar, class UserVector>
76 CubatureTensorSorted<Scalar> & outRule,
77 std::vector<int> index,
78 AdaptiveSparseGridInterface<Scalar,UserVector> & problem_data) {
79
80 int numPoints = 0;
81 int dimension = problem_data.getDimension();
82 bool isNormalized = problem_data.isNormalized();
83 std::vector<EIntrepidBurkardt> rule1D; problem_data.getRule(rule1D);
84 std::vector<EIntrepidGrowth> growth1D; problem_data.getGrowth(growth1D);
85
86 for (int i=0; i<dimension; i++) {
87 // Compute 1D rules
88 numPoints = growthRule1D(index[i],growth1D[i],rule1D[i]);
89 CubatureLineSorted<Scalar> diffRule1(rule1D[i],numPoints,isNormalized);
90
91 if (numPoints!=1) { // Compute differential rule
92 numPoints = growthRule1D(index[i]-1,growth1D[i],rule1D[i]);
93 CubatureLineSorted<Scalar> rule1(rule1D[i],numPoints,isNormalized);
94 diffRule1.update(-1.0,rule1,1.0);
95 }
96 // Build Tensor Rule
97 outRule = kron_prod<Scalar>(outRule,diffRule1);
98 }
99}
100
101template<class Scalar, class UserVector>
103 CubatureTensorSorted<Scalar> & outRule,
104 std::vector<int> index, int dimension,
105 std::vector<EIntrepidBurkardt> rule1D,
106 std::vector<EIntrepidGrowth> growth1D,
107 bool isNormalized) {
108
109 int numPoints = 0;
110 for (int i=0; i<dimension; i++) {
111 // Compute 1D rules
112 numPoints = growthRule1D(index[i],growth1D[i],rule1D[i]);
113 CubatureLineSorted<Scalar> diffRule1(rule1D[i],numPoints,isNormalized);
114
115 if (numPoints!=1) { // Differential Rule
116 numPoints = growthRule1D(index[i]-1,growth1D[i],rule1D[i]);
117 CubatureLineSorted<Scalar> rule1(rule1D[i],numPoints,isNormalized);
118 diffRule1.update(-1.0,rule1,1.0);
119 }
120 // Build Tensor Rule
121 outRule = kron_prod<Scalar>(outRule,diffRule1);
122 }
123}
124
125// Update Index Set - no knowledge of active or old indices
126template<class Scalar, class UserVector>
128 typename std::multimap<Scalar,std::vector<int> > & indexSet,
129 UserVector & integralValue,
130 AdaptiveSparseGridInterface<Scalar,UserVector> & problem_data) {
131
132 int dimension = problem_data.getDimension();
133 std::vector<EIntrepidBurkardt> rule1D; problem_data.getRule(rule1D);
134 std::vector<EIntrepidGrowth> growth1D; problem_data.getGrowth(growth1D);
135
136 // Copy Multimap into a Set for ease of use
137 typename std::multimap<Scalar,std::vector<int> >::iterator it;
138 std::set<std::vector<int> > oldSet;
139 std::set<std::vector<int> >::iterator it1(oldSet.begin());
140 for (it=indexSet.begin(); it!=indexSet.end(); it++) {
141 oldSet.insert(it1,it->second);
142 it1++;
143 }
144 indexSet.clear();
145
146 // Find Possible Active Points
147 int flag = 1;
148 std::vector<int> index(dimension,0);
149 typename std::multimap<Scalar,std::vector<int> > activeSet;
150 for (it1=oldSet.begin(); it1!=oldSet.end(); it1++) {
151 index = *it1;
152 for (int i=0; i<dimension; i++) {
153 index[i]++;
154 flag = (int)(!oldSet.count(index));
155 index[i]--;
156 if (flag) {
157 activeSet.insert(std::pair<Scalar,std::vector<int> >(1.0,index));
158 oldSet.erase(it1);
159 break;
160 }
161 }
162 }
163
164 // Compute local and global error indicators for active set
165 typename std::multimap<Scalar,std::vector<int> >::iterator it2;
166 Scalar eta = 0.0;
167 Scalar G = 0.0;
168 Teuchos::RCP<UserVector> s = integralValue.Create();
169 for (it2=activeSet.begin(); it2!=activeSet.end(); it2++) {
170 // Build Differential Quarature Rule
171 index = it2->second;
172 CubatureTensorSorted<Scalar> diffRule(0,dimension);
173 build_diffRule(diffRule,index,problem_data);
174
175 // Apply Rule to function
176 problem_data.eval_cubature(*s,diffRule);
177
178 // Update local error indicator and index set
179 G = problem_data.error_indicator(*s);
180 activeSet.erase(it2);
181 activeSet.insert(it2,std::pair<Scalar,std::vector<int> >(G,index));
182 eta += G;
183 }
184
185 // Refine Sparse Grid
186 eta = refine_grid(activeSet,oldSet,integralValue,eta,
187 dimension,rule1D,growth1D);
188
189 // Insert New Active and Old Index sets into indexSet
190 indexSet.insert(activeSet.begin(),activeSet.end());
191 for (it1=oldSet.begin(); it1!=oldSet.end(); it1++) {
192 index = *it1;
193 indexSet.insert(std::pair<Scalar,std::vector<int> >(-1.0,index));
194 }
195
196 return eta;
197}
198
199// Update index set and output integral
200template<class Scalar, class UserVector>
202 typename std::multimap<Scalar,std::vector<int> > & activeIndex,
203 std::set<std::vector<int> > & oldIndex,
204 UserVector & integralValue,
205 Scalar globalErrorIndicator,
206 AdaptiveSparseGridInterface<Scalar,UserVector> & problem_data) {
207
208 TEUCHOS_TEST_FOR_EXCEPTION((activeIndex.empty()),std::out_of_range,
209 ">>> ERROR (AdaptiveSparseGrid): Active Index set is empty.");
210
211 int dimension = problem_data.getDimension();
212 std::vector<EIntrepidBurkardt> rule1D; problem_data.getRule(rule1D);
213 std::vector<EIntrepidGrowth> growth1D; problem_data.getGrowth(growth1D);
214
215 // Initialize Flags
216 bool maxLevelFlag = true;
217 bool isAdmissibleFlag = true;
218
219 // Initialize Cubature Rule
220 Teuchos::RCP<UserVector> s = integralValue.Create();
221 // Initialize iterator at end of inOldIndex
222 std::set<std::vector<int> >::iterator it1(oldIndex.end());
223
224 // Initialize iterator at end of inActiveIndex
225 typename std::multimap<Scalar,std::vector<int> >::iterator it;
226
227 // Obtain Global Error Indicator as sum of key values of inActiveIndex
228 Scalar eta = globalErrorIndicator;
229
230 // Select Index to refine
231 it = activeIndex.end(); it--; // Decrement to position of final value
232 Scalar G = it->first; // Largest Error Indicator is at End
233 eta -= G; // Update global error indicator
234 std::vector<int> index = it->second; // Get Corresponding index
235 activeIndex.erase(it); // Erase Index from active index set
236 // Insert Index into old index set
237 oldIndex.insert(it1,index); it1 = oldIndex.end();
238
239 // Refinement process
240 for (int k=0; k<dimension; k++) {
241 index[k]++; // index + ek
242 // Check Max Level
243 maxLevelFlag = problem_data.max_level(index);
244 if (maxLevelFlag) {
245 // Check Admissibility
246 isAdmissibleFlag = isAdmissible(index,k,oldIndex,problem_data);
247 if (isAdmissibleFlag) { // If admissible
248 // Build Differential Quarature Rule
249 CubatureTensorSorted<Scalar> diffRule(0,dimension);
250 build_diffRule(diffRule,index,problem_data);
251
252 // Apply Rule to function
253 problem_data.eval_cubature(*s,diffRule);
254
255 // Update integral value
256 integralValue.Update(*s);
257
258 // Update local error indicator and index set
259 G = problem_data.error_indicator(*s);
260 if (activeIndex.end()!=activeIndex.begin())
261 activeIndex.insert(activeIndex.end()--,
262 std::pair<Scalar,std::vector<int> >(G,index));
263 else
264 activeIndex.insert(std::pair<Scalar,std::vector<int> >(G,index));
265
266 // Update global error indicators
267 eta += G;
268 }
269 }
270 else { // Max Level Exceeded
271 //std::cout << "Maximum Level Exceeded" << std::endl;
272 }
273 index[k]--;
274 }
275 return eta;
276}
277
278// Update index set and output integral/sparse grid
279template<class Scalar, class UserVector>
281 typename std::multimap<Scalar,std::vector<int> > & activeIndex,
282 std::set<std::vector<int> > & oldIndex,
283 UserVector & integralValue,
284 CubatureTensorSorted<Scalar> & cubRule,
285 Scalar globalErrorIndicator,
286 AdaptiveSparseGridInterface<Scalar,UserVector> & problem_data) {
287
288 TEUCHOS_TEST_FOR_EXCEPTION((activeIndex.empty()),std::out_of_range,
289 ">>> ERROR (AdaptiveSparseGrid): Active Index set is empty.");
290
291 int dimension = problem_data.getDimension();
292 std::vector<EIntrepidBurkardt> rule1D; problem_data.getRule(rule1D);
293 std::vector<EIntrepidGrowth> growth1D; problem_data.getGrowth(growth1D);
294
295 // Initialize Flags
296 bool maxLevelFlag = true;
297 bool isAdmissibleFlag = true;
298
299 // Initialize Cubature Rule
300 Teuchos::RCP<UserVector> s = integralValue.Create();
301
302 // Initialize iterator at end of inOldIndex
303 std::set<std::vector<int> >::iterator it1(oldIndex.end());
304
305 // Initialize iterator at end of inActiveIndex
306 typename std::multimap<Scalar,std::vector<int> >::iterator it;
307
308 // Obtain Global Error Indicator as sum of key values of inActiveIndex
309 Scalar eta = globalErrorIndicator;
310
311 // Select Index to refine
312 it = activeIndex.end(); it--; // Decrement to position of final value
313 Scalar G = it->first; // Largest Error Indicator is at End
314 eta -= G; // Update global error indicator
315 std::vector<int> index = it->second; // Get Corresponding index
316 activeIndex.erase(it); // Erase Index from active index set
317 // Insert Index into old index set
318 oldIndex.insert(it1,index); it1 = oldIndex.end();
319
320 // Refinement process
321 for (int k=0; k<dimension; k++) {
322 index[k]++; // index + ek
323 // Check Max Level
324 maxLevelFlag = problem_data.max_level(index);
325 if (maxLevelFlag) {
326 // Check Admissibility
327 isAdmissibleFlag = isAdmissible(index,k,oldIndex,problem_data);
328 if (isAdmissibleFlag) { // If admissible
329 // Build Differential Quarature Rule
330 CubatureTensorSorted<Scalar> diffRule(0,dimension);
331 build_diffRule(diffRule,index,problem_data);
332
333 // Apply Rule to function
334 problem_data.eval_cubature(*s,diffRule);
335
336 // Update integral value
337 integralValue.Update(*s);
338
339 // Update local error indicator and index set
340 G = problem_data.error_indicator(*s);
341 if (activeIndex.end()!=activeIndex.begin())
342 activeIndex.insert(activeIndex.end()--,
343 std::pair<Scalar,std::vector<int> >(G,index));
344 else
345 activeIndex.insert(std::pair<Scalar,std::vector<int> >(G,index));
346
347 // Update global error indicators
348 eta += G;
349
350 // Update adapted quadrature rule nodes and weights
351 cubRule.update(1.0,diffRule,1.0);
352 }
353 }
354 else { // Max Level Exceeded
355 //std::cout << "Maximum Level Exceeded" << std::endl;
356 }
357 index[k]--;
358 }
359 return eta;
360}
361
362template<class Scalar, class UserVector>
364 CubatureTensorSorted<Scalar> & output,
365 int dimension, int maxlevel,
366 std::vector<EIntrepidBurkardt> rule1D,
367 std::vector<EIntrepidGrowth> growth1D,
368 bool isNormalized) {
369
370 if (dimension == 2) {
371 std::vector<int> index(dimension,0);
372 for (int i=0; i<maxlevel; i++) {
373 for (int j=0; j<maxlevel; j++) {
374 if(i+j+dimension <= maxlevel+dimension-1) {
375 index[0] = i+1; index[1] = j+1;
376 CubatureTensorSorted<Scalar> diffRule(0,dimension);
377 build_diffRule(diffRule,index,dimension,rule1D,growth1D,isNormalized);
378 output.update(1.0,diffRule,1.0);
379 }
380 }
381 }
382 }
383 else if (dimension == 3) {
384 std::vector<int> index(dimension,0);
385 for (int i=0; i<maxlevel; i++) {
386 for (int j=0; j<maxlevel; j++) {
387 for (int k=0; k<maxlevel; k++) {
388 if(i+j+k+dimension <= maxlevel+dimension-1) {
389 index[0] = i+1; index[1] = j+1; index[2] = k+1;
390 CubatureTensorSorted<Scalar> diffRule(0,dimension);
391 build_diffRule(diffRule,index,dimension,rule1D,
392 growth1D,isNormalized);
393 output.update(1.0,diffRule,1.0);
394 }
395 }
396 }
397 }
398 }
399 else if (dimension == 4) {
400 std::vector<int> index(dimension,0);
401 for (int i=0; i<maxlevel; i++) {
402 for (int j=0; j<maxlevel; j++) {
403 for (int k=0; k<maxlevel; k++) {
404 for (int l=0; l<maxlevel; l++) {
405 if(i+j+k+l+dimension <= maxlevel+dimension-1) {
406 index[0] = i+1; index[1] = j+1; index[2] = k+1; index[3] = l+1;
407 CubatureTensorSorted<Scalar> diffRule(0,dimension);
408 build_diffRule(diffRule,index,dimension,rule1D,
409 growth1D,isNormalized);
410 output.update(1.0,diffRule,1.0);
411 }
412 }
413 }
414 }
415 }
416 }
417 else if (dimension == 5) {
418 std::vector<int> index(dimension,0);
419 for (int i=0; i<maxlevel; i++) {
420 for (int j=0; j<maxlevel; j++) {
421 for (int k=0; k<maxlevel; k++) {
422 for (int l=0; l<maxlevel; l++) {
423 for (int m=0; m<maxlevel; m++) {
424 if(i+j+k+l+m+dimension <= maxlevel+dimension-1) {
425 index[0] = i+1; index[1] = j+1; index[2] = k+1;
426 index[3] = l+1; index[4] = m+1;
427 CubatureTensorSorted<Scalar> diffRule(0,dimension);
428 build_diffRule(diffRule,index,dimension,rule1D,
429 growth1D,isNormalized);
430 output.update(1.0,diffRule,1.0);
431 }
432 }
433 }
434 }
435 }
436 }
437 }
438 else
439 std::cout << "Dimension Must Be Less Than 5\n";
440}
441
442} // end Intrepid namespace
static void buildSparseGrid(CubatureTensorSorted< Scalar > &output, int dimension, int maxlevel, std::vector< EIntrepidBurkardt > rule1D, std::vector< EIntrepidGrowth > growth1D, bool isNormalized)
Build a classic isotropic sparse grid.
static void build_diffRule(CubatureTensorSorted< Scalar > &outRule, std::vector< int > index, AdaptiveSparseGridInterface< Scalar, UserVector > &problem_data)
Given an index, build the corresponding differential cubature rule.
static bool isAdmissible(std::vector< int > index, int direction, std::set< std::vector< int > > inOldIndex, AdaptiveSparseGridInterface< Scalar, UserVector > &problem_data)
Check admissibility of an index set, outputs true if admissible.
static Scalar refine_grid(typename std::multimap< Scalar, std::vector< int > > &indexSet, UserVector &integralValue, AdaptiveSparseGridInterface< Scalar, UserVector > &problem_data)
Update adaptive sparse grid.