Rythmos - Transient Integration for Differential Equations  Version of the Day
 All Classes Functions Variables Typedefs Pages
Rythmos_QuadratureBase.hpp
1 //@HEADER
2 // ***********************************************************************
3 //
4 // Rythmos Package
5 // Copyright (2006) 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 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23 // USA
24 // Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
25 //
26 // ***********************************************************************
27 //@HEADER
28 
29 #ifndef Rythmos_QUADRATURE_BASE_H
30 #define Rythmos_QUADRATURE_BASE_H
31 
32 #include "Rythmos_TimeRange.hpp"
33 #include "Rythmos_Types.hpp"
34 
35 #include "Thyra_ModelEvaluator.hpp"
36 #include "Thyra_VectorStdOps.hpp"
37 
38 
39 namespace Rythmos {
40 
43 template<class Scalar>
44 class GaussQuadrature1D : virtual public Teuchos::Describable
45 {
46 public:
47 
48  virtual RCP<const Array<Scalar> > getPoints() const =0;
49  virtual RCP<const Array<Scalar> > getWeights() const =0;
50  virtual RCP<const TimeRange<Scalar> > getRange() const =0;
51  virtual int getOrder() const =0;
52 
53 };
54 
55 template<class Scalar>
56 class GaussLegendreQuadrature1D : virtual public GaussQuadrature1D<Scalar>
57 {
58  public:
59  GaussLegendreQuadrature1D(int numNodes);
60  virtual ~GaussLegendreQuadrature1D() {}
61 
62  RCP<const Array<Scalar> > getPoints() const { return points_; }
63  RCP<const Array<Scalar> > getWeights() const { return weights_; }
64  RCP<const TimeRange<Scalar> > getRange() const { return range_; }
65  int getOrder() const { return order_; }
66 
67  private:
68  int numNodes_;
69  void fixQuadrature_(int numNodes);
70  RCP<Array<Scalar> > points_;
71  RCP<Array<Scalar> > weights_;
72  int order_;
73  RCP<TimeRange<Scalar> > range_;
74 };
75 
76 template<class Scalar>
77 GaussLegendreQuadrature1D<Scalar>::GaussLegendreQuadrature1D(int numNodes) {
78  fixQuadrature_(numNodes);
79 }
80 
81 template<class Scalar>
82 void GaussLegendreQuadrature1D<Scalar>::fixQuadrature_(int numNodes) {
83  typedef Teuchos::ScalarTraits<Scalar> ST;
84  TEUCHOS_TEST_FOR_EXCEPTION( numNodes < 2, std::out_of_range, "Error, numNodes < 2" );
85  TEUCHOS_TEST_FOR_EXCEPTION( numNodes > 10, std::out_of_range, "Error, numNodes > 10" );
86  numNodes_ = numNodes;
87  range_ = Teuchos::rcp(new TimeRange<Scalar>(Scalar(-ST::one()),ST::one()));
88  order_ = 2*numNodes_;
89  points_ = rcp(new Array<Scalar>(numNodes_) );
90  weights_ = rcp(new Array<Scalar>(numNodes_) );
91 
92  // These numbers are from David Day's matlab script
93  if (numNodes_ == 2) {
94  (*points_)[0] = Scalar( -0.57735026918963 );
95  (*points_)[1] = Scalar( +0.57735026918963 );
96  (*weights_)[0] = ST::one();
97  (*weights_)[1] = ST::one();
98  } else if (numNodes_ == 3) {
99  (*points_)[0] = Scalar( -0.77459666924148 );
100  (*points_)[1] = ST::zero();
101  (*points_)[2] = Scalar( +0.77459666924148 );
102  (*weights_)[0] = Scalar( 0.55555555555556 );
103  (*weights_)[1] = Scalar( 0.88888888888889 );
104  (*weights_)[2] = Scalar( 0.55555555555556 );
105  } else if (numNodes_ == 4) {
106  (*points_)[0] = Scalar( -0.86113631159405 );
107  (*points_)[1] = Scalar( -0.33998104358486 );
108  (*points_)[2] = Scalar( +0.33998104358486 );
109  (*points_)[3] = Scalar( +0.86113631159405 );
110  (*weights_)[0] = Scalar( 0.34785484513745 );
111  (*weights_)[1] = Scalar( 0.65214515486255 );
112  (*weights_)[2] = Scalar( 0.65214515486255 );
113  (*weights_)[3] = Scalar( 0.34785484513745 );
114  } else if (numNodes_ == 5) {
115  (*points_)[0] = Scalar( -0.90617984593866 );
116  (*points_)[1] = Scalar( -0.53846931010568 );
117  (*points_)[2] = ST::zero();
118  (*points_)[3] = Scalar( +0.53846931010568 );
119  (*points_)[4] = Scalar( +0.90617984593866 );
120  (*weights_)[0] = Scalar( 0.23692688505619 );
121  (*weights_)[1] = Scalar( 0.47862867049937 );
122  (*weights_)[2] = Scalar( 0.56888888888889 );
123  (*weights_)[3] = Scalar( 0.47862867049937 );
124  (*weights_)[4] = Scalar( 0.23692688505619 );
125  } else if (numNodes_ == 6) {
126  (*points_)[0] = Scalar( -0.93246951420315 );
127  (*points_)[1] = Scalar( -0.66120938646626 );
128  (*points_)[2] = Scalar( -0.23861918608320 );
129  (*points_)[3] = Scalar( +0.23861918608320 );
130  (*points_)[4] = Scalar( +0.66120938646626 );
131  (*points_)[5] = Scalar( +0.93246951420315 );
132  (*weights_)[0] = Scalar( 0.17132449237917 );
133  (*weights_)[1] = Scalar( 0.36076157304814 );
134  (*weights_)[2] = Scalar( 0.46791393457269 );
135  (*weights_)[3] = Scalar( 0.46791393457269 );
136  (*weights_)[4] = Scalar( 0.36076157304814 );
137  (*weights_)[5] = Scalar( 0.17132449237917 );
138  } else if (numNodes_ == 7) {
139  (*points_)[0] = Scalar( -0.94910791234276 );
140  (*points_)[1] = Scalar( -0.74153118559939 );
141  (*points_)[2] = Scalar( -0.40584515137740 );
142  (*points_)[3] = ST::zero();
143  (*points_)[4] = Scalar( +0.40584515137740 );
144  (*points_)[5] = Scalar( +0.74153118559939 );
145  (*points_)[6] = Scalar( +0.94910791234276 );
146  (*weights_)[0] = Scalar( 0.12948496616887 );
147  (*weights_)[1] = Scalar( 0.27970539148928 );
148  (*weights_)[2] = Scalar( 0.38183005050512 );
149  (*weights_)[3] = Scalar( 0.41795918367347 );
150  (*weights_)[4] = Scalar( 0.38183005050512 );
151  (*weights_)[5] = Scalar( 0.27970539148928 );
152  (*weights_)[6] = Scalar( 0.12948496616887 );
153  } else if (numNodes_ == 8) {
154  (*points_)[0] = Scalar( -0.96028985649754 );
155  (*points_)[1] = Scalar( -0.79666647741363 );
156  (*points_)[2] = Scalar( -0.52553240991633 );
157  (*points_)[3] = Scalar( -0.18343464249565 );
158  (*points_)[4] = Scalar( +0.18343464249565 );
159  (*points_)[5] = Scalar( +0.52553240991633 );
160  (*points_)[6] = Scalar( +0.79666647741363 );
161  (*points_)[7] = Scalar( +0.96028985649754 );
162  (*weights_)[0] = Scalar( 0.10122853629038 );
163  (*weights_)[1] = Scalar( 0.22238103445337 );
164  (*weights_)[2] = Scalar( 0.31370664587789 );
165  (*weights_)[3] = Scalar( 0.36268378337836 );
166  (*weights_)[4] = Scalar( 0.36268378337836 );
167  (*weights_)[5] = Scalar( 0.31370664587789 );
168  (*weights_)[6] = Scalar( 0.22238103445337 );
169  (*weights_)[7] = Scalar( 0.10122853629038 );
170  } else if (numNodes_ == 9) {
171  (*points_)[0] = Scalar( -0.96816023950763 );
172  (*points_)[1] = Scalar( -0.83603110732664 );
173  (*points_)[2] = Scalar( -0.61337143270059 );
174  (*points_)[3] = Scalar( -0.32425342340381 );
175  (*points_)[4] = ST::zero();
176  (*points_)[5] = Scalar( +0.32425342340381 );
177  (*points_)[6] = Scalar( +0.61337143270059 );
178  (*points_)[7] = Scalar( +0.83603110732664 );
179  (*points_)[8] = Scalar( +0.96816023950763 );
180  (*weights_)[0] = Scalar( 0.08127438836157 );
181  (*weights_)[1] = Scalar( 0.18064816069486 );
182  (*weights_)[2] = Scalar( 0.26061069640294 );
183  (*weights_)[3] = Scalar( 0.31234707704000 );
184  (*weights_)[4] = Scalar( 0.33023935500126 );
185  (*weights_)[5] = Scalar( 0.31234707704000 );
186  (*weights_)[6] = Scalar( 0.26061069640294 );
187  (*weights_)[7] = Scalar( 0.18064816069486 );
188  (*weights_)[8] = Scalar( 0.08127438836157 );
189  } else if (numNodes_ == 10) {
190  (*points_)[0] = Scalar( -0.97390652851717 );
191  (*points_)[1] = Scalar( -0.86506336668898 );
192  (*points_)[2] = Scalar( -0.67940956829902 );
193  (*points_)[3] = Scalar( -0.43339539412925 );
194  (*points_)[4] = Scalar( -0.14887433898163 );
195  (*points_)[5] = Scalar( +0.14887433898163 );
196  (*points_)[6] = Scalar( +0.43339539412925 );
197  (*points_)[7] = Scalar( +0.67940956829902 );
198  (*points_)[8] = Scalar( +0.86506336668898 );
199  (*points_)[9] = Scalar( +0.97390652851717 );
200  (*weights_)[0] = Scalar( 0.06667134430869 );
201  (*weights_)[1] = Scalar( 0.14945134915058 );
202  (*weights_)[2] = Scalar( 0.21908636251598 );
203  (*weights_)[3] = Scalar( 0.26926671931000 );
204  (*weights_)[4] = Scalar( 0.29552422471475 );
205  (*weights_)[5] = Scalar( 0.29552422471475 );
206  (*weights_)[6] = Scalar( 0.26926671931000 );
207  (*weights_)[7] = Scalar( 0.21908636251598 );
208  (*weights_)[8] = Scalar( 0.14945134915058 );
209  (*weights_)[9] = Scalar( 0.06667134430869 );
210  }
211 }
212 
213 template<class Scalar>
214 class GaussLobattoQuadrature1D : virtual public GaussQuadrature1D<Scalar>
215 {
216  public:
217  GaussLobattoQuadrature1D(int numNodes);
218  virtual ~GaussLobattoQuadrature1D() {}
219 
220  RCP<const Array<Scalar> > getPoints() const { return points_; }
221  RCP<const Array<Scalar> > getWeights() const { return weights_; }
222  RCP<const TimeRange<Scalar> > getRange() const { return range_; }
223  int getOrder() const { return order_; }
224 
225  private:
226  int numNodes_;
227  void fixQuadrature_(int numNodes);
228  RCP<Array<Scalar> > points_;
229  RCP<Array<Scalar> > weights_;
230  int order_;
231  RCP<TimeRange<Scalar> > range_;
232 };
233 
234 template<class Scalar>
235 GaussLobattoQuadrature1D<Scalar>::GaussLobattoQuadrature1D(int numNodes) {
236  fixQuadrature_(numNodes);
237 }
238 
239 template<class Scalar>
240 void GaussLobattoQuadrature1D<Scalar>::fixQuadrature_(int numNodes) {
241  typedef Teuchos::ScalarTraits<Scalar> ST;
242  TEUCHOS_TEST_FOR_EXCEPTION( numNodes < 3, std::out_of_range, "Error, numNodes < 3" );
243  TEUCHOS_TEST_FOR_EXCEPTION( numNodes > 10, std::out_of_range, "Error, numNodes > 10" );
244  numNodes_ = numNodes;
245  range_ = Teuchos::rcp(new TimeRange<Scalar>(Scalar(-ST::one()),ST::one()));
246  order_ = 2*numNodes_-2;
247  points_ = rcp(new Array<Scalar>(numNodes_) );
248  weights_ = rcp(new Array<Scalar>(numNodes_) );
249 
250  // These numbers are from David Day's matlab script
251  if (numNodes_ == 3) {
252  (*points_)[0] = Scalar(-ST::one());
253  (*points_)[1] = ST::zero();
254  (*points_)[2] = ST::one();
255  (*weights_)[0] = Scalar( ST::one()/(3*ST::one()) );
256  (*weights_)[1] = Scalar( 4*ST::one()/(3*ST::one()) );
257  (*weights_)[2] = Scalar( ST::one()/(3*ST::one()) );
258  } else if (numNodes_ == 4) {
259  (*points_)[0] = Scalar(-ST::one());
260  (*points_)[1] = Scalar( -0.44721359549996);
261  (*points_)[2] = Scalar( +0.44721359549996);
262  (*points_)[3] = ST::one();
263  (*weights_)[0] = Scalar( ST::one()/(6*ST::one()) );
264  (*weights_)[1] = Scalar( 5*ST::one()/(6*ST::one()) );
265  (*weights_)[2] = Scalar( 5*ST::one()/(6*ST::one()) );
266  (*weights_)[3] = Scalar( ST::one()/(6*ST::one()) );
267  } else if (numNodes_ == 5) {
268  (*points_)[0] = Scalar(-ST::one());
269  (*points_)[1] = Scalar( -0.65465367070798 );
270  (*points_)[2] = ST::zero();
271  (*points_)[3] = Scalar( +0.65465367070798 );
272  (*points_)[4] = ST::one();
273  (*weights_)[0] = Scalar( ST::one()/(10*ST::one()) );
274  (*weights_)[1] = Scalar( 49*ST::one()/(90*ST::one()) );
275  (*weights_)[2] = Scalar( 32*ST::one()/(45*ST::one()) );
276  (*weights_)[3] = Scalar( 49*ST::one()/(90*ST::one()) );
277  (*weights_)[4] = Scalar( ST::one()/(10*ST::one()) );
278  } else if (numNodes_ == 6) {
279  (*points_)[0] = Scalar(-ST::one());
280  (*points_)[1] = Scalar( -0.76505532392946 );
281  (*points_)[2] = Scalar( -0.28523151648064 );
282  (*points_)[3] = Scalar( +0.28523151648064 );
283  (*points_)[4] = Scalar( +0.76505532392946 );
284  (*points_)[5] = ST::one();
285  (*weights_)[0] = Scalar( 0.06666666666667 );
286  (*weights_)[1] = Scalar( 0.37847495629785 );
287  (*weights_)[2] = Scalar( 0.55485837703549 );
288  (*weights_)[3] = Scalar( 0.55485837703549 );
289  (*weights_)[4] = Scalar( 0.37847495629785 );
290  (*weights_)[5] = Scalar( 0.06666666666667 );
291  } else if (numNodes_ == 7) {
292  (*points_)[0] = Scalar(-ST::one());
293  (*points_)[1] = Scalar( -0.83022389627857 );
294  (*points_)[2] = Scalar( -0.46884879347071 );
295  (*points_)[3] = ST::zero();
296  (*points_)[4] = Scalar( +0.46884879347071 );
297  (*points_)[5] = Scalar( +0.83022389627857 );
298  (*points_)[6] = ST::one();
299  (*weights_)[0] = Scalar( 0.04761904761905 );
300  (*weights_)[1] = Scalar( 0.27682604736157 );
301  (*weights_)[2] = Scalar( 0.43174538120986 );
302  (*weights_)[3] = Scalar( 0.48761904761905 );
303  (*weights_)[4] = Scalar( 0.43174538120986 );
304  (*weights_)[5] = Scalar( 0.27682604736157 );
305  (*weights_)[6] = Scalar( 0.04761904761905 );
306  } else if (numNodes_ == 8) {
307  (*points_)[0] = Scalar(-ST::one());
308  (*points_)[1] = Scalar( -0.87174014850961 );
309  (*points_)[2] = Scalar( -0.59170018143314 );
310  (*points_)[3] = Scalar( -0.20929921790248 );
311  (*points_)[4] = Scalar( +0.20929921790248 );
312  (*points_)[5] = Scalar( +0.59170018143314 );
313  (*points_)[6] = Scalar( +0.87174014850961 );
314  (*points_)[7] = ST::one();
315  (*weights_)[0] = Scalar( 0.03571428571429 );
316  (*weights_)[1] = Scalar( 0.21070422714351 );
317  (*weights_)[2] = Scalar( 0.34112269248350 );
318  (*weights_)[3] = Scalar( 0.41245879465870 );
319  (*weights_)[4] = Scalar( 0.41245879465870 );
320  (*weights_)[5] = Scalar( 0.34112269248350 );
321  (*weights_)[6] = Scalar( 0.21070422714351 );
322  (*weights_)[7] = Scalar( 0.03571428571429 );
323  } else if (numNodes_ == 9) {
324  (*points_)[0] = Scalar(-ST::one());
325  (*points_)[1] = Scalar( -0.89975799541146 );
326  (*points_)[2] = Scalar( -0.67718627951074 );
327  (*points_)[3] = Scalar( -0.36311746382618 );
328  (*points_)[4] = ST::zero();
329  (*points_)[5] = Scalar( +0.36311746382618 );
330  (*points_)[6] = Scalar( +0.67718627951074 );
331  (*points_)[7] = Scalar( +0.89975799541146 );
332  (*points_)[8] = ST::one();
333  (*weights_)[0] = Scalar( 0.02777777777778 );
334  (*weights_)[1] = Scalar( 0.16549536156081 );
335  (*weights_)[2] = Scalar( 0.27453871250016 );
336  (*weights_)[3] = Scalar( 0.34642851097305 );
337  (*weights_)[4] = Scalar( 0.37151927437642 );
338  (*weights_)[5] = Scalar( 0.34642851097305 );
339  (*weights_)[6] = Scalar( 0.27453871250016 );
340  (*weights_)[7] = Scalar( 0.16549536156081 );
341  (*weights_)[8] = Scalar( 0.02777777777778 );
342  } else if (numNodes_ == 10) {
343  (*points_)[0] = Scalar(-ST::one());
344  (*points_)[1] = Scalar( -0.91953390816646 );
345  (*points_)[2] = Scalar( -0.73877386510551 );
346  (*points_)[3] = Scalar( -0.47792494981044 );
347  (*points_)[4] = Scalar( -0.16527895766639 );
348  (*points_)[5] = Scalar( +0.16527895766639 );
349  (*points_)[6] = Scalar( +0.47792494981044 );
350  (*points_)[7] = Scalar( +0.73877386510551 );
351  (*points_)[8] = Scalar( +0.91953390816646 );
352  (*points_)[9] = ST::one();
353  (*weights_)[0] = Scalar( 0.02222222222222 );
354  (*weights_)[1] = Scalar( 0.13330599085107 );
355  (*weights_)[2] = Scalar( 0.22488934206313 );
356  (*weights_)[3] = Scalar( 0.29204268367968 );
357  (*weights_)[4] = Scalar( 0.32753976118390 );
358  (*weights_)[5] = Scalar( 0.32753976118390 );
359  (*weights_)[6] = Scalar( 0.29204268367968 );
360  (*weights_)[7] = Scalar( 0.22488934206313 );
361  (*weights_)[8] = Scalar( 0.13330599085107 );
362  (*weights_)[9] = Scalar( 0.02222222222222 );
363  }
364 }
365 
366 
367 template<class Scalar>
368 RCP<Thyra::VectorBase<Scalar> > eval_f_t(
369  const Thyra::ModelEvaluator<Scalar>& me,
370  Scalar t
371  ) {
372  typedef Teuchos::ScalarTraits<Scalar> ST;
373  typedef Thyra::ModelEvaluatorBase MEB;
374  MEB::InArgs<Scalar> inArgs = me.createInArgs();
375  inArgs.set_t(t);
376  MEB::OutArgs<Scalar> outArgs = me.createOutArgs();
377  RCP<Thyra::VectorBase<Scalar> > f_out = Thyra::createMember(me.get_f_space());
378  V_S(outArg(*f_out),ST::zero());
379  outArgs.set_f(f_out);
380  me.evalModel(inArgs,outArgs);
381  return f_out;
382 }
383 
384 template<class Scalar>
385 Scalar translateTimeRange(
386  Scalar t,
387  const TimeRange<Scalar>& sourceRange,
388  const TimeRange<Scalar>& destinationRange
389  ) {
390  Scalar r = destinationRange.length()/sourceRange.length();
391  return r*t+destinationRange.lower()-r*sourceRange.lower();
392 }
393 
394 template<class Scalar>
395 RCP<Thyra::VectorBase<Scalar> > computeArea(
396  const Thyra::ModelEvaluator<Scalar>& me,
397  const TimeRange<Scalar>& tr,
398  const GaussQuadrature1D<Scalar>& gq
399  ) {
400  typedef Teuchos::ScalarTraits<Scalar> ST;
401  RCP<Thyra::VectorBase<Scalar> > area = Thyra::createMember(me.get_x_space());
402  V_S(outArg(*area),ST::zero());
403  RCP<const TimeRange<Scalar> > sourceRange = gq.getRange();
404  RCP<const Array<Scalar> > sourcePts = gq.getPoints();
405  RCP<const Array<Scalar> > sourceWts = gq.getWeights();
406  Array<Scalar> destPts(*sourcePts);
407  for (unsigned int i=0 ; i<sourcePts->size() ; ++i) {
408  destPts[i] = translateTimeRange<Scalar>((*sourcePts)[i],*sourceRange,tr);
409  }
410  Scalar r = tr.length()/sourceRange->length();
411  for (unsigned int i=0 ; i<destPts.size() ; ++i) {
412  RCP<Thyra::VectorBase<Scalar> > tmpVec = eval_f_t<Scalar>(me,destPts[i]);
413  Vp_StV(outArg(*area),r*(*sourceWts)[i],*tmpVec);
414  }
415  return area;
416 }
417 
418 
419 } // namespace Rythmos
420 
421 #endif //Rythmos_QUADRATURE_BASE_H
Rythmos::GaussQuadrature1D
Specific implementation of 1D Gaussian based quadrature formulas.
Definition: Rythmos_QuadratureBase.hpp:44
Generated by   doxygen 1.8.5