|
Rythmos - Transient Integration for Differential Equations
Version of the Day
|
|
Rythmos - Transient Integration for Differential Equations
Version of the Day
|
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., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 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
1.7.6.1
