src

Rythmos_QuadratureBase.hpp

00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //                           Rythmos Package
00005 //                 Copyright (2006) 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 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 //@HEADER
00028 
00029 #ifndef Rythmos_QUADRATURE_BASE_H
00030 #define Rythmos_QUADRATURE_BASE_H
00031 
00032 #include "Rythmos_TimeRange.hpp"
00033 
00034 #include "Teuchos_Describable.hpp"
00035 #include "Teuchos_RCP.hpp"
00036 #include "Teuchos_Array.hpp"
00037 #include "Thyra_ModelEvaluator.hpp"
00038 
00039 
00040 namespace Rythmos {
00041 
00042 using Teuchos::RCP;
00043 using Teuchos::Array;
00044 
00047 template<class Scalar> 
00048 class GaussQuadrature1D : virtual public Teuchos::Describable
00049 {
00050 public:
00051 
00052   virtual RCP<const Array<Scalar> > getPoints() const =0;
00053   virtual RCP<const Array<Scalar> > getWeights() const =0;
00054   virtual RCP<const TimeRange<Scalar> > getRange() const =0;
00055   virtual int getOrder() const =0;
00056   
00057 };
00058 
00059 template<class Scalar>
00060 class GaussLegendreQuadrature1D : virtual public GaussQuadrature1D<Scalar>
00061 {
00062   public:
00063     GaussLegendreQuadrature1D(int numNodes);
00064     virtual ~GaussLegendreQuadrature1D() {}
00065     
00066     RCP<const Array<Scalar> > getPoints() const { return points_; }
00067     RCP<const Array<Scalar> > getWeights() const { return weights_; }
00068     RCP<const TimeRange<Scalar> > getRange() const { return range_; }
00069     int getOrder() const { return order_; }
00070 
00071   private:
00072     int numNodes_;
00073     void fixQuadrature_(int numNodes);
00074     RCP<Array<Scalar> > points_;
00075     RCP<Array<Scalar> > weights_;
00076     int order_;
00077     RCP<TimeRange<Scalar> > range_;
00078 };
00079 
00080 template<class Scalar>
00081 GaussLegendreQuadrature1D<Scalar>::GaussLegendreQuadrature1D(int numNodes) {
00082   fixQuadrature_(numNodes);
00083 }
00084 
00085 template<class Scalar>
00086 void GaussLegendreQuadrature1D<Scalar>::fixQuadrature_(int numNodes) {
00087   typedef Teuchos::ScalarTraits<Scalar> ST;
00088   TEST_FOR_EXCEPTION( numNodes < 2, std::out_of_range, "Error, numNodes < 2" );
00089   TEST_FOR_EXCEPTION( numNodes > 10, std::out_of_range, "Error, numNodes > 10" );
00090   numNodes_ = numNodes;
00091   range_ = Teuchos::rcp(new TimeRange<Scalar>(Scalar(-ST::one()),ST::one()));
00092   order_ = 2*numNodes_;
00093   points_ = rcp(new Array<Scalar>(numNodes_) );
00094   weights_ = rcp(new Array<Scalar>(numNodes_) );
00095 
00096   // These numbers are from David Day's matlab script
00097   if (numNodes_ == 2) {
00098     (*points_)[0] = Scalar( -0.57735026918963 );
00099     (*points_)[1] = Scalar( +0.57735026918963 );
00100     (*weights_)[0] = ST::one();
00101     (*weights_)[1] = ST::one();
00102   } else if (numNodes_ == 3) {
00103     (*points_)[0] = Scalar( -0.77459666924148 );
00104     (*points_)[1] = ST::zero();
00105     (*points_)[2] = Scalar( +0.77459666924148 );
00106     (*weights_)[0] = Scalar( 0.55555555555556 );
00107     (*weights_)[1] = Scalar( 0.88888888888889 );
00108     (*weights_)[2] = Scalar( 0.55555555555556 );
00109   } else if (numNodes_ == 4) {
00110     (*points_)[0] = Scalar( -0.86113631159405 );
00111     (*points_)[1] = Scalar( -0.33998104358486 );
00112     (*points_)[2] = Scalar( +0.33998104358486 );
00113     (*points_)[3] = Scalar( +0.86113631159405 );
00114     (*weights_)[0] = Scalar( 0.34785484513745 );
00115     (*weights_)[1] = Scalar( 0.65214515486255 );
00116     (*weights_)[2] = Scalar( 0.65214515486255 );
00117     (*weights_)[3] = Scalar( 0.34785484513745 );
00118   } else if (numNodes_ == 5) {
00119     (*points_)[0] = Scalar( -0.90617984593866 );
00120     (*points_)[1] = Scalar( -0.53846931010568 );
00121     (*points_)[2] = ST::zero();
00122     (*points_)[3] = Scalar( +0.53846931010568 );
00123     (*points_)[4] = Scalar( +0.90617984593866 );
00124     (*weights_)[0] = Scalar( 0.23692688505619 );
00125     (*weights_)[1] = Scalar( 0.47862867049937 );
00126     (*weights_)[2] = Scalar( 0.56888888888889 );
00127     (*weights_)[3] = Scalar( 0.47862867049937 );
00128     (*weights_)[4] = Scalar( 0.23692688505619 );
00129   } else if (numNodes_ == 6) {
00130     (*points_)[0] = Scalar( -0.93246951420315 );
00131     (*points_)[1] = Scalar( -0.66120938646626 );
00132     (*points_)[2] = Scalar( -0.23861918608320 );
00133     (*points_)[3] = Scalar( +0.23861918608320 );
00134     (*points_)[4] = Scalar( +0.66120938646626 );
00135     (*points_)[5] = Scalar( +0.93246951420315 );
00136     (*weights_)[0] = Scalar( 0.17132449237917 );
00137     (*weights_)[1] = Scalar( 0.36076157304814 );
00138     (*weights_)[2] = Scalar( 0.46791393457269 );
00139     (*weights_)[3] = Scalar( 0.46791393457269 );
00140     (*weights_)[4] = Scalar( 0.36076157304814 );
00141     (*weights_)[5] = Scalar( 0.17132449237917 );
00142   } else if (numNodes_ == 7) {
00143     (*points_)[0] = Scalar( -0.94910791234276 );
00144     (*points_)[1] = Scalar( -0.74153118559939 );
00145     (*points_)[2] = Scalar( -0.40584515137740 );
00146     (*points_)[3] = ST::zero();
00147     (*points_)[4] = Scalar( +0.40584515137740 );
00148     (*points_)[5] = Scalar( +0.74153118559939 );
00149     (*points_)[6] = Scalar( +0.94910791234276 );
00150     (*weights_)[0] = Scalar( 0.12948496616887 );
00151     (*weights_)[1] = Scalar( 0.27970539148928 );
00152     (*weights_)[2] = Scalar( 0.38183005050512 );
00153     (*weights_)[3] = Scalar( 0.41795918367347 );
00154     (*weights_)[4] = Scalar( 0.38183005050512 );
00155     (*weights_)[5] = Scalar( 0.27970539148928 );
00156     (*weights_)[6] = Scalar( 0.12948496616887 );
00157   } else if (numNodes_ == 8) {
00158     (*points_)[0] = Scalar( -0.96028985649754 );
00159     (*points_)[1] = Scalar( -0.79666647741363 );
00160     (*points_)[2] = Scalar( -0.52553240991633 );
00161     (*points_)[3] = Scalar( -0.18343464249565 );
00162     (*points_)[4] = Scalar( +0.18343464249565 );
00163     (*points_)[5] = Scalar( +0.52553240991633 );
00164     (*points_)[6] = Scalar( +0.79666647741363 );
00165     (*points_)[7] = Scalar( +0.96028985649754 );
00166     (*weights_)[0] = Scalar( 0.10122853629038 );
00167     (*weights_)[1] = Scalar( 0.22238103445337 );
00168     (*weights_)[2] = Scalar( 0.31370664587789 );
00169     (*weights_)[3] = Scalar( 0.36268378337836 );
00170     (*weights_)[4] = Scalar( 0.36268378337836 );
00171     (*weights_)[5] = Scalar( 0.31370664587789 );
00172     (*weights_)[6] = Scalar( 0.22238103445337 );
00173     (*weights_)[7] = Scalar( 0.10122853629038 );
00174   } else if (numNodes_ == 9) {
00175     (*points_)[0] = Scalar( -0.96816023950763 );
00176     (*points_)[1] = Scalar( -0.83603110732664 );
00177     (*points_)[2] = Scalar( -0.61337143270059 );
00178     (*points_)[3] = Scalar( -0.32425342340381 );
00179     (*points_)[4] = ST::zero();
00180     (*points_)[5] = Scalar( +0.32425342340381 );
00181     (*points_)[6] = Scalar( +0.61337143270059 );
00182     (*points_)[7] = Scalar( +0.83603110732664 );
00183     (*points_)[8] = Scalar( +0.96816023950763 );
00184     (*weights_)[0] = Scalar( 0.08127438836157 );
00185     (*weights_)[1] = Scalar( 0.18064816069486 );
00186     (*weights_)[2] = Scalar( 0.26061069640294 );
00187     (*weights_)[3] = Scalar( 0.31234707704000 );
00188     (*weights_)[4] = Scalar( 0.33023935500126 );
00189     (*weights_)[5] = Scalar( 0.31234707704000 );
00190     (*weights_)[6] = Scalar( 0.26061069640294 );
00191     (*weights_)[7] = Scalar( 0.18064816069486 );
00192     (*weights_)[8] = Scalar( 0.08127438836157 );
00193   } else if (numNodes_ == 10) {
00194     (*points_)[0] = Scalar( -0.97390652851717 );
00195     (*points_)[1] = Scalar( -0.86506336668898 );
00196     (*points_)[2] = Scalar( -0.67940956829902 );
00197     (*points_)[3] = Scalar( -0.43339539412925 );
00198     (*points_)[4] = Scalar( -0.14887433898163 );
00199     (*points_)[5] = Scalar( +0.14887433898163 );
00200     (*points_)[6] = Scalar( +0.43339539412925 );
00201     (*points_)[7] = Scalar( +0.67940956829902 );
00202     (*points_)[8] = Scalar( +0.86506336668898 );
00203     (*points_)[9] = Scalar( +0.97390652851717 );
00204     (*weights_)[0] = Scalar( 0.06667134430869 );
00205     (*weights_)[1] = Scalar( 0.14945134915058 );
00206     (*weights_)[2] = Scalar( 0.21908636251598 );
00207     (*weights_)[3] = Scalar( 0.26926671931000 );
00208     (*weights_)[4] = Scalar( 0.29552422471475 );
00209     (*weights_)[5] = Scalar( 0.29552422471475 );
00210     (*weights_)[6] = Scalar( 0.26926671931000 );
00211     (*weights_)[7] = Scalar( 0.21908636251598 );
00212     (*weights_)[8] = Scalar( 0.14945134915058 );
00213     (*weights_)[9] = Scalar( 0.06667134430869 );
00214   }
00215 }
00216 
00217 template<class Scalar>
00218 class GaussLobattoQuadrature1D : virtual public GaussQuadrature1D<Scalar>
00219 {
00220   public:
00221     GaussLobattoQuadrature1D(int numNodes);
00222     virtual ~GaussLobattoQuadrature1D() {}
00223     
00224     RCP<const Array<Scalar> > getPoints() const { return points_; }
00225     RCP<const Array<Scalar> > getWeights() const { return weights_; }
00226     RCP<const TimeRange<Scalar> > getRange() const { return range_; }
00227     int getOrder() const { return order_; }
00228 
00229   private:
00230     int numNodes_;
00231     void fixQuadrature_(int numNodes);
00232     RCP<Array<Scalar> > points_;
00233     RCP<Array<Scalar> > weights_;
00234     int order_;
00235     RCP<TimeRange<Scalar> > range_;
00236 };
00237 
00238 template<class Scalar>
00239 GaussLobattoQuadrature1D<Scalar>::GaussLobattoQuadrature1D(int numNodes) {
00240   fixQuadrature_(numNodes);
00241 }
00242 
00243 template<class Scalar>
00244 void GaussLobattoQuadrature1D<Scalar>::fixQuadrature_(int numNodes) {
00245   typedef Teuchos::ScalarTraits<Scalar> ST;
00246   TEST_FOR_EXCEPTION( numNodes < 3, std::out_of_range, "Error, numNodes < 3" );
00247   TEST_FOR_EXCEPTION( numNodes > 10, std::out_of_range, "Error, numNodes > 10" );
00248   numNodes_ = numNodes;
00249   range_ = Teuchos::rcp(new TimeRange<Scalar>(Scalar(-ST::one()),ST::one()));
00250   order_ = 2*numNodes_-2;
00251   points_ = rcp(new Array<Scalar>(numNodes_) );
00252   weights_ = rcp(new Array<Scalar>(numNodes_) );
00253 
00254   // These numbers are from David Day's matlab script
00255   if (numNodes_ == 3) {
00256     (*points_)[0] = Scalar(-ST::one());
00257     (*points_)[1] = ST::zero();
00258     (*points_)[2] = ST::one();
00259     (*weights_)[0] = Scalar( ST::one()/(3*ST::one()) ); 
00260     (*weights_)[1] = Scalar( 4*ST::one()/(3*ST::one()) );
00261     (*weights_)[2] = Scalar( ST::one()/(3*ST::one()) ); 
00262   } else if (numNodes_ == 4) {
00263     (*points_)[0] = Scalar(-ST::one());
00264     (*points_)[1] = Scalar( -0.44721359549996);
00265     (*points_)[2] = Scalar( +0.44721359549996);
00266     (*points_)[3] = ST::one();
00267     (*weights_)[0] = Scalar( ST::one()/(6*ST::one()) );
00268     (*weights_)[1] = Scalar( 5*ST::one()/(6*ST::one()) );
00269     (*weights_)[2] = Scalar( 5*ST::one()/(6*ST::one()) );
00270     (*weights_)[3] = Scalar( ST::one()/(6*ST::one()) );
00271   } else if (numNodes_ == 5) {
00272     (*points_)[0] = Scalar(-ST::one());
00273     (*points_)[1] = Scalar( -0.65465367070798 );
00274     (*points_)[2] = ST::zero();
00275     (*points_)[3] = Scalar( +0.65465367070798 );
00276     (*points_)[4] = ST::one();
00277     (*weights_)[0] = Scalar( ST::one()/(10*ST::one()) );
00278     (*weights_)[1] = Scalar( 49*ST::one()/(90*ST::one()) );
00279     (*weights_)[2] = Scalar( 32*ST::one()/(45*ST::one()) );
00280     (*weights_)[3] = Scalar( 49*ST::one()/(90*ST::one()) );
00281     (*weights_)[4] = Scalar( ST::one()/(10*ST::one()) );
00282   } else if (numNodes_ == 6) {
00283     (*points_)[0] = Scalar(-ST::one());
00284     (*points_)[1] = Scalar( -0.76505532392946 );
00285     (*points_)[2] = Scalar( -0.28523151648064 );
00286     (*points_)[3] = Scalar( +0.28523151648064 );
00287     (*points_)[4] = Scalar( +0.76505532392946 );
00288     (*points_)[5] = ST::one();
00289     (*weights_)[0] = Scalar( 0.06666666666667 );
00290     (*weights_)[1] = Scalar( 0.37847495629785 );
00291     (*weights_)[2] = Scalar( 0.55485837703549 );
00292     (*weights_)[3] = Scalar( 0.55485837703549 );
00293     (*weights_)[4] = Scalar( 0.37847495629785 );
00294     (*weights_)[5] = Scalar( 0.06666666666667 );
00295   } else if (numNodes_ == 7) {
00296     (*points_)[0] = Scalar(-ST::one());
00297     (*points_)[1] = Scalar( -0.83022389627857 );
00298     (*points_)[2] = Scalar( -0.46884879347071 );
00299     (*points_)[3] = ST::zero();
00300     (*points_)[4] = Scalar( +0.46884879347071 );
00301     (*points_)[5] = Scalar( +0.83022389627857 );
00302     (*points_)[6] = ST::one();
00303     (*weights_)[0] = Scalar( 0.04761904761905 );
00304     (*weights_)[1] = Scalar( 0.27682604736157 );
00305     (*weights_)[2] = Scalar( 0.43174538120986 );
00306     (*weights_)[3] = Scalar( 0.48761904761905 );
00307     (*weights_)[4] = Scalar( 0.43174538120986 );
00308     (*weights_)[5] = Scalar( 0.27682604736157 );
00309     (*weights_)[6] = Scalar( 0.04761904761905 );
00310   } else if (numNodes_ == 8) {
00311     (*points_)[0] = Scalar(-ST::one());
00312     (*points_)[1] = Scalar( -0.87174014850961 );
00313     (*points_)[2] = Scalar( -0.59170018143314 );
00314     (*points_)[3] = Scalar( -0.20929921790248 );
00315     (*points_)[4] = Scalar( +0.20929921790248 );
00316     (*points_)[5] = Scalar( +0.59170018143314 );
00317     (*points_)[6] = Scalar( +0.87174014850961 );
00318     (*points_)[7] = ST::one();
00319     (*weights_)[0] = Scalar( 0.03571428571429 );
00320     (*weights_)[1] = Scalar( 0.21070422714351 );
00321     (*weights_)[2] = Scalar( 0.34112269248350 );
00322     (*weights_)[3] = Scalar( 0.41245879465870 );
00323     (*weights_)[4] = Scalar( 0.41245879465870 );
00324     (*weights_)[5] = Scalar( 0.34112269248350 );
00325     (*weights_)[6] = Scalar( 0.21070422714351 );
00326     (*weights_)[7] = Scalar( 0.03571428571429 );
00327   } else if (numNodes_ == 9) {
00328     (*points_)[0] = Scalar(-ST::one());
00329     (*points_)[1] = Scalar( -0.89975799541146 );
00330     (*points_)[2] = Scalar( -0.67718627951074 );
00331     (*points_)[3] = Scalar( -0.36311746382618 );
00332     (*points_)[4] = ST::zero();
00333     (*points_)[5] = Scalar( +0.36311746382618 );
00334     (*points_)[6] = Scalar( +0.67718627951074 );
00335     (*points_)[7] = Scalar( +0.89975799541146 );
00336     (*points_)[8] = ST::one();
00337     (*weights_)[0] = Scalar( 0.02777777777778 );
00338     (*weights_)[1] = Scalar( 0.16549536156081 );
00339     (*weights_)[2] = Scalar( 0.27453871250016 );
00340     (*weights_)[3] = Scalar( 0.34642851097305 );
00341     (*weights_)[4] = Scalar( 0.37151927437642 );
00342     (*weights_)[5] = Scalar( 0.34642851097305 );
00343     (*weights_)[6] = Scalar( 0.27453871250016 );
00344     (*weights_)[7] = Scalar( 0.16549536156081 );
00345     (*weights_)[8] = Scalar( 0.02777777777778 );
00346   } else if (numNodes_ == 10) {
00347     (*points_)[0] = Scalar(-ST::one());
00348     (*points_)[1] = Scalar( -0.91953390816646 );
00349     (*points_)[2] = Scalar( -0.73877386510551 );
00350     (*points_)[3] = Scalar( -0.47792494981044 );
00351     (*points_)[4] = Scalar( -0.16527895766639 );
00352     (*points_)[5] = Scalar( +0.16527895766639 );
00353     (*points_)[6] = Scalar( +0.47792494981044 );
00354     (*points_)[7] = Scalar( +0.73877386510551 );
00355     (*points_)[8] = Scalar( +0.91953390816646 );
00356     (*points_)[9] = ST::one();
00357     (*weights_)[0] = Scalar( 0.02222222222222 );
00358     (*weights_)[1] = Scalar( 0.13330599085107 );
00359     (*weights_)[2] = Scalar( 0.22488934206313 );
00360     (*weights_)[3] = Scalar( 0.29204268367968 );
00361     (*weights_)[4] = Scalar( 0.32753976118390 );
00362     (*weights_)[5] = Scalar( 0.32753976118390 );
00363     (*weights_)[6] = Scalar( 0.29204268367968 );
00364     (*weights_)[7] = Scalar( 0.22488934206313 );
00365     (*weights_)[8] = Scalar( 0.13330599085107 );
00366     (*weights_)[9] = Scalar( 0.02222222222222 );
00367   }
00368 }
00369 
00370 
00371 template<class Scalar>
00372 RCP<Thyra::VectorBase<Scalar> > eval_f_t(
00373     const Thyra::ModelEvaluator<Scalar>& me,
00374     Scalar t
00375     ) {
00376   typedef Teuchos::ScalarTraits<Scalar> ST;
00377   typedef Thyra::ModelEvaluatorBase MEB;
00378   MEB::InArgs<Scalar> inArgs = me.createInArgs();
00379   inArgs.set_t(t);
00380   MEB::OutArgs<Scalar> outArgs = me.createOutArgs();
00381   RCP<Thyra::VectorBase<Scalar> > f_out = Thyra::createMember(me.get_f_space());
00382   V_S(outArg(*f_out),ST::zero());
00383   outArgs.set_f(f_out);
00384   me.evalModel(inArgs,outArgs);
00385   return f_out;
00386 }
00387 
00388 template<class Scalar>
00389 Scalar translateTimeRange(
00390     Scalar t,
00391     const TimeRange<Scalar>& sourceRange,
00392     const TimeRange<Scalar>& destinationRange
00393     ) {
00394   Scalar r = destinationRange.length()/sourceRange.length();
00395   return r*t+destinationRange.lower()-r*sourceRange.lower();
00396 }
00397 
00398 template<class Scalar>
00399 RCP<Thyra::VectorBase<Scalar> > computeArea(
00400     const Thyra::ModelEvaluator<Scalar>& me, 
00401     const TimeRange<Scalar>& tr, 
00402     const GaussQuadrature1D<Scalar>& gq
00403     ) {
00404   typedef Teuchos::ScalarTraits<Scalar> ST;
00405   RCP<Thyra::VectorBase<Scalar> > area = Thyra::createMember(me.get_x_space());
00406   V_S(outArg(*area),ST::zero());
00407   RCP<const TimeRange<Scalar> > sourceRange = gq.getRange();
00408   RCP<const Array<Scalar> > sourcePts = gq.getPoints();
00409   RCP<const Array<Scalar> > sourceWts = gq.getWeights();
00410   Array<Scalar> destPts(*sourcePts);
00411   for (unsigned int i=0 ; i<sourcePts->size() ; ++i) {
00412     destPts[i] = translateTimeRange<Scalar>((*sourcePts)[i],*sourceRange,tr);
00413   }
00414   Scalar r = tr.length()/sourceRange->length();
00415   for (unsigned int i=0 ; i<destPts.size() ; ++i) {
00416     RCP<Thyra::VectorBase<Scalar> > tmpVec = eval_f_t<Scalar>(me,destPts[i]);
00417     Vp_StV(outArg(*area),r*(*sourceWts)[i],*tmpVec);
00418   }
00419   return area;
00420 }
00421 
00422 
00423 } // namespace Rythmos
00424 
00425 #endif //Rythmos_QUADRATURE_BASE_H
Generated on Wed May 12 21:25:43 2010 for Rythmos - Transient Integration for Differential Equations by  doxygen 1.4.7