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 
49 namespace Intrepid {
50 
51 template<class Scalar, class UserVector>
53  std::vector<int> index,
54  int direction,
55  std::set<std::vector<int> > inOldIndex,
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 
74 template<class Scalar, class UserVector>
77  std::vector<int> index,
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 
101 template<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
126 template<class Scalar, class UserVector>
128  typename std::multimap<Scalar,std::vector<int> > & indexSet,
129  UserVector & integralValue,
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
200 template<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,
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
279 template<class Scalar, class UserVector>
281  typename std::multimap<Scalar,std::vector<int> > & activeIndex,
282  std::set<std::vector<int> > & oldIndex,
283  UserVector & integralValue,
285  Scalar globalErrorIndicator,
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 
362 template<class Scalar, class UserVector>
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.
void getRule(std::vector< EIntrepidBurkardt > &rule1D)
Return user defined 1D quadrature rules.
static Scalar refine_grid(typename std::multimap< Scalar, std::vector< int > > &indexSet, UserVector &integralValue, AdaptiveSparseGridInterface< Scalar, UserVector > &problem_data)
Update adaptive sparse grid.
Utilizes cubature (integration) rules contained in the library sandia_rules (John Burkardt...
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.
void update(Scalar alpha2, CubatureLineSorted< Scalar > &cubRule2, Scalar alpha1)
Replace CubatureLineSorted values with &quot;this = alpha1*this+alpha2*cubRule2&quot;.
virtual Scalar error_indicator(UserVector &input)=0
User defined error indicator function.
int getDimension()
Return dimension of integration domain.
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.
void update(Scalar alpha2, CubatureTensorSorted< Scalar > &cubRule2, Scalar alpha1)
Replace CubatureLineSorted values with &quot;this = alpha1*this+alpha2*cubRule2&quot;.
bool isNormalized()
Return whether or not cubature weights are normalized.
virtual bool max_level(std::vector< int > index)
User defined test for maximum level of cubature.
virtual void eval_cubature(UserVector &output, CubatureTensorSorted< Scalar > &cubRule)
Evaluate the cubature rule.
void getGrowth(std::vector< EIntrepidGrowth > &growth1D)
Return user defined 1D growth rules.
Utilizes 1D cubature (integration) rules contained in the library sandia_rules (John Burkardt...