Intrepid
http://trilinos.sandia.gov/packages/docs/r10.10/packages/intrepid/src/Discretization/Integration/Intrepid_CubatureLineSortedDef.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ************************************************************************
00003 //
00004 //                           Intrepid Package
00005 //                 Copyright (2007) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Pavel Bochev  (pbboche@sandia.gov)
00038 //                    Denis Ridzal  (dridzal@sandia.gov), or
00039 //                    Kara Peterson (kjpeter@sandia.gov)
00040 //
00041 // ************************************************************************
00042 // @HEADER
00043 
00049 namespace Intrepid {
00050 
00051 template <class Scalar, class ArrayPoint, class ArrayWeight> 
00052 CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::CubatureLineSorted(
00053                       int degree, EIntrepidBurkardt rule, bool isNormalized ) {
00054   TEUCHOS_TEST_FOR_EXCEPTION((degree < 0),std::out_of_range,
00055     ">>> ERROR (CubatureLineSorted): No rule implemented for desired polynomial degree.");
00056   degree_    = degree;
00057   rule_type_ = rule;
00058  
00059   if (rule==BURK_CHEBYSHEV1||rule==BURK_CHEBYSHEV2||
00060       rule==BURK_LAGUERRE  ||rule==BURK_LEGENDRE  ||
00061       rule==BURK_HERMITE) {
00062     numPoints_ = (degree+1)/2+1;
00063   }
00064   else if (rule==BURK_CLENSHAWCURTIS||rule==BURK_FEJER2) {
00065     numPoints_ = degree+1;
00066   }
00067   else if (rule==BURK_TRAPEZOIDAL) {
00068     numPoints_ = 2;
00069   }
00070   else if (rule==BURK_PATTERSON) {
00071     int l = 0, o = (degree-0.5)/1.5;
00072     for (int i=0; i<8; i++) {
00073       l = (int)pow(2.0,(double)i+1.0)-1;
00074       if (l>=o) {
00075         numPoints_ = l;
00076         break;
00077       }
00078     }
00079   }
00080   else if (rule==BURK_GENZKEISTER) {
00081     int o_ghk[8] = {1,3,9,19,35,37,41,43}; 
00082     int o = (degree-0.5)/1.5;
00083     for (int i=0; i<8; i++) {
00084       if (o_ghk[i]>=o) {
00085         numPoints_ = o_ghk[i];
00086         break;
00087       }
00088     }
00089   }
00090 
00091   Teuchos::Array<Scalar> nodes(numPoints_), weights(numPoints_);
00092 
00093   if (rule==BURK_CHEBYSHEV1) { // Gauss-Chebyshev Type 1
00094     IntrepidBurkardtRules::chebyshev1_compute<Scalar>(numPoints_,
00095                                     nodes.getRawPtr(),weights.getRawPtr());
00096   }
00097   else if (rule==BURK_CHEBYSHEV2) { // Gauss-Chebyshev Type 2
00098     IntrepidBurkardtRules::chebyshev2_compute<Scalar>(numPoints_,
00099                                     nodes.getRawPtr(),weights.getRawPtr());
00100   }
00101   else if (rule==BURK_CLENSHAWCURTIS) { // Clenshaw-Curtis    
00102     IntrepidBurkardtRules::clenshaw_curtis_compute<Scalar>(numPoints_,
00103                                     nodes.getRawPtr(),weights.getRawPtr());
00104   }
00105   else if (rule==BURK_FEJER2) { // Fejer Type 2
00106     IntrepidBurkardtRules::fejer2_compute<Scalar>(numPoints_,
00107                                     nodes.getRawPtr(),weights.getRawPtr());
00108   }
00109   else if (rule==BURK_LEGENDRE) { // Gauss-Legendre
00110     IntrepidBurkardtRules::legendre_compute<Scalar>(numPoints_,
00111                                     nodes.getRawPtr(),weights.getRawPtr());
00112   }
00113   else if (rule==BURK_PATTERSON) { // Gauss-Patterson
00114     IntrepidBurkardtRules::patterson_lookup<Scalar>(numPoints_,
00115                                     nodes.getRawPtr(),weights.getRawPtr());
00116   }
00117   else if (rule==BURK_TRAPEZOIDAL) { // Trapezoidal Rule
00118     IntrepidBurkardtRules::trapezoidal_compute<Scalar>(numPoints_,
00119                                     nodes.getRawPtr(),weights.getRawPtr());
00120   }
00121   else if (rule==BURK_HERMITE) { // Gauss-Hermite
00122     IntrepidBurkardtRules::hermite_compute<Scalar>(numPoints_, 
00123                                     nodes.getRawPtr(),weights.getRawPtr());
00124   }
00125   else if (rule==BURK_GENZKEISTER) { // Hermite-Genz-Keister
00126     IntrepidBurkardtRules::hermite_genz_keister_lookup<Scalar>(numPoints_,
00127                                     nodes.getRawPtr(),weights.getRawPtr());
00128     if (numPoints_>=37) {
00129       for (int i=0; i<numPoints_; i++) {
00130         weights[i] *= sqrt(M_PI);
00131       }
00132     }
00133   }
00134   else if (rule==BURK_LAGUERRE) { // Gauss-Laguerre
00135     IntrepidBurkardtRules::laguerre_compute<Scalar>(numPoints_,
00136                                     nodes.getRawPtr(),weights.getRawPtr());
00137   }
00138 
00139   if (isNormalized) {
00140     Scalar sum = 0.0;
00141     for (int i=0; i<numPoints_; i++) {
00142       sum += weights[i];
00143     }
00144     for (int i=0; i<numPoints_; i++) {
00145       weights[i] /= sum;
00146     }
00147   }
00148 
00149   points_.clear(); weights_.clear();
00150   typename std::map<Scalar,int>::iterator it(points_.begin());
00151   for (int i=0; i<numPoints_; i++) {
00152     points_.insert(it,std::pair<Scalar,int>(nodes[i],i));
00153     weights_.push_back(weights[i]);
00154     it = points_.end();
00155   }
00156 } // end constructor
00157   
00158 template <class Scalar, class ArrayPoint, class ArrayWeight> 
00159 CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::CubatureLineSorted(
00160                    EIntrepidBurkardt rule, int numPoints, bool isNormalized ) {
00161   TEUCHOS_TEST_FOR_EXCEPTION((numPoints < 0),std::out_of_range, 
00162      ">>> ERROR (CubatureLineSorted): No rule implemented for desired number of points.");
00163   numPoints_ = numPoints;
00164   rule_type_ = rule;
00165 
00166   Teuchos::Array<Scalar> nodes(numPoints_), weights(numPoints_);
00167   if (rule==BURK_CHEBYSHEV1) { // Gauss-Chebyshev Type 1
00168     degree_ = 2*numPoints-1;
00169     IntrepidBurkardtRules::chebyshev1_compute<Scalar>(numPoints_,
00170                                         nodes.getRawPtr(),weights.getRawPtr());
00171   }
00172   else if (rule==BURK_CHEBYSHEV2) { // Gauss-Chebyshev Type 2
00173     degree_ = 2*numPoints-1;
00174     IntrepidBurkardtRules::chebyshev2_compute<Scalar>(numPoints_,
00175                                         nodes.getRawPtr(),weights.getRawPtr());
00176   }
00177   else if (rule==BURK_CLENSHAWCURTIS) { // Clenshaw-Curtis    
00178     degree_ = numPoints-1;
00179     IntrepidBurkardtRules::clenshaw_curtis_compute<Scalar>(numPoints_,
00180                                         nodes.getRawPtr(),weights.getRawPtr());
00181   }
00182   else if (rule==BURK_FEJER2) { // Fejer Type 2
00183     degree_ = numPoints-1;
00184     IntrepidBurkardtRules::fejer2_compute<Scalar>(numPoints_,
00185                                         nodes.getRawPtr(),weights.getRawPtr());
00186   }
00187   else if (rule==BURK_LEGENDRE) { // Gauss-Legendre
00188     degree_ = 2*numPoints-1;
00189     IntrepidBurkardtRules::legendre_compute<Scalar>(numPoints_,
00190                                         nodes.getRawPtr(),weights.getRawPtr());
00191   }
00192   else if (rule==BURK_PATTERSON) { // Gauss-Patterson
00193     bool correctNumPoints = false;
00194     for (int i=0; i<8; i++) {
00195       int l = (int)pow(2.0,(double)i+1.0)-1;
00196       if (numPoints==l) {
00197         correctNumPoints = true;
00198         break;
00199       }
00200     }
00201     TEUCHOS_TEST_FOR_EXCEPTION((correctNumPoints==false),std::out_of_range,
00202         ">>> ERROR (CubatureLineSorted): Number of points must be numPoints = 1, 3, 7, 15, 31, 63, 127, 255.");
00203     Scalar degree = 1.5*(double)numPoints+0.5;
00204     degree_ = (int)degree;
00205     IntrepidBurkardtRules::patterson_lookup<Scalar>(numPoints_,
00206                                         nodes.getRawPtr(),weights.getRawPtr());
00207   }
00208   else if (rule==BURK_TRAPEZOIDAL) { // Trapezoidal Rule
00209     degree_ = 2;
00210     IntrepidBurkardtRules::trapezoidal_compute<Scalar>(numPoints_,
00211                                         nodes.getRawPtr(),weights.getRawPtr());
00212   }
00213   else if (rule==BURK_HERMITE) { // Gauss-Hermite
00214     degree_ = 2*numPoints-1;
00215     IntrepidBurkardtRules::hermite_compute<Scalar>(numPoints_,
00216                                         nodes.getRawPtr(),weights.getRawPtr());
00217   }
00218   else if (rule==BURK_GENZKEISTER) { // Hermite-Genz-Keister
00219     bool correctNumPoints = false;
00220     int o_ghk[8] = {1,3,9,19,35,37,41,43};
00221     for (int i=0; i<8; i++) {
00222       if (o_ghk[i]==numPoints) {
00223         correctNumPoints = true;
00224         break;
00225       }
00226     }
00227     TEUCHOS_TEST_FOR_EXCEPTION((correctNumPoints==false),std::out_of_range,
00228        ">>> ERROR (CubatureLineSorted): Number of points must be numPoints = 1, 3, 9, 35, 37, 41, 43.");
00229     Scalar degree = 1.5*(double)numPoints+0.5;
00230     degree_ = (int)degree;
00231     IntrepidBurkardtRules::hermite_genz_keister_lookup<Scalar>(numPoints_,
00232                                         nodes.getRawPtr(),weights.getRawPtr());
00233     if (numPoints_>=37) {
00234       for (int i=0; i<numPoints_; i++) {
00235         weights[i] *= sqrt(M_PI);
00236       }
00237     }
00238   }
00239   else if (rule==BURK_LAGUERRE) { // Gauss-Laguerre
00240     degree_ = 2*numPoints-1;
00241     IntrepidBurkardtRules::laguerre_compute<Scalar>(numPoints_,
00242                                         nodes.getRawPtr(),weights.getRawPtr());
00243   }
00244   
00245   if (isNormalized) {
00246     Scalar sum = 0.0;
00247     for (int i=0; i<numPoints_; i++) {
00248       sum += weights[i];
00249     }
00250     for (int i=0; i<numPoints_; i++) {
00251       weights[i] /= sum;
00252     }
00253   }
00254   points_.clear(); weights_.clear();
00255   typename std::map<Scalar,int>::iterator it(points_.begin());
00256   for (int i=0; i<numPoints; i++) {
00257     points_.insert(it,std::pair<Scalar,int>(nodes[i],i));
00258     weights_.push_back(weights[i]);
00259     it = points_.end(); 
00260   }
00261 } // end constructor
00262   
00263 template <class Scalar, class ArrayPoint, class ArrayWeight>
00264 CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::CubatureLineSorted(
00265                  std::vector<Scalar> & points, std::vector<Scalar> & weights) {
00266 
00267   int size = (int)weights.size();
00268   TEUCHOS_TEST_FOR_EXCEPTION(((int)points.size()!=size),std::out_of_range,
00269              ">>> ERROR (CubatureLineSorted): Input dimension mismatch.");
00270   points_.clear(); weights.clear();
00271   for (int loc=0; loc<size; loc++) {
00272     points_.insert(std::pair<Scalar,int>(points[loc],loc));
00273     weights_.push_back(weights[loc]);
00274   }
00275   numPoints_ = size;
00276 }
00277 
00278 template <class Scalar, class ArrayPoint, class ArrayWeight>
00279 const char* CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::getName() const {
00280   return cubature_name_;
00281 } // end getName
00282 
00283 template <class Scalar, class ArrayPoint, class ArrayWeight>
00284 int CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::getNumPoints() const {
00285   return numPoints_;
00286 } // end getNumPoints
00287 
00288 template <class Scalar, class ArrayPoint, class ArrayWeight>
00289 void CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::getAccuracy(
00290                                            std::vector<int> & accuracy) const {
00291   accuracy.assign(1, degree_);
00292 } // end getAccuracy
00293 
00294 template <class Scalar, class ArrayPoint, class ArrayWeight>
00295 const char* CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::cubature_name_ = "INTREPID_CUBATURE_LINESORTED";
00296 
00297 template <class Scalar, class ArrayPoint, class ArrayWeight>
00298 void CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::getCubature(
00299                       ArrayPoint & cubPoints, ArrayWeight & cubWeights) const {
00300   typename std::map<Scalar,int>::const_iterator it;
00301   int i = 0;
00302   for (it = points_.begin(); it!=points_.end(); it++) {
00303     cubPoints(i)  = it->first;
00304     cubWeights(i) = weights_[it->second];
00305     i++;
00306   }
00307 } // end getCubature
00308 
00309 template <class Scalar, class ArrayPoint, class ArrayWeight>
00310 int CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::getDimension() const {
00311   return 1;
00312 } // end getDimension
00313 
00314 template <class Scalar, class ArrayPoint, class ArrayWeight>
00315 Scalar CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::getNode(
00316                                   typename std::map<Scalar,int>::iterator it) { 
00317   return it->first;
00318 } // end getNode
00319 
00320 template <class Scalar, class ArrayPoint, class ArrayWeight>
00321 Scalar CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::getWeight(int node) {
00322   return weights_[node];
00323 } // end getWeight
00324 
00325 template <class Scalar, class ArrayPoint, class ArrayWeight>
00326 Scalar CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::getWeight(
00327                                                                 Scalar point) {
00328   return weights_[points_[point]];
00329 } // end getWeight
00330 
00331 
00332 template <class Scalar, class ArrayPoint, class ArrayWeight>
00333 typename std::map<Scalar,int>::iterator CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::begin(void) {
00334   return points_.begin();
00335 } // end begin
00336 
00337 template <class Scalar, class ArrayPoint, class ArrayWeight> 
00338 typename std::map<Scalar,int>::iterator CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::end(void) {
00339   return points_.end();
00340 } // end end
00341 
00342 template <class Scalar, class ArrayPoint, class ArrayWeight>
00343 void CubatureLineSorted<Scalar,ArrayPoint,ArrayWeight>::update(
00344          Scalar alpha2, CubatureLineSorted<Scalar> & cubRule2, Scalar alpha1) {
00345 
00346   // Initialize an iterator on std::map<Scalar,Scalar>
00347   typename std::map<Scalar,int>::iterator it;
00348 
00349   // Temporary Container for updated rule
00350   typename std::map<Scalar,int> newPoints;
00351   std::vector<Scalar> newWeights;
00352  
00353   int loc = 0;
00354   Scalar node = 0.0;
00355 
00356   // Set Intersection rule1 and rule2
00357   typename std::map<Scalar,int> inter; 
00358   std::set_intersection(points_.begin(),points_.end(),
00359                         cubRule2.begin(),cubRule2.end(),
00360                         inserter(inter,inter.begin()),inter.value_comp());
00361   for (it=inter.begin(); it!=inter.end(); it++) {
00362     node = it->first;
00363     newWeights.push_back( alpha1*weights_[it->second]
00364                          +alpha2*cubRule2.getWeight(node));
00365     newPoints.insert(std::pair<Scalar,int>(node,loc));
00366     loc++;
00367   }
00368   int isize = inter.size();
00369 
00370   // Set Difference rule1 \ rule2
00371   int size = weights_.size();
00372   if (isize!=size) {
00373     typename std::map<Scalar,int> diff1; 
00374     std::set_difference(points_.begin(),points_.end(),
00375                         cubRule2.begin(),cubRule2.end(),
00376                         inserter(diff1,diff1.begin()),diff1.value_comp());
00377     for (it=diff1.begin(); it!=diff1.end(); it++) {
00378       node = it->first;
00379       newWeights.push_back(alpha1*weights_[it->second]);
00380       newPoints.insert(std::pair<Scalar,int>(node,loc));
00381       loc++;
00382     }
00383   }
00384 
00385   // Set Difference rule2 \ rule1
00386   size = cubRule2.getNumPoints();
00387   if(isize!=size) {    
00388     typename std::map<Scalar,int> diff2; 
00389     std::set_difference(cubRule2.begin(),cubRule2.end(),
00390                         points_.begin(),points_.end(),
00391                         inserter(diff2,diff2.begin()),diff2.value_comp());
00392     for (it=diff2.begin(); it!=diff2.end(); it++) {
00393       node = it->first;
00394       newWeights.push_back(alpha2*cubRule2.getWeight(it->second));
00395       newPoints.insert(std::pair<Scalar,int>(node,loc));
00396       loc++;
00397     }
00398   }
00399 
00400   points_.clear();  points_.insert(newPoints.begin(),newPoints.end());
00401   weights_.clear(); weights_.assign(newWeights.begin(),newWeights.end());
00402   numPoints_ = (int)points_.size(); 
00403 }
00404 
00405 int growthRule1D(int index, EIntrepidGrowth growth, EIntrepidBurkardt rule) {
00406   //
00407   //  Compute the growth sequence for 1D quadrature rules according to growth.
00408   //  For more information on growth rules, see 
00409   //  
00410   //  J. Burkardt. 1D Quadrature Rules For Sparse Grids.
00411   //  http://people.sc.fsu.edu/~jburkardt/presentations/sgmga_1d_rules.pdf.
00412   //  
00413   //  J. Burkardt. SGMGA: Sparse Grid Mixed Growth Anisotropic Rules.
00414   //  http://people.sc.fsu.edu/~jburkardt/cpp_src/sgmga/sgmga.html.
00415   //  
00416   //  Drew P. Kouri
00417   //  Sandia National Laboratories - CSRI
00418   //  May 27, 2011
00419   //
00420 
00421   int level = index-1;
00422   //int level = index;
00423   if (rule==BURK_CLENSHAWCURTIS) { // Clenshaw-Curtis
00424     if (growth==GROWTH_SLOWLIN) {
00425       return level+1;
00426     }
00427     else if (growth==GROWTH_SLOWLINODD) {
00428       return 2*((level+1)/2)+1;
00429     }
00430     else if (growth==GROWTH_MODLIN) {
00431       return 2*level+1;
00432     }
00433     else if (growth==GROWTH_SLOWEXP) {
00434       if (level==0) {
00435         return 1;
00436       }
00437       else { 
00438         int o = 2;
00439         while(o<2*level+1) {
00440           o = 2*(o-1)+1;
00441         }
00442         return o;
00443       }
00444     }
00445     else if (growth==GROWTH_MODEXP||growth==GROWTH_DEFAULT) {
00446       if (level==0) {
00447         return 1;
00448       }
00449       else {
00450         int o = 2;
00451         while (o<4*level+1) {
00452           o = 2*(o-1)+1;
00453         }
00454         return o;
00455       }
00456     }
00457     else if (growth==GROWTH_FULLEXP) {
00458       if (level==0) {
00459         return 1;
00460       }
00461       else {
00462         return (int)pow(2.0,(double)level)+1;
00463       }
00464     }
00465   }
00466   else if (rule==BURK_FEJER2) { // Fejer Type 2
00467     if (growth==GROWTH_SLOWLIN) {
00468       return level+1;
00469     }
00470     else if (growth==GROWTH_SLOWLINODD) {
00471       return 2*((level+1)/2)+1;
00472     }
00473     else if (growth==GROWTH_MODLIN) {
00474       return 2*level+1;
00475     }
00476     else if (growth==GROWTH_SLOWEXP) {
00477       int o = 1;
00478       while (o<2*level+1) {
00479         o = 2*o+1;
00480       }
00481       return o;
00482     }
00483     else if (growth==GROWTH_MODEXP||growth==GROWTH_DEFAULT) {
00484       int o = 1;
00485       while (o<4*level+1) {
00486         o = 2*o+1;
00487       }
00488       return o;
00489     }
00490     else if (growth==GROWTH_FULLEXP) {
00491       return (int)pow(2.0,(double)level+1.0)-1;
00492     }
00493   }
00494 
00495   else if (rule==BURK_PATTERSON) { // Gauss-Patterson
00496     if (growth==GROWTH_SLOWLIN||
00497         growth==GROWTH_SLOWLINODD||
00498         growth==GROWTH_MODLIN) {
00499       std::cout << "Specified Growth Rule Not Allowed!\n";
00500       return 0;
00501     }
00502     else if (growth==GROWTH_SLOWEXP) {
00503       if (level==0) {
00504         return 1;
00505       }
00506       else {
00507         int p = 5;
00508         int o = 3;
00509         while (p<2*level+1) {
00510           p = 2*p+1;
00511           o = 2*o+1;
00512         }
00513         return o;
00514       }
00515     }
00516     else if (growth==GROWTH_MODEXP||growth==GROWTH_DEFAULT) {
00517       if (level==0) {
00518         return 1;
00519       }
00520       else {
00521         int p = 5;
00522         int o = 3;
00523         while (p<4*level+1) {
00524           p = 2*p+1;
00525           o = 2*o+1;
00526         }
00527         return o;
00528       }
00529     }
00530     else if (growth==GROWTH_FULLEXP) {
00531       return (int)pow(2.0,(double)level+1.0)-1;
00532     }
00533   }
00534 
00535   else if (rule==BURK_LEGENDRE) { // Gauss-Legendre
00536     if (growth==GROWTH_SLOWLIN) {
00537       return level+1;
00538     }
00539     else if (growth==GROWTH_SLOWLINODD) {
00540       return 2*((level+1)/2)+1;
00541     }
00542     else if (growth==GROWTH_MODLIN||growth==GROWTH_DEFAULT) {
00543       return 2*level+1;
00544     }
00545     else if (growth==GROWTH_SLOWEXP) {
00546       int o = 1;
00547       while (2*o-1<2*level+1) {
00548         o = 2*o+1;
00549       }
00550       return o;
00551     }
00552     else if (growth==GROWTH_MODEXP) {
00553       int o = 1;
00554       while (2*o-1<4*level+1) {
00555         o = 2*o+1;
00556       }
00557       return o;
00558     }
00559     else if (growth==GROWTH_FULLEXP) {
00560       return (int)pow(2.0,(double)level+1.0)-1;
00561     }
00562   }
00563 
00564   else if (rule==BURK_HERMITE) { // Gauss-Hermite
00565     if (growth==GROWTH_SLOWLIN) {
00566       return level+1;
00567     }
00568     else if (growth==GROWTH_SLOWLINODD) {
00569       return 2*((level+1)/2)+1;
00570     }
00571     else if (growth==GROWTH_MODLIN||growth==GROWTH_DEFAULT) {
00572       return 2*level+1;
00573     }
00574     else if (growth==GROWTH_SLOWEXP) {
00575       int o = 1;
00576       while (2*o-1<2*level+1) {
00577         o = 2*o+1;
00578       }
00579       return o;
00580     }
00581     else if (growth==GROWTH_MODEXP) {
00582       int o = 1;
00583       while (2*o-1<4*level+1) {
00584         o = 2*o+1;
00585       }
00586       return o;
00587     }
00588     else if (growth==GROWTH_FULLEXP) {
00589       return (int)pow(2.0,(double)level+1.0)-1;
00590     }
00591   }
00592   
00593   else if (rule==BURK_LAGUERRE) { // Gauss-Laguerre
00594     if (growth==GROWTH_SLOWLIN) {
00595       return level+1;
00596     }
00597     else if (growth==GROWTH_SLOWLINODD) {
00598       return 2*((level+1)/2)+1;
00599     }
00600     else if (growth==GROWTH_MODLIN||growth==GROWTH_DEFAULT) {
00601       return 2*level+1;
00602     }
00603     else if (growth==GROWTH_SLOWEXP) {
00604       int o = 1;
00605       while (2*o-1<2*level+1) {
00606         o = 2*o+1;
00607       }
00608       return o;
00609     }
00610     else if (growth==GROWTH_MODEXP) {
00611       int o = 1;
00612       while (2*o-1<4*level+1) {
00613         o = 2*o+1;
00614       }
00615       return o;
00616     }
00617     else if (growth==GROWTH_FULLEXP) {
00618       return (int)pow(2.0,(double)level+1.0)-1;
00619     }
00620   }
00621 
00622   else if (rule==BURK_CHEBYSHEV1) { // Gauss-Chebyshev Type 1
00623     if (growth==GROWTH_SLOWLIN) {
00624       return level+1;
00625     }
00626     else if (growth==GROWTH_SLOWLINODD) {
00627       return 2*((level+1)/2)+1;
00628     }
00629     else if (growth==GROWTH_MODLIN||growth==GROWTH_DEFAULT) {
00630       return 2*level+1;
00631     }
00632     else if (growth==GROWTH_SLOWEXP) {
00633       int o = 1;
00634       while (2*o-1<2*level+1) {
00635         o = 2*o+1;
00636       }
00637       return o;
00638     }
00639     else if (growth==GROWTH_MODEXP) {
00640       int o = 1;
00641       while (2*o-1<4*level+1) {
00642         o = 2*o+1;
00643       }
00644       return o;
00645     }
00646     else if (growth==GROWTH_FULLEXP) {
00647       return (int)pow(2.0,(double)level+1.0)-1;
00648     }
00649   }
00650 
00651 
00652   else if (rule==BURK_CHEBYSHEV2) { // Gauss-Chebyshev Type 2
00653     if (growth==GROWTH_SLOWLIN) {
00654       return level+1;
00655     }
00656     else if (growth==GROWTH_SLOWLINODD) {
00657       return 2*((level+1)/2)+1;
00658     }
00659     else if (growth==GROWTH_MODLIN||growth==GROWTH_DEFAULT) {
00660       return 2*level+1;
00661     }
00662     else if (growth==GROWTH_SLOWEXP) {
00663       int o = 1;
00664       while (2*o-1<2*level+1) {
00665         o = 2*o+1;
00666       }
00667       return o;
00668     }
00669     else if (growth==GROWTH_MODEXP) {
00670       int o = 1;
00671       while (2*o-1<4*level+1) {
00672         o = 2*o+1;
00673       }
00674       return o;
00675     }
00676     else if (growth==GROWTH_FULLEXP) {
00677       return (int)pow(2.0,(double)level+1.0)-1;
00678     }
00679   }
00680   
00681   else if (rule==BURK_GENZKEISTER) { // Hermite-Genz-Keister  
00682     static int o_hgk[5] = { 1, 3, 9, 19, 35 };
00683     static int p_hgk[5] = { 1, 5, 15, 29, 51 };
00684     if (growth==GROWTH_SLOWLIN||
00685         growth==GROWTH_SLOWLINODD||
00686         growth==GROWTH_MODLIN) {
00687       std::cout << "Specified Growth Rule Not Allowed!\n";
00688       return 0;
00689     }
00690     else if (growth==GROWTH_SLOWEXP) { 
00691       int l = 0, p = p_hgk[l], o = o_hgk[l];
00692       while (p<2*level+1 && l<4) {
00693         l++;
00694         p = p_hgk[l];
00695         o = o_hgk[l];
00696       }
00697       return o;
00698     }
00699     else if (growth==GROWTH_MODEXP||growth==GROWTH_DEFAULT) {
00700       int l = 0, p = p_hgk[l], o = o_hgk[l];
00701       while (p<4*level+1 && l<4) {
00702         l++;
00703         p = p_hgk[l];
00704         o = o_hgk[l];
00705       }
00706       return o;
00707     }
00708     else if (growth==GROWTH_FULLEXP) {
00709       int l = level; l = std::max(l,0); l = std::min(l,4);
00710       return o_hgk[l];
00711     }
00712   }  
00713 
00714   else if (rule==BURK_TRAPEZOIDAL) { // Trapezoidal
00715     if (growth==GROWTH_SLOWLIN) {
00716       return level+1;
00717     }
00718     else if (growth==GROWTH_SLOWLINODD) {
00719       return 2*((level+1)/2)+1;
00720     }
00721     else if (growth==GROWTH_MODLIN) {
00722       return 2*level+1;
00723     }
00724     else if (growth==GROWTH_SLOWEXP) {
00725       if (level==0) {
00726         return 1;
00727       }
00728       else { 
00729         int o = 2;
00730         while(o<2*level+1) {
00731           o = 2*(o-1)+1;
00732         }
00733         return o;
00734       }
00735     }
00736     else if (growth==GROWTH_MODEXP||growth==GROWTH_DEFAULT) {
00737       if (level==0) {
00738         return 1;
00739       }
00740       else {
00741         int o = 2;
00742         while (o<4*level+1) {
00743           o = 2*(o-1)+1;
00744         }
00745         return o;
00746       }
00747     }
00748     else if (growth==GROWTH_FULLEXP) {
00749       if (level==0) {
00750         return 1;
00751       }
00752       else {
00753         return (int)pow(2.0,(double)level)+1;
00754       }
00755     }
00756   }
00757   return 0;
00758 } // end growthRule1D
00759 
00760 } // Intrepid namespace