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
48
49namespace Intrepid {
50
51template <class Scalar, class ArrayPoint, class ArrayWeight>
52CubatureLineSorted<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
158template <class Scalar, class ArrayPoint, class ArrayWeight>
159CubatureLineSorted<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
263template <class Scalar, class ArrayPoint, class ArrayWeight>
264CubatureLineSorted<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
278template <class Scalar, class ArrayPoint, class ArrayWeight>
279const char* CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::getName() const {
280 return cubature_name_;
281} // end getName
282
283template <class Scalar, class ArrayPoint, class ArrayWeight>
284int CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::getNumPoints() const {
285 return numPoints_;
286} // end getNumPoints
287
288template <class Scalar, class ArrayPoint, class ArrayWeight>
289void CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::getAccuracy(
290 std::vector<int> & accuracy) const {
291 accuracy.assign(1, degree_);
292} // end getAccuracy
293
294template <class Scalar, class ArrayPoint, class ArrayWeight>
295const char* CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::cubature_name_ = "INTREPID_CUBATURE_LINESORTED";
296
297template <class Scalar, class ArrayPoint, class ArrayWeight>
298void 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
309template <class Scalar, class ArrayPoint, class ArrayWeight>
310void 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
318template <class Scalar, class ArrayPoint, class ArrayWeight>
319int CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::getDimension() const {
320 return 1;
321} // end getDimension
322
323template <class Scalar, class ArrayPoint, class ArrayWeight>
324Scalar CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::getNode(
325 typename std::map<Scalar,int>::iterator it) {
326 return it->first;
327} // end getNode
328
329template <class Scalar, class ArrayPoint, class ArrayWeight>
330Scalar CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::getWeight(int node) {
331 return weights_[node];
332} // end getWeight
333
334template <class Scalar, class ArrayPoint, class ArrayWeight>
335Scalar CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::getWeight(
336 Scalar point) {
337 return weights_[points_[point]];
338} // end getWeight
339
340
341template <class Scalar, class ArrayPoint, class ArrayWeight>
342typename std::map<Scalar,int>::iterator CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::begin(void) {
343 return points_.begin();
344} // end begin
345
346template <class Scalar, class ArrayPoint, class ArrayWeight>
347typename std::map<Scalar,int>::iterator CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::end(void) {
348 return points_.end();
349} // end end
350
351template <class Scalar, class ArrayPoint, class ArrayWeight>
352void 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
414int 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
770
771#if defined(Intrepid_SHOW_DEPRECATED_WARNINGS)
772#ifdef __GNUC__
773#warning "The Intrepid package is deprecated"
774#endif
775#endif
776
Intrepid::CubatureLineSorted::begin
std::map< Scalar, int >::iterator begin(void)
Initiate iterator at the beginning of data.
Definition Intrepid_CubatureLineSortedDef.hpp:342
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::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::end
std::map< Scalar, int >::iterator end(void)
Initiate iterator at the end of data.
Definition Intrepid_CubatureLineSortedDef.hpp:347
Intrepid::CubatureLineSorted::getWeight
Scalar getWeight(int weight)
Get a specific weight described by the integer location.
Definition Intrepid_CubatureLineSortedDef.hpp:330
Intrepid::CubatureLineSorted::CubatureLineSorted
CubatureLineSorted(int degree=0, EIntrepidBurkardt rule=BURK_CLENSHAWCURTIS, bool isNormalized=false)
Constructor.
Definition Intrepid_CubatureLineSortedDef.hpp:52
Intrepid::CubatureLineSorted::getName
const char * getName() const
Returns cubature name.
Definition Intrepid_CubatureLineSortedDef.hpp:279
Intrepid::CubatureLineSorted::cubature_name_
static const char * cubature_name_
Cubature name.
Definition Intrepid_CubatureLineSorted.hpp:116
Intrepid::CubatureLineSorted::numPoints_
int numPoints_
Contains the number of nodes for this cubature rule.
Definition Intrepid_CubatureLineSorted.hpp:103
Intrepid::CubatureLineSorted::rule_type_
EIntrepidBurkardt rule_type_
Type of integration points.
Definition Intrepid_CubatureLineSorted.hpp:112
Intrepid::CubatureLineSorted::weights_
std::vector< Scalar > weights_
Contains points of this cubature rule.
Definition Intrepid_CubatureLineSorted.hpp:99
Intrepid::CubatureLineSorted::getDimension
int getDimension() const
Returns dimension of domain of integration.
Definition Intrepid_CubatureLineSortedDef.hpp:319
Intrepid::CubatureLineSorted::update
void update(Scalar alpha2, CubatureLineSorted< Scalar > &cubRule2, Scalar alpha1)
Replace CubatureLineSorted values with "this = alpha1*this+alpha2*cubRule2".
Definition Intrepid_CubatureLineSortedDef.hpp:352
Intrepid::CubatureLineSorted::degree_
int degree_
The degree of polynomials that are integrated exactly by this cubature rule.
Definition Intrepid_CubatureLineSorted.hpp:108
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::points_
std::map< Scalar, int > points_
Contains points of this cubature rule.
Definition Intrepid_CubatureLineSorted.hpp:95
Intrepid::CubatureLineSorted::getNumPoints
int getNumPoints() const
Returns the number of cubature points.
Definition Intrepid_CubatureLineSortedDef.hpp:284
Intrepid::IntrepidBurkardtRules::patterson_lookup
static void patterson_lookup(int n, Scalar x[], Scalar w[])
Gauss-Patterson; returns points and weights.
Definition Intrepid_BurkardtRulesDef.hpp:5114
Intrepid::IntrepidBurkardtRules::hermite_genz_keister_lookup
static void hermite_genz_keister_lookup(int n, Scalar x[], Scalar w[])
Hermite-Genz-Keister; returns points and weights.
Definition Intrepid_BurkardtRulesDef.hpp:1042
Intrepid::IntrepidBurkardtRules::trapezoidal_compute
static void trapezoidal_compute(int n, Scalar x[], Scalar w[])
Trapezoidal rule; returns points and weights.
Definition Intrepid_BurkardtRulesDef.hpp:6322
Intrepid::IntrepidBurkardtRules::fejer2_compute
static void fejer2_compute(int order, Scalar x[], Scalar w[])
Fejer type 2; returns points and weights.
Definition Intrepid_BurkardtRulesDef.hpp:653
Intrepid::IntrepidBurkardtRules::chebyshev1_compute
static void chebyshev1_compute(int order, Scalar x[], Scalar w[])
Gauss-Chebyshev of Type 1; returns points and weights.
Definition Intrepid_BurkardtRulesDef.hpp:62
Intrepid::IntrepidBurkardtRules::chebyshev2_compute
static void chebyshev2_compute(int order, Scalar x[], Scalar w[])
Gauss-Chebyshev of Type 2; returns points and weights.
Definition Intrepid_BurkardtRulesDef.hpp:239
Intrepid::IntrepidBurkardtRules::hermite_compute
static void hermite_compute(int order, Scalar x[], Scalar w[])
Gauss-Hermite; returns points and weights.
Definition Intrepid_BurkardtRulesDef.hpp:877
Intrepid::IntrepidBurkardtRules::laguerre_compute
static void laguerre_compute(int n, Scalar x[], Scalar w[])
Gauss-Laguerre; returns points and weights.
Definition Intrepid_BurkardtRulesDef.hpp:2605
Intrepid::IntrepidBurkardtRules::legendre_compute
static void legendre_compute(int n, Scalar x[], Scalar w[])
Gauss-Legendre; returns points and weights.
Definition Intrepid_BurkardtRulesDef.hpp:3491
Intrepid::IntrepidBurkardtRules::clenshaw_curtis_compute
static void clenshaw_curtis_compute(int order, Scalar x[], Scalar w[])
Clenshaw-Curtis; returns points and weights.
Definition Intrepid_BurkardtRulesDef.hpp:425
Generated by doxygen 1.15.0