Rythmos - Transient Integration for Differential Equations Version of the Day
Loading...
Searching...
No Matches
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
39namespace Rythmos {
40
43template<class Scalar>
44class GaussQuadrature1D : virtual public Teuchos::Describable
45{
46public:
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
55template<class Scalar>
56class 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
76template<class Scalar>
77GaussLegendreQuadrature1D<Scalar>::GaussLegendreQuadrature1D(int numNodes) {
78 fixQuadrature_(numNodes);
79}
80
81template<class Scalar>
82void 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
213template<class Scalar>
214class 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
234template<class Scalar>
235GaussLobattoQuadrature1D<Scalar>::GaussLobattoQuadrature1D(int numNodes) {
236 fixQuadrature_(numNodes);
237}
238
239template<class Scalar>
240void 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
367template<class Scalar>
368RCP<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
384template<class Scalar>
385Scalar 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
394template<class Scalar>
395RCP<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
Specific implementation of 1D Gaussian based quadrature formulas.