Intrepid
Intrepid_CubatureLineSortedDef.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 ArrayPoint, class ArrayWeight>
52 CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::CubatureLineSorted(
53  int degree, EIntrepidBurkardt rule, bool isNormalized ) {
54  TEUCHOS_TEST_FOR_EXCEPTION((degree < 0),std::out_of_range,
55  ">>> ERROR (CubatureLineSorted): No rule implemented for desired polynomial degree.");
56  degree_ = degree;
57  rule_type_ = rule;
58 
59  if (rule==BURK_CHEBYSHEV1||rule==BURK_CHEBYSHEV2||
60  rule==BURK_LAGUERRE ||rule==BURK_LEGENDRE ||
61  rule==BURK_HERMITE) {
62  numPoints_ = (degree+1)/2+1;
63  }
64  else if (rule==BURK_CLENSHAWCURTIS||rule==BURK_FEJER2) {
65  numPoints_ = degree+1;
66  }
67  else if (rule==BURK_TRAPEZOIDAL) {
68  numPoints_ = 2;
69  }
70  else if (rule==BURK_PATTERSON) {
71  int l = 0, o = (degree-0.5)/1.5;
72  for (int i=0; i<8; i++) {
73  l = (int)pow(2.0,(double)i+1.0)-1;
74  if (l>=o) {
75  numPoints_ = l;
76  break;
77  }
78  }
79  }
80  else if (rule==BURK_GENZKEISTER) {
81  int o_ghk[8] = {1,3,9,19,35,37,41,43};
82  int o = (degree-0.5)/1.5;
83  for (int i=0; i<8; i++) {
84  if (o_ghk[i]>=o) {
85  numPoints_ = o_ghk[i];
86  break;
87  }
88  }
89  }
90 
91  Teuchos::Array<Scalar> nodes(numPoints_), weights(numPoints_);
92 
93  if (rule==BURK_CHEBYSHEV1) { // Gauss-Chebyshev Type 1
94  IntrepidBurkardtRules::chebyshev1_compute<Scalar>(numPoints_,
95  nodes.getRawPtr(),weights.getRawPtr());
96  }
97  else if (rule==BURK_CHEBYSHEV2) { // Gauss-Chebyshev Type 2
98  IntrepidBurkardtRules::chebyshev2_compute<Scalar>(numPoints_,
99  nodes.getRawPtr(),weights.getRawPtr());
100  }
101  else if (rule==BURK_CLENSHAWCURTIS) { // Clenshaw-Curtis
102  IntrepidBurkardtRules::clenshaw_curtis_compute<Scalar>(numPoints_,
103  nodes.getRawPtr(),weights.getRawPtr());
104  }
105  else if (rule==BURK_FEJER2) { // Fejer Type 2
106  IntrepidBurkardtRules::fejer2_compute<Scalar>(numPoints_,
107  nodes.getRawPtr(),weights.getRawPtr());
108  }
109  else if (rule==BURK_LEGENDRE) { // Gauss-Legendre
110  IntrepidBurkardtRules::legendre_compute<Scalar>(numPoints_,
111  nodes.getRawPtr(),weights.getRawPtr());
112  }
113  else if (rule==BURK_PATTERSON) { // Gauss-Patterson
114  IntrepidBurkardtRules::patterson_lookup<Scalar>(numPoints_,
115  nodes.getRawPtr(),weights.getRawPtr());
116  }
117  else if (rule==BURK_TRAPEZOIDAL) { // Trapezoidal Rule
118  IntrepidBurkardtRules::trapezoidal_compute<Scalar>(numPoints_,
119  nodes.getRawPtr(),weights.getRawPtr());
120  }
121  else if (rule==BURK_HERMITE) { // Gauss-Hermite
122  IntrepidBurkardtRules::hermite_compute<Scalar>(numPoints_,
123  nodes.getRawPtr(),weights.getRawPtr());
124  }
125  else if (rule==BURK_GENZKEISTER) { // Hermite-Genz-Keister
126  IntrepidBurkardtRules::hermite_genz_keister_lookup<Scalar>(numPoints_,
127  nodes.getRawPtr(),weights.getRawPtr());
128  if (numPoints_>=37) {
129  for (int i=0; i<numPoints_; i++) {
130  weights[i] *= sqrt(M_PI);
131  }
132  }
133  }
134  else if (rule==BURK_LAGUERRE) { // Gauss-Laguerre
135  IntrepidBurkardtRules::laguerre_compute<Scalar>(numPoints_,
136  nodes.getRawPtr(),weights.getRawPtr());
137  }
138 
139  if (isNormalized) {
140  Scalar sum = 0.0;
141  for (int i=0; i<numPoints_; i++) {
142  sum += weights[i];
143  }
144  for (int i=0; i<numPoints_; i++) {
145  weights[i] /= sum;
146  }
147  }
148 
149  points_.clear(); weights_.clear();
150  typename std::map<Scalar,int>::iterator it(points_.begin());
151  for (int i=0; i<numPoints_; i++) {
152  points_.insert(it,std::pair<Scalar,int>(nodes[i],i));
153  weights_.push_back(weights[i]);
154  it = points_.end();
155  }
156 } // end constructor
157 
158 template <class Scalar, class ArrayPoint, class ArrayWeight>
159 CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::CubatureLineSorted(
160  EIntrepidBurkardt rule, int numPoints, bool isNormalized ) {
161  TEUCHOS_TEST_FOR_EXCEPTION((numPoints < 0),std::out_of_range,
162  ">>> ERROR (CubatureLineSorted): No rule implemented for desired number of points.");
163  numPoints_ = numPoints;
164  rule_type_ = rule;
165 
166  Teuchos::Array<Scalar> nodes(numPoints_), weights(numPoints_);
167  if (rule==BURK_CHEBYSHEV1) { // Gauss-Chebyshev Type 1
168  degree_ = 2*numPoints-1;
169  IntrepidBurkardtRules::chebyshev1_compute<Scalar>(numPoints_,
170  nodes.getRawPtr(),weights.getRawPtr());
171  }
172  else if (rule==BURK_CHEBYSHEV2) { // Gauss-Chebyshev Type 2
173  degree_ = 2*numPoints-1;
174  IntrepidBurkardtRules::chebyshev2_compute<Scalar>(numPoints_,
175  nodes.getRawPtr(),weights.getRawPtr());
176  }
177  else if (rule==BURK_CLENSHAWCURTIS) { // Clenshaw-Curtis
178  degree_ = numPoints-1;
179  IntrepidBurkardtRules::clenshaw_curtis_compute<Scalar>(numPoints_,
180  nodes.getRawPtr(),weights.getRawPtr());
181  }
182  else if (rule==BURK_FEJER2) { // Fejer Type 2
183  degree_ = numPoints-1;
184  IntrepidBurkardtRules::fejer2_compute<Scalar>(numPoints_,
185  nodes.getRawPtr(),weights.getRawPtr());
186  }
187  else if (rule==BURK_LEGENDRE) { // Gauss-Legendre
188  degree_ = 2*numPoints-1;
189  IntrepidBurkardtRules::legendre_compute<Scalar>(numPoints_,
190  nodes.getRawPtr(),weights.getRawPtr());
191  }
192  else if (rule==BURK_PATTERSON) { // Gauss-Patterson
193  bool correctNumPoints = false;
194  for (int i=0; i<8; i++) {
195  int l = (int)pow(2.0,(double)i+1.0)-1;
196  if (numPoints==l) {
197  correctNumPoints = true;
198  break;
199  }
200  }
201  TEUCHOS_TEST_FOR_EXCEPTION((correctNumPoints==false),std::out_of_range,
202  ">>> ERROR (CubatureLineSorted): Number of points must be numPoints = 1, 3, 7, 15, 31, 63, 127, 255.");
203  Scalar degree = 1.5*(double)numPoints+0.5;
204  degree_ = (int)degree;
205  IntrepidBurkardtRules::patterson_lookup<Scalar>(numPoints_,
206  nodes.getRawPtr(),weights.getRawPtr());
207  }
208  else if (rule==BURK_TRAPEZOIDAL) { // Trapezoidal Rule
209  degree_ = 2;
210  IntrepidBurkardtRules::trapezoidal_compute<Scalar>(numPoints_,
211  nodes.getRawPtr(),weights.getRawPtr());
212  }
213  else if (rule==BURK_HERMITE) { // Gauss-Hermite
214  degree_ = 2*numPoints-1;
215  IntrepidBurkardtRules::hermite_compute<Scalar>(numPoints_,
216  nodes.getRawPtr(),weights.getRawPtr());
217  }
218  else if (rule==BURK_GENZKEISTER) { // Hermite-Genz-Keister
219  bool correctNumPoints = false;
220  int o_ghk[8] = {1,3,9,19,35,37,41,43};
221  for (int i=0; i<8; i++) {
222  if (o_ghk[i]==numPoints) {
223  correctNumPoints = true;
224  break;
225  }
226  }
227  TEUCHOS_TEST_FOR_EXCEPTION((correctNumPoints==false),std::out_of_range,
228  ">>> ERROR (CubatureLineSorted): Number of points must be numPoints = 1, 3, 9, 35, 37, 41, 43.");
229  Scalar degree = 1.5*(double)numPoints+0.5;
230  degree_ = (int)degree;
231  IntrepidBurkardtRules::hermite_genz_keister_lookup<Scalar>(numPoints_,
232  nodes.getRawPtr(),weights.getRawPtr());
233  if (numPoints_>=37) {
234  for (int i=0; i<numPoints_; i++) {
235  weights[i] *= sqrt(M_PI);
236  }
237  }
238  }
239  else if (rule==BURK_LAGUERRE) { // Gauss-Laguerre
240  degree_ = 2*numPoints-1;
241  IntrepidBurkardtRules::laguerre_compute<Scalar>(numPoints_,
242  nodes.getRawPtr(),weights.getRawPtr());
243  }
244 
245  if (isNormalized) {
246  Scalar sum = 0.0;
247  for (int i=0; i<numPoints_; i++) {
248  sum += weights[i];
249  }
250  for (int i=0; i<numPoints_; i++) {
251  weights[i] /= sum;
252  }
253  }
254  points_.clear(); weights_.clear();
255  typename std::map<Scalar,int>::iterator it(points_.begin());
256  for (int i=0; i<numPoints; i++) {
257  points_.insert(it,std::pair<Scalar,int>(nodes[i],i));
258  weights_.push_back(weights[i]);
259  it = points_.end();
260  }
261 } // end constructor
262 
263 template <class Scalar, class ArrayPoint, class ArrayWeight>
264 CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::CubatureLineSorted(
265  std::vector<Scalar> & points, std::vector<Scalar> & weights) {
266 
267  int size = (int)weights.size();
268  TEUCHOS_TEST_FOR_EXCEPTION(((int)points.size()!=size),std::out_of_range,
269  ">>> ERROR (CubatureLineSorted): Input dimension mismatch.");
270  points_.clear(); weights.clear();
271  for (int loc=0; loc<size; loc++) {
272  points_.insert(std::pair<Scalar,int>(points[loc],loc));
273  weights_.push_back(weights[loc]);
274  }
275  numPoints_ = size;
276 }
277 
278 template <class Scalar, class ArrayPoint, class ArrayWeight>
279 const char* CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::getName() const {
280  return cubature_name_;
281 } // end getName
282 
283 template <class Scalar, class ArrayPoint, class ArrayWeight>
284 int CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::getNumPoints() const {
285  return numPoints_;
286 } // end getNumPoints
287 
288 template <class Scalar, class ArrayPoint, class ArrayWeight>
289 void CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::getAccuracy(
290  std::vector<int> & accuracy) const {
291  accuracy.assign(1, degree_);
292 } // end getAccuracy
293 
294 template <class Scalar, class ArrayPoint, class ArrayWeight>
295 const char* CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::cubature_name_ = "INTREPID_CUBATURE_LINESORTED";
296 
297 template <class Scalar, class ArrayPoint, class ArrayWeight>
298 void CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::getCubature(
299  ArrayPoint & cubPoints, ArrayWeight & cubWeights) const {
300  typename std::map<Scalar,int>::const_iterator it;
301  int i = 0;
302  for (it = points_.begin(); it!=points_.end(); it++) {
303  cubPoints(i) = it->first;
304  cubWeights(i) = weights_[it->second];
305  i++;
306  }
307 } // end getCubature
308 
309 template <class Scalar, class ArrayPoint, class ArrayWeight>
310 void CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::getCubature(ArrayPoint& cubPoints,
311  ArrayWeight& cubWeights,
312  ArrayPoint& cellCoords) const
313 {
314  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
315  ">>> ERROR (CubatureLineSorted): Cubature defined in reference space calling method for physical space cubature.");
316 }
317 
318 template <class Scalar, class ArrayPoint, class ArrayWeight>
319 int CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::getDimension() const {
320  return 1;
321 } // end getDimension
322 
323 template <class Scalar, class ArrayPoint, class ArrayWeight>
324 Scalar CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::getNode(
325  typename std::map<Scalar,int>::iterator it) {
326  return it->first;
327 } // end getNode
328 
329 template <class Scalar, class ArrayPoint, class ArrayWeight>
330 Scalar CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::getWeight(int node) {
331  return weights_[node];
332 } // end getWeight
333 
334 template <class Scalar, class ArrayPoint, class ArrayWeight>
335 Scalar CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::getWeight(
336  Scalar point) {
337  return weights_[points_[point]];
338 } // end getWeight
339 
340 
341 template <class Scalar, class ArrayPoint, class ArrayWeight>
342 typename std::map<Scalar,int>::iterator CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::begin(void) {
343  return points_.begin();
344 } // end begin
345 
346 template <class Scalar, class ArrayPoint, class ArrayWeight>
347 typename std::map<Scalar,int>::iterator CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::end(void) {
348  return points_.end();
349 } // end end
350 
351 template <class Scalar, class ArrayPoint, class ArrayWeight>
352 void CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::update(
353  Scalar alpha2, CubatureLineSorted<Scalar> & cubRule2, Scalar alpha1) {
354 
355  // Initialize an iterator on std::map<Scalar,Scalar>
356  typename std::map<Scalar,int>::iterator it;
357 
358  // Temporary Container for updated rule
359  typename std::map<Scalar,int> newPoints;
360  std::vector<Scalar> newWeights;
361 
362  int loc = 0;
363  Scalar node = 0.0;
364 
365  // Set Intersection rule1 and rule2
366  typename std::map<Scalar,int> inter;
367  std::set_intersection(points_.begin(),points_.end(),
368  cubRule2.begin(),cubRule2.end(),
369  inserter(inter,inter.begin()),inter.value_comp());
370  for (it=inter.begin(); it!=inter.end(); it++) {
371  node = it->first;
372  newWeights.push_back( alpha1*weights_[it->second]
373  +alpha2*cubRule2.getWeight(node));
374  newPoints.insert(std::pair<Scalar,int>(node,loc));
375  loc++;
376  }
377  int isize = inter.size();
378 
379  // Set Difference rule1 \ rule2
380  int size = weights_.size();
381  if (isize!=size) {
382  typename std::map<Scalar,int> diff1;
383  std::set_difference(points_.begin(),points_.end(),
384  cubRule2.begin(),cubRule2.end(),
385  inserter(diff1,diff1.begin()),diff1.value_comp());
386  for (it=diff1.begin(); it!=diff1.end(); it++) {
387  node = it->first;
388  newWeights.push_back(alpha1*weights_[it->second]);
389  newPoints.insert(std::pair<Scalar,int>(node,loc));
390  loc++;
391  }
392  }
393 
394  // Set Difference rule2 \ rule1
395  size = cubRule2.getNumPoints();
396  if(isize!=size) {
397  typename std::map<Scalar,int> diff2;
398  std::set_difference(cubRule2.begin(),cubRule2.end(),
399  points_.begin(),points_.end(),
400  inserter(diff2,diff2.begin()),diff2.value_comp());
401  for (it=diff2.begin(); it!=diff2.end(); it++) {
402  node = it->first;
403  newWeights.push_back(alpha2*cubRule2.getWeight(it->second));
404  newPoints.insert(std::pair<Scalar,int>(node,loc));
405  loc++;
406  }
407  }
408 
409  points_.clear(); points_.insert(newPoints.begin(),newPoints.end());
410  weights_.clear(); weights_.assign(newWeights.begin(),newWeights.end());
411  numPoints_ = (int)points_.size();
412 }
413 
414 int growthRule1D(int index, EIntrepidGrowth growth, EIntrepidBurkardt rule) {
415  //
416  // Compute the growth sequence for 1D quadrature rules according to growth.
417  // For more information on growth rules, see
418  //
419  // J. Burkardt. 1D Quadrature Rules For Sparse Grids.
420  // http://people.sc.fsu.edu/~jburkardt/presentations/sgmga_1d_rules.pdf.
421  //
422  // J. Burkardt. SGMGA: Sparse Grid Mixed Growth Anisotropic Rules.
423  // http://people.sc.fsu.edu/~jburkardt/cpp_src/sgmga/sgmga.html.
424  //
425  // Drew P. Kouri
426  // Sandia National Laboratories - CSRI
427  // May 27, 2011
428  //
429 
430  int level = index-1;
431  //int level = index;
432  if (rule==BURK_CLENSHAWCURTIS) { // Clenshaw-Curtis
433  if (growth==GROWTH_SLOWLIN) {
434  return level+1;
435  }
436  else if (growth==GROWTH_SLOWLINODD) {
437  return 2*((level+1)/2)+1;
438  }
439  else if (growth==GROWTH_MODLIN) {
440  return 2*level+1;
441  }
442  else if (growth==GROWTH_SLOWEXP) {
443  if (level==0) {
444  return 1;
445  }
446  else {
447  int o = 2;
448  while(o<2*level+1) {
449  o = 2*(o-1)+1;
450  }
451  return o;
452  }
453  }
454  else if (growth==GROWTH_MODEXP||growth==GROWTH_DEFAULT) {
455  if (level==0) {
456  return 1;
457  }
458  else {
459  int o = 2;
460  while (o<4*level+1) {
461  o = 2*(o-1)+1;
462  }
463  return o;
464  }
465  }
466  else if (growth==GROWTH_FULLEXP) {
467  if (level==0) {
468  return 1;
469  }
470  else {
471  return (int)pow(2.0,(double)level)+1;
472  }
473  }
474  }
475  else if (rule==BURK_FEJER2) { // Fejer Type 2
476  if (growth==GROWTH_SLOWLIN) {
477  return level+1;
478  }
479  else if (growth==GROWTH_SLOWLINODD) {
480  return 2*((level+1)/2)+1;
481  }
482  else if (growth==GROWTH_MODLIN) {
483  return 2*level+1;
484  }
485  else if (growth==GROWTH_SLOWEXP) {
486  int o = 1;
487  while (o<2*level+1) {
488  o = 2*o+1;
489  }
490  return o;
491  }
492  else if (growth==GROWTH_MODEXP||growth==GROWTH_DEFAULT) {
493  int o = 1;
494  while (o<4*level+1) {
495  o = 2*o+1;
496  }
497  return o;
498  }
499  else if (growth==GROWTH_FULLEXP) {
500  return (int)pow(2.0,(double)level+1.0)-1;
501  }
502  }
503 
504  else if (rule==BURK_PATTERSON) { // Gauss-Patterson
505  if (growth==GROWTH_SLOWLIN||
506  growth==GROWTH_SLOWLINODD||
507  growth==GROWTH_MODLIN) {
508  std::cout << "Specified Growth Rule Not Allowed!\n";
509  return 0;
510  }
511  else if (growth==GROWTH_SLOWEXP) {
512  if (level==0) {
513  return 1;
514  }
515  else {
516  int p = 5;
517  int o = 3;
518  while (p<2*level+1) {
519  p = 2*p+1;
520  o = 2*o+1;
521  }
522  return o;
523  }
524  }
525  else if (growth==GROWTH_MODEXP||growth==GROWTH_DEFAULT) {
526  if (level==0) {
527  return 1;
528  }
529  else {
530  int p = 5;
531  int o = 3;
532  while (p<4*level+1) {
533  p = 2*p+1;
534  o = 2*o+1;
535  }
536  return o;
537  }
538  }
539  else if (growth==GROWTH_FULLEXP) {
540  return (int)pow(2.0,(double)level+1.0)-1;
541  }
542  }
543 
544  else if (rule==BURK_LEGENDRE) { // Gauss-Legendre
545  if (growth==GROWTH_SLOWLIN) {
546  return level+1;
547  }
548  else if (growth==GROWTH_SLOWLINODD) {
549  return 2*((level+1)/2)+1;
550  }
551  else if (growth==GROWTH_MODLIN||growth==GROWTH_DEFAULT) {
552  return 2*level+1;
553  }
554  else if (growth==GROWTH_SLOWEXP) {
555  int o = 1;
556  while (2*o-1<2*level+1) {
557  o = 2*o+1;
558  }
559  return o;
560  }
561  else if (growth==GROWTH_MODEXP) {
562  int o = 1;
563  while (2*o-1<4*level+1) {
564  o = 2*o+1;
565  }
566  return o;
567  }
568  else if (growth==GROWTH_FULLEXP) {
569  return (int)pow(2.0,(double)level+1.0)-1;
570  }
571  }
572 
573  else if (rule==BURK_HERMITE) { // Gauss-Hermite
574  if (growth==GROWTH_SLOWLIN) {
575  return level+1;
576  }
577  else if (growth==GROWTH_SLOWLINODD) {
578  return 2*((level+1)/2)+1;
579  }
580  else if (growth==GROWTH_MODLIN||growth==GROWTH_DEFAULT) {
581  return 2*level+1;
582  }
583  else if (growth==GROWTH_SLOWEXP) {
584  int o = 1;
585  while (2*o-1<2*level+1) {
586  o = 2*o+1;
587  }
588  return o;
589  }
590  else if (growth==GROWTH_MODEXP) {
591  int o = 1;
592  while (2*o-1<4*level+1) {
593  o = 2*o+1;
594  }
595  return o;
596  }
597  else if (growth==GROWTH_FULLEXP) {
598  return (int)pow(2.0,(double)level+1.0)-1;
599  }
600  }
601 
602  else if (rule==BURK_LAGUERRE) { // Gauss-Laguerre
603  if (growth==GROWTH_SLOWLIN) {
604  return level+1;
605  }
606  else if (growth==GROWTH_SLOWLINODD) {
607  return 2*((level+1)/2)+1;
608  }
609  else if (growth==GROWTH_MODLIN||growth==GROWTH_DEFAULT) {
610  return 2*level+1;
611  }
612  else if (growth==GROWTH_SLOWEXP) {
613  int o = 1;
614  while (2*o-1<2*level+1) {
615  o = 2*o+1;
616  }
617  return o;
618  }
619  else if (growth==GROWTH_MODEXP) {
620  int o = 1;
621  while (2*o-1<4*level+1) {
622  o = 2*o+1;
623  }
624  return o;
625  }
626  else if (growth==GROWTH_FULLEXP) {
627  return (int)pow(2.0,(double)level+1.0)-1;
628  }
629  }
630 
631  else if (rule==BURK_CHEBYSHEV1) { // Gauss-Chebyshev Type 1
632  if (growth==GROWTH_SLOWLIN) {
633  return level+1;
634  }
635  else if (growth==GROWTH_SLOWLINODD) {
636  return 2*((level+1)/2)+1;
637  }
638  else if (growth==GROWTH_MODLIN||growth==GROWTH_DEFAULT) {
639  return 2*level+1;
640  }
641  else if (growth==GROWTH_SLOWEXP) {
642  int o = 1;
643  while (2*o-1<2*level+1) {
644  o = 2*o+1;
645  }
646  return o;
647  }
648  else if (growth==GROWTH_MODEXP) {
649  int o = 1;
650  while (2*o-1<4*level+1) {
651  o = 2*o+1;
652  }
653  return o;
654  }
655  else if (growth==GROWTH_FULLEXP) {
656  return (int)pow(2.0,(double)level+1.0)-1;
657  }
658  }
659 
660 
661  else if (rule==BURK_CHEBYSHEV2) { // Gauss-Chebyshev Type 2
662  if (growth==GROWTH_SLOWLIN) {
663  return level+1;
664  }
665  else if (growth==GROWTH_SLOWLINODD) {
666  return 2*((level+1)/2)+1;
667  }
668  else if (growth==GROWTH_MODLIN||growth==GROWTH_DEFAULT) {
669  return 2*level+1;
670  }
671  else if (growth==GROWTH_SLOWEXP) {
672  int o = 1;
673  while (2*o-1<2*level+1) {
674  o = 2*o+1;
675  }
676  return o;
677  }
678  else if (growth==GROWTH_MODEXP) {
679  int o = 1;
680  while (2*o-1<4*level+1) {
681  o = 2*o+1;
682  }
683  return o;
684  }
685  else if (growth==GROWTH_FULLEXP) {
686  return (int)pow(2.0,(double)level+1.0)-1;
687  }
688  }
689 
690  else if (rule==BURK_GENZKEISTER) { // Hermite-Genz-Keister
691  static int o_hgk[5] = { 1, 3, 9, 19, 35 };
692  static int p_hgk[5] = { 1, 5, 15, 29, 51 };
693  if (growth==GROWTH_SLOWLIN||
694  growth==GROWTH_SLOWLINODD||
695  growth==GROWTH_MODLIN) {
696  std::cout << "Specified Growth Rule Not Allowed!\n";
697  return 0;
698  }
699  else if (growth==GROWTH_SLOWEXP) {
700  int l = 0, p = p_hgk[l], o = o_hgk[l];
701  while (p<2*level+1 && l<4) {
702  l++;
703  p = p_hgk[l];
704  o = o_hgk[l];
705  }
706  return o;
707  }
708  else if (growth==GROWTH_MODEXP||growth==GROWTH_DEFAULT) {
709  int l = 0, p = p_hgk[l], o = o_hgk[l];
710  while (p<4*level+1 && l<4) {
711  l++;
712  p = p_hgk[l];
713  o = o_hgk[l];
714  }
715  return o;
716  }
717  else if (growth==GROWTH_FULLEXP) {
718  int l = level; l = std::max(l,0); l = std::min(l,4);
719  return o_hgk[l];
720  }
721  }
722 
723  else if (rule==BURK_TRAPEZOIDAL) { // Trapezoidal
724  if (growth==GROWTH_SLOWLIN) {
725  return level+1;
726  }
727  else if (growth==GROWTH_SLOWLINODD) {
728  return 2*((level+1)/2)+1;
729  }
730  else if (growth==GROWTH_MODLIN) {
731  return 2*level+1;
732  }
733  else if (growth==GROWTH_SLOWEXP) {
734  if (level==0) {
735  return 1;
736  }
737  else {
738  int o = 2;
739  while(o<2*level+1) {
740  o = 2*(o-1)+1;
741  }
742  return o;
743  }
744  }
745  else if (growth==GROWTH_MODEXP||growth==GROWTH_DEFAULT) {
746  if (level==0) {
747  return 1;
748  }
749  else {
750  int o = 2;
751  while (o<4*level+1) {
752  o = 2*(o-1)+1;
753  }
754  return o;
755  }
756  }
757  else if (growth==GROWTH_FULLEXP) {
758  if (level==0) {
759  return 1;
760  }
761  else {
762  return (int)pow(2.0,(double)level)+1;
763  }
764  }
765  }
766  return 0;
767 } // end growthRule1D
768 
769 } // Intrepid namespace
Intrepid::CubatureLineSorted::CubatureLineSorted
CubatureLineSorted(int degree=0, EIntrepidBurkardt rule=BURK_CLENSHAWCURTIS, bool isNormalized=false)
Constructor.
Definition: Intrepid_CubatureLineSortedDef.hpp:52
Intrepid::CubatureLineSorted
Utilizes cubature (integration) rules contained in the library sandia_rules (John Burkardt...
Definition: Intrepid_CubatureLineSorted.hpp:89
Intrepid::CubatureLineSorted::begin
std::map< Scalar, int >::iterator begin(void)
Initiate iterator at the beginning of data.
Definition: Intrepid_CubatureLineSortedDef.hpp:342
Intrepid::CubatureLineSorted::update
void update(Scalar alpha2, CubatureLineSorted< Scalar > &cubRule2, Scalar alpha1)
Replace CubatureLineSorted values with &quot;this = alpha1*this+alpha2*cubRule2&quot;.
Definition: Intrepid_CubatureLineSortedDef.hpp:352
Intrepid::CubatureLineSorted::getAccuracy
void getAccuracy(std::vector< int > &accuracy) const
Returns max. degree of polynomials that are integrated exactly. The return vector has size 1...
Definition: Intrepid_CubatureLineSortedDef.hpp:289
Intrepid::CubatureLineSorted::getCubature
void getCubature(ArrayPoint &cubPoints, ArrayWeight &cubWeights) const
Returns cubature points and weights (return arrays must be pre-sized/pre-allocated).
Definition: Intrepid_CubatureLineSortedDef.hpp:298
Intrepid::CubatureLineSorted::getNumPoints
int getNumPoints() const
Returns the number of cubature points.
Definition: Intrepid_CubatureLineSortedDef.hpp:284
Intrepid::CubatureLineSorted::getDimension
int getDimension() const
Returns dimension of domain of integration.
Definition: Intrepid_CubatureLineSortedDef.hpp:319
Intrepid::CubatureLineSorted::getNode
Scalar getNode(typename std::map< Scalar, int >::iterator it)
Get a specific node described by the iterator location.
Definition: Intrepid_CubatureLineSortedDef.hpp:324
Intrepid::CubatureLineSorted::getName
const char * getName() const
Returns cubature name.
Definition: Intrepid_CubatureLineSortedDef.hpp:279
Intrepid::CubatureLineSorted::getWeight
Scalar getWeight(int weight)
Get a specific weight described by the integer location.
Definition: Intrepid_CubatureLineSortedDef.hpp:330
Intrepid::CubatureLineSorted::end
std::map< Scalar, int >::iterator end(void)
Initiate iterator at the end of data.
Definition: Intrepid_CubatureLineSortedDef.hpp:347
Generated by   doxygen 1.8.5