Rythmos - Transient Integration for Differential Equations Version of the Day
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 #include "Rythmos_Types.hpp"
00034 
00035 #include "Thyra_ModelEvaluator.hpp"
00036 #include "Thyra_VectorStdOps.hpp"
00037 
00038 
00039 namespace Rythmos {
00040 
00043 template<class Scalar> 
00044 class GaussQuadrature1D : virtual public Teuchos::Describable
00045 {
00046 public:
00047 
00048   virtual RCP<const Array<Scalar> > getPoints() const =0;
00049   virtual RCP<const Array<Scalar> > getWeights() const =0;
00050   virtual RCP<const TimeRange<Scalar> > getRange() const =0;
00051   virtual int getOrder() const =0;
00052   
00053 };
00054 
00055 template<class Scalar>
00056 class GaussLegendreQuadrature1D : virtual public GaussQuadrature1D<Scalar>
00057 {
00058   public:
00059     GaussLegendreQuadrature1D(int numNodes);
00060     virtual ~GaussLegendreQuadrature1D() {}
00061     
00062     RCP<const Array<Scalar> > getPoints() const { return points_; }
00063     RCP<const Array<Scalar> > getWeights() const { return weights_; }
00064     RCP<const TimeRange<Scalar> > getRange() const { return range_; }
00065     int getOrder() const { return order_; }
00066 
00067   private:
00068     int numNodes_;
00069     void fixQuadrature_(int numNodes);
00070     RCP<Array<Scalar> > points_;
00071     RCP<Array<Scalar> > weights_;
00072     int order_;
00073     RCP<TimeRange<Scalar> > range_;
00074 };
00075 
00076 template<class Scalar>
00077 GaussLegendreQuadrature1D<Scalar>::GaussLegendreQuadrature1D(int numNodes) {
00078   fixQuadrature_(numNodes);
00079 }
00080 
00081 template<class Scalar>
00082 void GaussLegendreQuadrature1D<Scalar>::fixQuadrature_(int numNodes) {
00083   typedef Teuchos::ScalarTraits<Scalar> ST;
00084   TEUCHOS_TEST_FOR_EXCEPTION( numNodes < 2, std::out_of_range, "Error, numNodes < 2" );
00085   TEUCHOS_TEST_FOR_EXCEPTION( numNodes > 10, std::out_of_range, "Error, numNodes > 10" );
00086   numNodes_ = numNodes;
00087   range_ = Teuchos::rcp(new TimeRange<Scalar>(Scalar(-ST::one()),ST::one()));
00088   order_ = 2*numNodes_;
00089   points_ = rcp(new Array<Scalar>(numNodes_) );
00090   weights_ = rcp(new Array<Scalar>(numNodes_) );
00091 
00092   // These numbers are from David Day's matlab script
00093   if (numNodes_ == 2) {
00094     (*points_)[0] = Scalar( -0.57735026918963 );
00095     (*points_)[1] = Scalar( +0.57735026918963 );
00096     (*weights_)[0] = ST::one();
00097     (*weights_)[1] = ST::one();
00098   } else if (numNodes_ == 3) {
00099     (*points_)[0] = Scalar( -0.77459666924148 );
00100     (*points_)[1] = ST::zero();
00101     (*points_)[2] = Scalar( +0.77459666924148 );
00102     (*weights_)[0] = Scalar( 0.55555555555556 );
00103     (*weights_)[1] = Scalar( 0.88888888888889 );
00104     (*weights_)[2] = Scalar( 0.55555555555556 );
00105   } else if (numNodes_ == 4) {
00106     (*points_)[0] = Scalar( -0.86113631159405 );
00107     (*points_)[1] = Scalar( -0.33998104358486 );
00108     (*points_)[2] = Scalar( +0.33998104358486 );
00109     (*points_)[3] = Scalar( +0.86113631159405 );
00110     (*weights_)[0] = Scalar( 0.34785484513745 );
00111     (*weights_)[1] = Scalar( 0.65214515486255 );
00112     (*weights_)[2] = Scalar( 0.65214515486255 );
00113     (*weights_)[3] = Scalar( 0.34785484513745 );
00114   } else if (numNodes_ == 5) {
00115     (*points_)[0] = Scalar( -0.90617984593866 );
00116     (*points_)[1] = Scalar( -0.53846931010568 );
00117     (*points_)[2] = ST::zero();
00118     (*points_)[3] = Scalar( +0.53846931010568 );
00119     (*points_)[4] = Scalar( +0.90617984593866 );
00120     (*weights_)[0] = Scalar( 0.23692688505619 );
00121     (*weights_)[1] = Scalar( 0.47862867049937 );
00122     (*weights_)[2] = Scalar( 0.56888888888889 );
00123     (*weights_)[3] = Scalar( 0.47862867049937 );
00124     (*weights_)[4] = Scalar( 0.23692688505619 );
00125   } else if (numNodes_ == 6) {
00126     (*points_)[0] = Scalar( -0.93246951420315 );
00127     (*points_)[1] = Scalar( -0.66120938646626 );
00128     (*points_)[2] = Scalar( -0.23861918608320 );
00129     (*points_)[3] = Scalar( +0.23861918608320 );
00130     (*points_)[4] = Scalar( +0.66120938646626 );
00131     (*points_)[5] = Scalar( +0.93246951420315 );
00132     (*weights_)[0] = Scalar( 0.17132449237917 );
00133     (*weights_)[1] = Scalar( 0.36076157304814 );
00134     (*weights_)[2] = Scalar( 0.46791393457269 );
00135     (*weights_)[3] = Scalar( 0.46791393457269 );
00136     (*weights_)[4] = Scalar( 0.36076157304814 );
00137     (*weights_)[5] = Scalar( 0.17132449237917 );
00138   } else if (numNodes_ == 7) {
00139     (*points_)[0] = Scalar( -0.94910791234276 );
00140     (*points_)[1] = Scalar( -0.74153118559939 );
00141     (*points_)[2] = Scalar( -0.40584515137740 );
00142     (*points_)[3] = ST::zero();
00143     (*points_)[4] = Scalar( +0.40584515137740 );
00144     (*points_)[5] = Scalar( +0.74153118559939 );
00145     (*points_)[6] = Scalar( +0.94910791234276 );
00146     (*weights_)[0] = Scalar( 0.12948496616887 );
00147     (*weights_)[1] = Scalar( 0.27970539148928 );
00148     (*weights_)[2] = Scalar( 0.38183005050512 );
00149     (*weights_)[3] = Scalar( 0.41795918367347 );
00150     (*weights_)[4] = Scalar( 0.38183005050512 );
00151     (*weights_)[5] = Scalar( 0.27970539148928 );
00152     (*weights_)[6] = Scalar( 0.12948496616887 );
00153   } else if (numNodes_ == 8) {
00154     (*points_)[0] = Scalar( -0.96028985649754 );
00155     (*points_)[1] = Scalar( -0.79666647741363 );
00156     (*points_)[2] = Scalar( -0.52553240991633 );
00157     (*points_)[3] = Scalar( -0.18343464249565 );
00158     (*points_)[4] = Scalar( +0.18343464249565 );
00159     (*points_)[5] = Scalar( +0.52553240991633 );
00160     (*points_)[6] = Scalar( +0.79666647741363 );
00161     (*points_)[7] = Scalar( +0.96028985649754 );
00162     (*weights_)[0] = Scalar( 0.10122853629038 );
00163     (*weights_)[1] = Scalar( 0.22238103445337 );
00164     (*weights_)[2] = Scalar( 0.31370664587789 );
00165     (*weights_)[3] = Scalar( 0.36268378337836 );
00166     (*weights_)[4] = Scalar( 0.36268378337836 );
00167     (*weights_)[5] = Scalar( 0.31370664587789 );
00168     (*weights_)[6] = Scalar( 0.22238103445337 );
00169     (*weights_)[7] = Scalar( 0.10122853629038 );
00170   } else if (numNodes_ == 9) {
00171     (*points_)[0] = Scalar( -0.96816023950763 );
00172     (*points_)[1] = Scalar( -0.83603110732664 );
00173     (*points_)[2] = Scalar( -0.61337143270059 );
00174     (*points_)[3] = Scalar( -0.32425342340381 );
00175     (*points_)[4] = ST::zero();
00176     (*points_)[5] = Scalar( +0.32425342340381 );
00177     (*points_)[6] = Scalar( +0.61337143270059 );
00178     (*points_)[7] = Scalar( +0.83603110732664 );
00179     (*points_)[8] = Scalar( +0.96816023950763 );
00180     (*weights_)[0] = Scalar( 0.08127438836157 );
00181     (*weights_)[1] = Scalar( 0.18064816069486 );
00182     (*weights_)[2] = Scalar( 0.26061069640294 );
00183     (*weights_)[3] = Scalar( 0.31234707704000 );
00184     (*weights_)[4] = Scalar( 0.33023935500126 );
00185     (*weights_)[5] = Scalar( 0.31234707704000 );
00186     (*weights_)[6] = Scalar( 0.26061069640294 );
00187     (*weights_)[7] = Scalar( 0.18064816069486 );
00188     (*weights_)[8] = Scalar( 0.08127438836157 );
00189   } else if (numNodes_ == 10) {
00190     (*points_)[0] = Scalar( -0.97390652851717 );
00191     (*points_)[1] = Scalar( -0.86506336668898 );
00192     (*points_)[2] = Scalar( -0.67940956829902 );
00193     (*points_)[3] = Scalar( -0.43339539412925 );
00194     (*points_)[4] = Scalar( -0.14887433898163 );
00195     (*points_)[5] = Scalar( +0.14887433898163 );
00196     (*points_)[6] = Scalar( +0.43339539412925 );
00197     (*points_)[7] = Scalar( +0.67940956829902 );
00198     (*points_)[8] = Scalar( +0.86506336668898 );
00199     (*points_)[9] = Scalar( +0.97390652851717 );
00200     (*weights_)[0] = Scalar( 0.06667134430869 );
00201     (*weights_)[1] = Scalar( 0.14945134915058 );
00202     (*weights_)[2] = Scalar( 0.21908636251598 );
00203     (*weights_)[3] = Scalar( 0.26926671931000 );
00204     (*weights_)[4] = Scalar( 0.29552422471475 );
00205     (*weights_)[5] = Scalar( 0.29552422471475 );
00206     (*weights_)[6] = Scalar( 0.26926671931000 );
00207     (*weights_)[7] = Scalar( 0.21908636251598 );
00208     (*weights_)[8] = Scalar( 0.14945134915058 );
00209     (*weights_)[9] = Scalar( 0.06667134430869 );
00210   }
00211 }
00212 
00213 template<class Scalar>
00214 class GaussLobattoQuadrature1D : virtual public GaussQuadrature1D<Scalar>
00215 {
00216   public:
00217     GaussLobattoQuadrature1D(int numNodes);
00218     virtual ~GaussLobattoQuadrature1D() {}
00219     
00220     RCP<const Array<Scalar> > getPoints() const { return points_; }
00221     RCP<const Array<Scalar> > getWeights() const { return weights_; }
00222     RCP<const TimeRange<Scalar> > getRange() const { return range_; }
00223     int getOrder() const { return order_; }
00224 
00225   private:
00226     int numNodes_;
00227     void fixQuadrature_(int numNodes);
00228     RCP<Array<Scalar> > points_;
00229     RCP<Array<Scalar> > weights_;
00230     int order_;
00231     RCP<TimeRange<Scalar> > range_;
00232 };
00233 
00234 template<class Scalar>
00235 GaussLobattoQuadrature1D<Scalar>::GaussLobattoQuadrature1D(int numNodes) {
00236   fixQuadrature_(numNodes);
00237 }
00238 
00239 template<class Scalar>
00240 void GaussLobattoQuadrature1D<Scalar>::fixQuadrature_(int numNodes) {
00241   typedef Teuchos::ScalarTraits<Scalar> ST;
00242   TEUCHOS_TEST_FOR_EXCEPTION( numNodes < 3, std::out_of_range, "Error, numNodes < 3" );
00243   TEUCHOS_TEST_FOR_EXCEPTION( numNodes > 10, std::out_of_range, "Error, numNodes > 10" );
00244   numNodes_ = numNodes;
00245   range_ = Teuchos::rcp(new TimeRange<Scalar>(Scalar(-ST::one()),ST::one()));
00246   order_ = 2*numNodes_-2;
00247   points_ = rcp(new Array<Scalar>(numNodes_) );
00248   weights_ = rcp(new Array<Scalar>(numNodes_) );
00249 
00250   // These numbers are from David Day's matlab script
00251   if (numNodes_ == 3) {
00252     (*points_)[0] = Scalar(-ST::one());
00253     (*points_)[1] = ST::zero();
00254     (*points_)[2] = ST::one();
00255     (*weights_)[0] = Scalar( ST::one()/(3*ST::one()) ); 
00256     (*weights_)[1] = Scalar( 4*ST::one()/(3*ST::one()) );
00257     (*weights_)[2] = Scalar( ST::one()/(3*ST::one()) ); 
00258   } else if (numNodes_ == 4) {
00259     (*points_)[0] = Scalar(-ST::one());
00260     (*points_)[1] = Scalar( -0.44721359549996);
00261     (*points_)[2] = Scalar( +0.44721359549996);
00262     (*points_)[3] = ST::one();
00263     (*weights_)[0] = Scalar( ST::one()/(6*ST::one()) );
00264     (*weights_)[1] = Scalar( 5*ST::one()/(6*ST::one()) );
00265     (*weights_)[2] = Scalar( 5*ST::one()/(6*ST::one()) );
00266     (*weights_)[3] = Scalar( ST::one()/(6*ST::one()) );
00267   } else if (numNodes_ == 5) {
00268     (*points_)[0] = Scalar(-ST::one());
00269     (*points_)[1] = Scalar( -0.65465367070798 );
00270     (*points_)[2] = ST::zero();
00271     (*points_)[3] = Scalar( +0.65465367070798 );
00272     (*points_)[4] = ST::one();
00273     (*weights_)[0] = Scalar( ST::one()/(10*ST::one()) );
00274     (*weights_)[1] = Scalar( 49*ST::one()/(90*ST::one()) );
00275     (*weights_)[2] = Scalar( 32*ST::one()/(45*ST::one()) );
00276     (*weights_)[3] = Scalar( 49*ST::one()/(90*ST::one()) );
00277     (*weights_)[4] = Scalar( ST::one()/(10*ST::one()) );
00278   } else if (numNodes_ == 6) {
00279     (*points_)[0] = Scalar(-ST::one());
00280     (*points_)[1] = Scalar( -0.76505532392946 );
00281     (*points_)[2] = Scalar( -0.28523151648064 );
00282     (*points_)[3] = Scalar( +0.28523151648064 );
00283     (*points_)[4] = Scalar( +0.76505532392946 );
00284     (*points_)[5] = ST::one();
00285     (*weights_)[0] = Scalar( 0.06666666666667 );
00286     (*weights_)[1] = Scalar( 0.37847495629785 );
00287     (*weights_)[2] = Scalar( 0.55485837703549 );
00288     (*weights_)[3] = Scalar( 0.55485837703549 );
00289     (*weights_)[4] = Scalar( 0.37847495629785 );
00290     (*weights_)[5] = Scalar( 0.06666666666667 );
00291   } else if (numNodes_ == 7) {
00292     (*points_)[0] = Scalar(-ST::one());
00293     (*points_)[1] = Scalar( -0.83022389627857 );
00294     (*points_)[2] = Scalar( -0.46884879347071 );
00295     (*points_)[3] = ST::zero();
00296     (*points_)[4] = Scalar( +0.46884879347071 );
00297     (*points_)[5] = Scalar( +0.83022389627857 );
00298     (*points_)[6] = ST::one();
00299     (*weights_)[0] = Scalar( 0.04761904761905 );
00300     (*weights_)[1] = Scalar( 0.27682604736157 );
00301     (*weights_)[2] = Scalar( 0.43174538120986 );
00302     (*weights_)[3] = Scalar( 0.48761904761905 );
00303     (*weights_)[4] = Scalar( 0.43174538120986 );
00304     (*weights_)[5] = Scalar( 0.27682604736157 );
00305     (*weights_)[6] = Scalar( 0.04761904761905 );
00306   } else if (numNodes_ == 8) {
00307     (*points_)[0] = Scalar(-ST::one());
00308     (*points_)[1] = Scalar( -0.87174014850961 );
00309     (*points_)[2] = Scalar( -0.59170018143314 );
00310     (*points_)[3] = Scalar( -0.20929921790248 );
00311     (*points_)[4] = Scalar( +0.20929921790248 );
00312     (*points_)[5] = Scalar( +0.59170018143314 );
00313     (*points_)[6] = Scalar( +0.87174014850961 );
00314     (*points_)[7] = ST::one();
00315     (*weights_)[0] = Scalar( 0.03571428571429 );
00316     (*weights_)[1] = Scalar( 0.21070422714351 );
00317     (*weights_)[2] = Scalar( 0.34112269248350 );
00318     (*weights_)[3] = Scalar( 0.41245879465870 );
00319     (*weights_)[4] = Scalar( 0.41245879465870 );
00320     (*weights_)[5] = Scalar( 0.34112269248350 );
00321     (*weights_)[6] = Scalar( 0.21070422714351 );
00322     (*weights_)[7] = Scalar( 0.03571428571429 );
00323   } else if (numNodes_ == 9) {
00324     (*points_)[0] = Scalar(-ST::one());
00325     (*points_)[1] = Scalar( -0.89975799541146 );
00326     (*points_)[2] = Scalar( -0.67718627951074 );
00327     (*points_)[3] = Scalar( -0.36311746382618 );
00328     (*points_)[4] = ST::zero();
00329     (*points_)[5] = Scalar( +0.36311746382618 );
00330     (*points_)[6] = Scalar( +0.67718627951074 );
00331     (*points_)[7] = Scalar( +0.89975799541146 );
00332     (*points_)[8] = ST::one();
00333     (*weights_)[0] = Scalar( 0.02777777777778 );
00334     (*weights_)[1] = Scalar( 0.16549536156081 );
00335     (*weights_)[2] = Scalar( 0.27453871250016 );
00336     (*weights_)[3] = Scalar( 0.34642851097305 );
00337     (*weights_)[4] = Scalar( 0.37151927437642 );
00338     (*weights_)[5] = Scalar( 0.34642851097305 );
00339     (*weights_)[6] = Scalar( 0.27453871250016 );
00340     (*weights_)[7] = Scalar( 0.16549536156081 );
00341     (*weights_)[8] = Scalar( 0.02777777777778 );
00342   } else if (numNodes_ == 10) {
00343     (*points_)[0] = Scalar(-ST::one());
00344     (*points_)[1] = Scalar( -0.91953390816646 );
00345     (*points_)[2] = Scalar( -0.73877386510551 );
00346     (*points_)[3] = Scalar( -0.47792494981044 );
00347     (*points_)[4] = Scalar( -0.16527895766639 );
00348     (*points_)[5] = Scalar( +0.16527895766639 );
00349     (*points_)[6] = Scalar( +0.47792494981044 );
00350     (*points_)[7] = Scalar( +0.73877386510551 );
00351     (*points_)[8] = Scalar( +0.91953390816646 );
00352     (*points_)[9] = ST::one();
00353     (*weights_)[0] = Scalar( 0.02222222222222 );
00354     (*weights_)[1] = Scalar( 0.13330599085107 );
00355     (*weights_)[2] = Scalar( 0.22488934206313 );
00356     (*weights_)[3] = Scalar( 0.29204268367968 );
00357     (*weights_)[4] = Scalar( 0.32753976118390 );
00358     (*weights_)[5] = Scalar( 0.32753976118390 );
00359     (*weights_)[6] = Scalar( 0.29204268367968 );
00360     (*weights_)[7] = Scalar( 0.22488934206313 );
00361     (*weights_)[8] = Scalar( 0.13330599085107 );
00362     (*weights_)[9] = Scalar( 0.02222222222222 );
00363   }
00364 }
00365 
00366 
00367 template<class Scalar>
00368 RCP<Thyra::VectorBase<Scalar> > eval_f_t(
00369     const Thyra::ModelEvaluator<Scalar>& me,
00370     Scalar t
00371     ) {
00372   typedef Teuchos::ScalarTraits<Scalar> ST;
00373   typedef Thyra::ModelEvaluatorBase MEB;
00374   MEB::InArgs<Scalar> inArgs = me.createInArgs();
00375   inArgs.set_t(t);
00376   MEB::OutArgs<Scalar> outArgs = me.createOutArgs();
00377   RCP<Thyra::VectorBase<Scalar> > f_out = Thyra::createMember(me.get_f_space());
00378   V_S(outArg(*f_out),ST::zero());
00379   outArgs.set_f(f_out);
00380   me.evalModel(inArgs,outArgs);
00381   return f_out;
00382 }
00383 
00384 template<class Scalar>
00385 Scalar translateTimeRange(
00386     Scalar t,
00387     const TimeRange<Scalar>& sourceRange,
00388     const TimeRange<Scalar>& destinationRange
00389     ) {
00390   Scalar r = destinationRange.length()/sourceRange.length();
00391   return r*t+destinationRange.lower()-r*sourceRange.lower();
00392 }
00393 
00394 template<class Scalar>
00395 RCP<Thyra::VectorBase<Scalar> > computeArea(
00396     const Thyra::ModelEvaluator<Scalar>& me, 
00397     const TimeRange<Scalar>& tr, 
00398     const GaussQuadrature1D<Scalar>& gq
00399     ) {
00400   typedef Teuchos::ScalarTraits<Scalar> ST;
00401   RCP<Thyra::VectorBase<Scalar> > area = Thyra::createMember(me.get_x_space());
00402   V_S(outArg(*area),ST::zero());
00403   RCP<const TimeRange<Scalar> > sourceRange = gq.getRange();
00404   RCP<const Array<Scalar> > sourcePts = gq.getPoints();
00405   RCP<const Array<Scalar> > sourceWts = gq.getWeights();
00406   Array<Scalar> destPts(*sourcePts);
00407   for (unsigned int i=0 ; i<sourcePts->size() ; ++i) {
00408     destPts[i] = translateTimeRange<Scalar>((*sourcePts)[i],*sourceRange,tr);
00409   }
00410   Scalar r = tr.length()/sourceRange->length();
00411   for (unsigned int i=0 ; i<destPts.size() ; ++i) {
00412     RCP<Thyra::VectorBase<Scalar> > tmpVec = eval_f_t<Scalar>(me,destPts[i]);
00413     Vp_StV(outArg(*area),r*(*sourceWts)[i],*tmpVec);
00414   }
00415   return area;
00416 }
00417 
00418 
00419 } // namespace Rythmos
00420 
00421 #endif //Rythmos_QUADRATURE_BASE_H
 All Classes Functions Variables Typedefs Friends