queso-0.56.0
1DQuadrature.C
Go to the documentation of this file.
1 //-----------------------------------------------------------------------bl-
2 //--------------------------------------------------------------------------
3 //
4 // QUESO - a library to support the Quantification of Uncertainty
5 // for Estimation, Simulation and Optimization
6 //
7 // Copyright (C) 2008-2015 The PECOS Development Team
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the Version 2.1 GNU Lesser General
11 // Public License as published by the Free Software Foundation.
12 //
13 // This library is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // Lesser General Public License for more details.
17 //
18 // You should have received a copy of the GNU Lesser General Public
19 // License along with this library; if not, write to the Free Software
20 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
21 // Boston, MA 02110-1301 USA
22 //
23 //-----------------------------------------------------------------------el-
24 
25 #include <queso/1DQuadrature.h>
26 
27 namespace QUESO {
28 
29 //*****************************************************
30 // Base 1D quadrature class
31 //*****************************************************
32 Base1DQuadrature::Base1DQuadrature(
33  double minDomainValue,
34  double maxDomainValue,
35  unsigned int order)
36  :
37  m_minDomainValue(minDomainValue),
38  m_maxDomainValue(maxDomainValue),
39  m_order (order),
40  m_positions (0),
41  m_weights (0)
42 {
43  queso_require_less_msg(m_minDomainValue, m_maxDomainValue, "min >= max");
44 
45  // UQ_UNAVAILABLE_RANK,
46  // "Base1DQuadrature::constructor()",
47  // "order = 0");
48 }
49 
50 Base1DQuadrature::~Base1DQuadrature()
51 {
52 }
53 
54 double
55 Base1DQuadrature::minDomainValue() const
56 {
57  return m_minDomainValue;
58 }
59 
60 double
61 Base1DQuadrature::maxDomainValue() const
62 {
63  return m_maxDomainValue;
64 }
65 
66 unsigned int
67 Base1DQuadrature::order() const
68 {
69  return m_order;
70 }
71 
72 const std::vector<double>&
73 Base1DQuadrature::positions() const
74 {
75  queso_require_not_equal_to_msg(m_positions.size(), 0, "size = 0");
76  return m_positions;
77 }
78 
79 const std::vector<double>&
80 Base1DQuadrature::weights() const
81 {
82  queso_require_not_equal_to_msg(m_weights.size(), 0, "size = 0");
83  return m_weights;
84 }
85 
86 //*****************************************************
87 // Generic 1D quadrature class
88 //*****************************************************
89 Generic1DQuadrature::Generic1DQuadrature(
90  double minDomainValue,
91  double maxDomainValue,
92  const std::vector<double>& positions,
93  const std::vector<double>& weights)
94  :
95  Base1DQuadrature(minDomainValue,maxDomainValue,positions.size()-1)
96 {
97  m_positions = positions;
98  m_weights = weights;
99 
100  queso_require_not_equal_to_msg(m_positions.size(), 0, "invalid positions");
101 
102  queso_require_equal_to_msg(m_positions.size(), m_weights.size(), "inconsistent positions and weight");
103 }
104 
105 Generic1DQuadrature::~Generic1DQuadrature()
106 {
107 }
108 
109 void
110 Generic1DQuadrature::dumbRoutine() const
111 {
112  return;
113 }
114 
115 //*****************************************************
116 // UniformLegendre 1D quadrature class
117 //*****************************************************
118 UniformLegendre1DQuadrature::UniformLegendre1DQuadrature(
119  double minDomainValue,
120  double maxDomainValue,
121  unsigned int order,
122  bool densityIsNormalized)
123  :
124  Base1DQuadrature(minDomainValue,maxDomainValue,order)
125 {
126  m_positions.resize(m_order+1,0.); // Yes, '+1'
127  m_weights.resize (m_order+1,0.); // Yes, '+1'
128 
129  // http://www.holoborodko.com/pavel/?page_id=679
130  switch (m_order) { // eep2011
131  case 1:
132  m_weights [0] = 1.;
133  m_weights [1] = 1.;
134 
135  m_positions[0] = -1./sqrt(3.);
136  m_positions[1] = 1./sqrt(3.);
137  break;
138 
139  case 2:
140  m_weights [0] = 5./9.;
141  m_weights [1] = 8./9.;
142  m_weights [2] = 5./9.;
143 
144  m_positions[0] = -sqrt(.6);
145  m_positions[1] = 0.;
146  m_positions[2] = sqrt(.6);
147  break;
148 
149  case 3:
150  m_weights [0] = 0.5 - sqrt(30.)/36.;
151  m_weights [1] = 0.5 + sqrt(30.)/36.;
152  m_weights [2] = 0.5 + sqrt(30.)/36.;
153  m_weights [3] = 0.5 - sqrt(30.)/36.;
154 
155  m_positions[0] = -sqrt(3.+2.*sqrt(1.2))/sqrt(7.);
156  m_positions[1] = -sqrt(3.-2.*sqrt(1.2))/sqrt(7.);
157  m_positions[2] = sqrt(3.-2.*sqrt(1.2))/sqrt(7.);
158  m_positions[3] = sqrt(3.+2.*sqrt(1.2))/sqrt(7.);
159  break;
160 
161  case 4:
162  m_weights [0] = (322.-13.*sqrt(70.))/900.; // 0.236926885
163  m_weights [1] = (322.+13.*sqrt(70.))/900.; // 0.478628670
164  m_weights [2] = 128./225.; // 0.568888889
165  m_weights [3] = (322.+13.*sqrt(70.))/900.;
166  m_weights [4] = (322.-13.*sqrt(70.))/900.;
167 
168  m_positions[0] = -sqrt(5.+2.*sqrt(10./7.))/3.;
169  m_positions[1] = -sqrt(5.-2.*sqrt(10./7.))/3.;
170  m_positions[2] = 0.;
171  m_positions[3] = sqrt(5.-2.*sqrt(10./7.))/3.;
172  m_positions[4] = sqrt(5.+2.*sqrt(10./7.))/3.;
173  break;
174 
175  case 5:
176  m_weights [0] = 0.1713244923791703450402961;
177  m_weights [1] = 0.3607615730481386075698335;
178  m_weights [2] = 0.4679139345726910473898703;
179  m_weights [3] = 0.4679139345726910473898703;
180  m_weights [4] = 0.3607615730481386075698335;
181  m_weights [5] = 0.1713244923791703450402961;
182 
183  m_positions[0] = -0.9324695142031520278123016;
184  m_positions[1] = -0.6612093864662645136613996;
185  m_positions[2] = -0.2386191860831969086305017;
186  m_positions[3] = 0.2386191860831969086305017;
187  m_positions[4] = 0.6612093864662645136613996;
188  m_positions[5] = 0.9324695142031520278123016;
189  break;
190 
191  case 6:
192  m_weights [0] = 0.1294849661688696932706114;
193  m_weights [1] = 0.2797053914892766679014678;
194  m_weights [2] = 0.3818300505051189449503698;
195  m_weights [3] = 0.4179591836734693877551020;
196  m_weights [4] = 0.3818300505051189449503698;
197  m_weights [5] = 0.2797053914892766679014678;
198  m_weights [6] = 0.1294849661688696932706114;
199 
200  m_positions[0] = -0.9491079123427585245261897;
201  m_positions[1] = -0.7415311855993944398638648;
202  m_positions[2] = -0.4058451513773971669066064;
203  m_positions[3] = 0.;
204  m_positions[4] = 0.4058451513773971669066064;
205  m_positions[5] = 0.7415311855993944398638648;
206  m_positions[6] = 0.9491079123427585245261897;
207  break;
208 
209  case 7:
210  m_weights [0] = 0.10122854;
211  m_weights [1] = 0.22238103;
212  m_weights [2] = 0.31370665;
213  m_weights [3] = 0.36268378;
214  m_weights [4] = 0.36268378;
215  m_weights [5] = 0.31370665;
216  m_weights [6] = 0.22238103;
217  m_weights [7] = 0.10122854;
218 
219  m_positions[0] = -0.96028986;
220  m_positions[1] = -0.79666648;
221  m_positions[2] = -0.52553241;
222  m_positions[3] = -0.18343464;
223  m_positions[4] = 0.18343464;
224  m_positions[5] = 0.52553241;
225  m_positions[6] = 0.79666648;
226  m_positions[7] = 0.96028986;
227  break;
228 
229  case 10:
230  m_weights [ 0] = 0.0556685671161736664827537;
231  m_weights [ 1] = 0.1255803694649046246346943;
232  m_weights [ 2] = 0.1862902109277342514260976;
233  m_weights [ 3] = 0.2331937645919904799185237;
234  m_weights [ 4] = 0.2628045445102466621806889;
235  m_weights [ 5] = 0.2729250867779006307144835;
236  m_weights [ 6] = 0.2628045445102466621806889;
237  m_weights [ 7] = 0.2331937645919904799185237;
238  m_weights [ 8] = 0.1862902109277342514260976;
239  m_weights [ 9] = 0.1255803694649046246346943;
240  m_weights [10] = 0.0556685671161736664827537;
241 
242  m_positions[ 0] = -0.9782286581460569928039380;
243  m_positions[ 1] = -0.8870625997680952990751578;
244  m_positions[ 2] = -0.7301520055740493240934163;
245  m_positions[ 3] = -0.5190961292068118159257257;
246  m_positions[ 4] = -0.2695431559523449723315320;
247  m_positions[ 5] = 0.;
248  m_positions[ 6] = 0.2695431559523449723315320;
249  m_positions[ 7] = 0.5190961292068118159257257;
250  m_positions[ 8] = 0.7301520055740493240934163;
251  m_positions[ 9] = 0.8870625997680952990751578;
252  m_positions[10] = 0.9782286581460569928039380;
253  break;
254 
255  case 11:
256  m_weights [ 0] = 0.0471753363865118271946160;
257  m_weights [ 1] = 0.1069393259953184309602547;
258  m_weights [ 2] = 0.1600783285433462263346525;
259  m_weights [ 3] = 0.2031674267230659217490645;
260  m_weights [ 4] = 0.2334925365383548087608499;
261  m_weights [ 5] = 0.2491470458134027850005624;
262  m_weights [ 6] = 0.2491470458134027850005624;
263  m_weights [ 7] = 0.2334925365383548087608499;
264  m_weights [ 8] = 0.2031674267230659217490645;
265  m_weights [ 9] = 0.1600783285433462263346525;
266  m_weights [10] = 0.1069393259953184309602547;
267  m_weights [11] = 0.0471753363865118271946160;
268 
269  m_positions[ 0] = -0.9815606342467192506905491;
270  m_positions[ 1] = -0.9041172563704748566784659;
271  m_positions[ 2] = -0.7699026741943046870368938;
272  m_positions[ 3] = -0.5873179542866174472967024;
273  m_positions[ 4] = -0.3678314989981801937526915;
274  m_positions[ 5] = -0.1252334085114689154724414;
275  m_positions[ 6] = 0.1252334085114689154724414;
276  m_positions[ 7] = 0.3678314989981801937526915;
277  m_positions[ 8] = 0.5873179542866174472967024;
278  m_positions[ 9] = 0.7699026741943046870368938;
279  m_positions[10] = 0.9041172563704748566784659;
280  m_positions[11] = 0.9815606342467192506905491;
281  break;
282 
283  case 12:
284  m_weights [ 0] = 0.0404840047653158795200216;
285  m_weights [ 1] = 0.0921214998377284479144218;
286  m_weights [ 2] = 0.1388735102197872384636018;
287  m_weights [ 3] = 0.1781459807619457382800467;
288  m_weights [ 4] = 0.2078160475368885023125232;
289  m_weights [ 5] = 0.2262831802628972384120902;
290  m_weights [ 6] = 0.2325515532308739101945895;
291  m_weights [ 7] = 0.2262831802628972384120902;
292  m_weights [ 8] = 0.2078160475368885023125232;
293  m_weights [ 9] = 0.1781459807619457382800467;
294  m_weights [10] = 0.1388735102197872384636018;
295  m_weights [11] = 0.0921214998377284479144218;
296  m_weights [12] = 0.0404840047653158795200216;
297 
298  m_positions[ 0] = -0.9841830547185881494728294;
299  m_positions[ 1] = -0.9175983992229779652065478;
300  m_positions[ 2] = -0.8015780907333099127942065;
301  m_positions[ 3] = -0.6423493394403402206439846;
302  m_positions[ 4] = -0.4484927510364468528779129;
303  m_positions[ 5] = -0.2304583159551347940655281;
304  m_positions[ 6] = 0.;
305  m_positions[ 7] = 0.2304583159551347940655281;
306  m_positions[ 8] = 0.4484927510364468528779129;
307  m_positions[ 9] = 0.6423493394403402206439846;
308  m_positions[10] = 0.8015780907333099127942065;
309  m_positions[11] = 0.9175983992229779652065478;
310  m_positions[12] = 0.9841830547185881494728294;
311  break;
312 #if 0
313  case 13:
314  m_weights [ 0] = ;
315  m_weights [ 1] = ;
316  m_weights [ 2] = ;
317  m_weights [ 3] = ;
318  m_weights [ 4] = ;
319  m_weights [ 5] = ;
320  m_weights [ 6] = ;
321  m_weights [ 7] = ;
322  m_weights [ 8] = ;
323  m_weights [ 9] = ;
324  m_weights [10] = ;
325  m_weights [11] = ;
326  m_weights [12] = ;
327  m_weights [13] = ;
328 
329  m_positions[ 0] = -;
330  m_positions[ 1] = -;
331  m_positions[ 2] = -;
332  m_positions[ 3] = -;
333  m_positions[ 4] = -;
334  m_positions[ 5] = -;
335  m_positions[ 6] = -;
336  m_positions[ 7] = ;
337  m_positions[ 8] = ;
338  m_positions[ 9] = ;
339  m_positions[10] = ;
340  m_positions[11] = ;
341  m_positions[12] = ;
342  m_positions[13] = ;
343  break;
344 #endif
345 #if 0
346  14 ±0.1080549487073436620662447 0.2152638534631577901958764
347  ±0.3191123689278897604356718 0.2051984637212956039659241
348  ±0.5152486363581540919652907 0.1855383974779378137417166
349  ±0.6872929048116854701480198 0.1572031671581935345696019
350  ±0.8272013150697649931897947 0.1215185706879031846894148
351  ±0.9284348836635735173363911 0.0801580871597602098056333
352  ±0.9862838086968123388415973 0.0351194603317518630318329
353 #endif
354 
355  case 16:
356  m_weights [ 0] = 0.0241483028685479319601100;
357  m_weights [ 1] = 0.0554595293739872011294402;
358  m_weights [ 2] = 0.0850361483171791808835354;
359  m_weights [ 3] = 0.1118838471934039710947884;
360  m_weights [ 4] = 0.1351363684685254732863200;
361  m_weights [ 5] = 0.1540457610768102880814316;
362  m_weights [ 6] = 0.1680041021564500445099707;
363  m_weights [ 7] = 0.1765627053669926463252710;
364  m_weights [ 8] = 0.1794464703562065254582656;
365  m_weights [ 9] = 0.1765627053669926463252710;
366  m_weights [10] = 0.1680041021564500445099707;
367  m_weights [11] = 0.1540457610768102880814316;
368  m_weights [12] = 0.1351363684685254732863200;
369  m_weights [13] = 0.1118838471934039710947884;
370  m_weights [14] = 0.0850361483171791808835354;
371  m_weights [15] = 0.0554595293739872011294402;
372  m_weights [16] = 0.0241483028685479319601100;
373 
374  m_positions[ 0] = -0.9905754753144173356754340;
375  m_positions[ 1] = -0.9506755217687677612227170;
376  m_positions[ 2] = -0.8802391537269859021229557;
377  m_positions[ 3] = -0.7815140038968014069252301;
378  m_positions[ 4] = -0.6576711592166907658503022;
379  m_positions[ 5] = -0.5126905370864769678862466;
380  m_positions[ 6] = -0.3512317634538763152971855;
381  m_positions[ 7] = -0.1784841814958478558506775;
382  m_positions[ 8] = 0.;
383  m_positions[ 9] = 0.1784841814958478558506775;
384  m_positions[10] = 0.3512317634538763152971855;
385  m_positions[11] = 0.5126905370864769678862466;
386  m_positions[12] = 0.6576711592166907658503022;
387  m_positions[13] = 0.7815140038968014069252301;
388  m_positions[14] = 0.8802391537269859021229557;
389  m_positions[15] = 0.9506755217687677612227170;
390  m_positions[16] = 0.9905754753144173356754340;
391  break;
392 
393  default:
394  std::cerr << "In UniformLegendre1DQuadrature::constructor()"
395  << ": m_order = " << m_order
396  << std::endl;
397  queso_error_msg("order not supported");
398  break;
399  }
400 
401  // Scale positions from the interval [-1, 1] to the interval [min,max]
402  for (unsigned int j = 0; j < m_positions.size(); ++j) {
403  m_positions[j] = .5*(m_maxDomainValue - m_minDomainValue)*m_positions[j] + .5*(m_maxDomainValue + m_minDomainValue);
404  if (densityIsNormalized) {
405  // Since \rho is "1/\Delta", we just multiply by ".5"
406  m_weights[j] *= .5;
407  }
408  else {
409  // Since \rho is "1", we multiply by ".5 * \Delta"
410  m_weights[j] *= .5*(m_maxDomainValue - m_minDomainValue);
411  }
412  }
413 }
414 
415 UniformLegendre1DQuadrature::~UniformLegendre1DQuadrature()
416 {
417 }
418 
419 void
420 UniformLegendre1DQuadrature::dumbRoutine() const
421 {
422  return;
423 }
424 
425 //*****************************************************
426 // GaussianHermite 1D quadrature class
427 //*****************************************************
428 GaussianHermite1DQuadrature::GaussianHermite1DQuadrature(
429  double mean,
430  double stddev,
431  unsigned int order)
432  :
433  Base1DQuadrature(-INFINITY,INFINITY,order),
434  m_mean (mean),
435  m_stddev(stddev)
436 {
437  // FIX ME: prepare code for mean != 0 and stddev != 1
438  m_positions.resize(m_order+1,0.); // Yes, '+1'
439  m_weights.resize (m_order+1,0.); // Yes, '+1'
440 
441  // http://www.efunda.com/math/num_integration/findgausshermite.cfm
442  switch (m_order) { // eep2011
443  case 1:
444  m_weights [0] = sqrt(M_PI)/2.;
445  m_weights [1] = sqrt(M_PI)/2.;
446 
447  m_positions[0] = -1./sqrt(2.);
448  m_positions[1] = 1./sqrt(2.);
449  break;
450 
451  case 2:
452  m_weights [0] = sqrt(M_PI)/6.;
453  m_weights [1] = 2.*sqrt(M_PI)/3.;
454  m_weights [2] = sqrt(M_PI)/6.;
455 
456  m_positions[0] = -sqrt(1.5);
457  m_positions[1] = 0.;
458  m_positions[2] = sqrt(1.5);
459  break;
460 
461  case 3:
462  m_weights [0] = sqrt(M_PI)/4./(3.-sqrt(6.));
463  m_weights [1] = sqrt(M_PI)/4./(3.+sqrt(6.));
464  m_weights [2] = sqrt(M_PI)/4./(3.+sqrt(6.));
465  m_weights [3] = sqrt(M_PI)/4./(3.-sqrt(6.));
466 
467  m_positions[0] = -sqrt(1.5+sqrt(1.5));
468  m_positions[1] = -sqrt(1.5-sqrt(1.5));
469  m_positions[2] = sqrt(1.5-sqrt(1.5));
470  m_positions[3] = sqrt(1.5+sqrt(1.5));
471  break;
472 
473  case 4:
474  m_weights [0] = 0.019953242049;
475  m_weights [1] = 0.393619323152;
476  m_weights [2] = 0.945308720483;
477  m_weights [3] = 0.393619323152;
478  m_weights [4] = 0.019953242059;
479 
480  m_positions[0] = -sqrt(2.5+sqrt(2.5));
481  m_positions[1] = -sqrt(2.5-sqrt(2.5));
482  m_positions[2] = 0.;
483  m_positions[3] = sqrt(2.5-sqrt(2.5));
484  m_positions[4] = sqrt(2.5+sqrt(2.5));
485  break;
486 
487  case 5:
488  m_weights [0] = 0.00453000990551;
489  m_weights [1] = 0.157067320323;
490  m_weights [2] = 0.724629595224;
491  m_weights [3] = 0.724629595224;
492  m_weights [4] = 0.157067320323;
493  m_weights [5] = 0.00453000990551;
494 
495  m_positions [0] = -2.35060497367;
496  m_positions [1] = -1.33584907401;
497  m_positions [2] = -0.436077411928;
498  m_positions [3] = 0.436077411928;
499  m_positions [4] = 1.33584907401;
500  m_positions [5] = 2.35060497367;
501  break;
502 
503  case 6:
504  m_weights [0] = 0.0009717812451;
505  m_weights [1] = 0.0545155828191;
506  m_weights [2] = 0.42560725261;
507  m_weights [3] = 0.810264617557;
508  m_weights [4] = 0.42560725261;
509  m_weights [5] = 0.0545155828191;
510  m_weights [6] = 0.0009717812451;
511 
512  m_positions [0] = -2.65196135684;
513  m_positions [1] = -1.67355162877;
514  m_positions [2] = -0.816287882859;
515  m_positions [3] = 0.;
516  m_positions [4] = 0.816287882859;
517  m_positions [5] = 1.67355162877;
518  m_positions [6] = 2.65196135684;
519  break;
520 
521  case 7:
522  m_weights [0] = 0.000199604072211;
523  m_weights [1] = 0.0170779830074;
524  m_weights [2] = 0.207802325815;
525  m_weights [3] = 0.661147012558;
526  m_weights [4] = 0.661147012558;
527  m_weights [5] = 0.207802325815;
528  m_weights [6] = 0.0170779830074;
529  m_weights [7] = 0.000199604072211;
530 
531  m_positions [0] = -2.93063742026;
532  m_positions [1] = -1.9816567567;
533  m_positions [2] = -1.15719371245;
534  m_positions [3] = -0.381186990207;
535  m_positions [4] = 0.381186990207;
536  m_positions [5] = 1.15719371245;
537  m_positions [6] = 1.9816567567;
538  m_positions [7] = 2.93063742026;
539  break;
540 
541  case 8:
542  m_weights [0] = 3.96069772633e-5;
543  m_weights [1] = 0.00494362427554;
544  m_weights [2] = 0.0884745273944;
545  m_weights [3] = 0.432651559003;
546  m_weights [4] = 0.720235215606;
547  m_weights [5] = 0.432651559003;
548  m_weights [6] = 0.0884745273944;
549  m_weights [7] = 0.00494362427554;
550  m_weights [8] = 3.96069772633e-5;
551 
552  m_positions [0] = -3.19099320178;
553  m_positions [1] = -2.26658058453;
554  m_positions [2] = -1.46855328922;
555  m_positions [3] = -0.723551018753;
556  m_positions [4] = 0.;
557  m_positions [5] = 0.723551018753;
558  m_positions [6] = 1.46855328922;
559  m_positions [7] = 2.26658058453;
560  m_positions [8] = 3.19099320178;
561  break;
562 
563  case 9:
564  m_weights [0] = 7.64043285523e-6;
565  m_weights [1] = 0.00134364574678;
566  m_weights [2] = 0.0338743944555;
567  m_weights [3] = 0.240138611082;
568  m_weights [4] = 0.610862633735;
569  m_weights [5] = 0.610862633735;
570  m_weights [6] = 0.240138611082;
571  m_weights [7] = 0.0338743944555;
572  m_weights [8] = 0.00134364574678;
573  m_weights [9] = 7.64043285523e-6;
574 
575  m_positions [0] = -3.43615911884;
576  m_positions [1] = -2.53273167423;
577  m_positions [2] = -1.7566836493;
578  m_positions [3] = -1.03661082979;
579  m_positions [4] = -0.342901327224;
580  m_positions [5] = 0.342901327224;
581  m_positions [6] = 1.03661082979;
582  m_positions [7] = 1.7566836493;
583  m_positions [8] = 2.53273167423;
584  m_positions [9] = 3.43615911884;
585  break;
586 
587  case 19:
588  m_weights [0] = 2.22939364554e-13;
589  m_weights [1] = 4.39934099226e-10;
590  m_weights [2] = 1.08606937077e-7;
591  m_weights [3] = 7.8025564785e-6;
592  m_weights [4] = 0.000228338636017;
593  m_weights [5] = 0.00324377334224;
594  m_weights [6] = 0.0248105208875;
595  m_weights [7] = 0.10901720602;
596  m_weights [8] = 0.286675505363;
597  m_weights [9] = 0.462243669601;
598  m_weights [10] = 0.462243669601;
599  m_weights [11] = 0.286675505363;
600  m_weights [12] = 0.10901720602;
601  m_weights [13] = 0.0248105208875;
602  m_weights [14] = 0.00324377334224;
603  m_weights [15] = 0.000228338636017;
604  m_weights [16] = 7.8025564785e-6;
605  m_weights [17] = 1.08606937077e-7;
606  m_weights [18] = 4.39934099226e-10;
607  m_weights [19] = 2.22939364554e-13;
608 
609  m_positions [0] = -5.38748089001;
610  m_positions [1] = -4.60368244955;
611  m_positions [2] = -3.94476404012;
612  m_positions [3] = -3.34785456738;
613  m_positions [4] = -2.78880605843;
614  m_positions [5] = -2.25497400209;
615  m_positions [6] = -1.73853771212;
616  m_positions [7] = -1.2340762154;
617  m_positions [8] = -0.737473728545;
618  m_positions [9] = -0.245340708301;
619  m_positions[10] = 0.245340708301;
620  m_positions[11] = 0.737473728545;
621  m_positions[12] = 1.2340762154;
622  m_positions[13] = 1.73853771212;
623  m_positions[14] = 2.25497400209;
624  m_positions[15] = 2.78880605843;
625  m_positions[16] = 3.34785456738;
626  m_positions[17] = 3.94476404012;
627  m_positions[18] = 4.60368244955;
628  m_positions[19] = 5.38748089001;
629  break;
630 
631  default:
632  std::cerr << "In GaussianHermite1DQuadrature::constructor()"
633  << ": m_order = " << m_order
634  << std::endl;
635  queso_error_msg("order not supported");
636  break;
637  }
638  for (unsigned int j = 0; j < (m_order+1); ++j) {
639  m_weights [j] *= sqrt(2.);
640  m_positions[j] *= sqrt(2.);
641  }
642 }
643 
644 GaussianHermite1DQuadrature::~GaussianHermite1DQuadrature()
645 {
646 }
647 
648 void
649 GaussianHermite1DQuadrature::dumbRoutine() const
650 {
651  return;
652 }
653 
654 //*****************************************************
655 // WignerInverseChebyshev1st 1D quadrature class
656 //*****************************************************
657 WignerInverseChebyshev1st1DQuadrature::WignerInverseChebyshev1st1DQuadrature(
658  double minDomainValue,
659  double maxDomainValue,
660  unsigned int order)
661  :
662  Base1DQuadrature(minDomainValue,maxDomainValue,order)
663 {
664  m_positions.resize(m_order+1,0.); // Yes, '+1'
665  m_weights.resize (m_order+1,0.); // Yes, '+1'
666 
667  // http://en.wikipedia.org/wiki/Chebyshev-Gauss_quadrature
668  switch (m_order) {
669  default:
670  queso_error_msg("order not supported");
671  break;
672  }
673 
674  // Scale positions from the interval [-1, 1] to the interval [min,max]
675  for (unsigned int j = 0; j < m_positions.size(); ++j) {
676  m_positions[j] = .5*(m_maxDomainValue - m_minDomainValue)*m_positions[j] + .5*(m_maxDomainValue + m_minDomainValue);
677  m_weights[j] *= .5*(m_maxDomainValue - m_minDomainValue);
678  }
679 }
680 
681 WignerInverseChebyshev1st1DQuadrature::~WignerInverseChebyshev1st1DQuadrature()
682 {
683 }
684 
685 void
686 WignerInverseChebyshev1st1DQuadrature::dumbRoutine() const
687 {
688  return;
689 }
690 
691 //*****************************************************
692 // WignerChebyshev2nd 1D quadrature class
693 //*****************************************************
694 WignerChebyshev2nd1DQuadrature::WignerChebyshev2nd1DQuadrature(
695  double minDomainValue,
696  double maxDomainValue,
697  unsigned int order)
698  :
699  Base1DQuadrature(minDomainValue,maxDomainValue,order)
700 {
701  m_positions.resize(m_order+1,0.); // Yes, '+1'
702  m_weights.resize (m_order+1,0.); // Yes, '+1'
703 
704  // http://en.wikipedia.org/wiki/Chebyshev-Gauss_quadrature
705  unsigned int n = m_order+1;
706  for (unsigned int i = 0; i < n; ++i) {
707  double angle = M_PI*((double)(i+1))/((double)(n+1));
708  double cosValue = cos(angle);
709  double sinValue = sin(angle);
710  m_positions[i] = cosValue;
711  m_weights[i] = ( M_PI/((double)(n+1)) )*sinValue*sinValue;
712  }
713 
714  // Scale positions from the interval [-1, 1] to the interval [min,max]
715  for (unsigned int j = 0; j < m_positions.size(); ++j) {
716  m_positions[j] = .5*(m_maxDomainValue - m_minDomainValue)*m_positions[j] + .5*(m_maxDomainValue + m_minDomainValue);
717  m_weights[j] *= .5*(m_maxDomainValue - m_minDomainValue);
718  }
719 }
720 
721 WignerChebyshev2nd1DQuadrature::~WignerChebyshev2nd1DQuadrature()
722 {
723 }
724 
725 void
726 WignerChebyshev2nd1DQuadrature::dumbRoutine() const
727 {
728  return;
729 }
730 
731 } // End namespace QUESO
QUESO::Base1DQuadrature::m_positions
std::vector< double > m_positions
Definition: 1DQuadrature.h:85
QUESO::GaussianHermite1DQuadrature::GaussianHermite1DQuadrature
GaussianHermite1DQuadrature(double mean, double stddev, unsigned int order)
Default constructor.
Definition: 1DQuadrature.C:428
QUESO::Generic1DQuadrature::Generic1DQuadrature
Generic1DQuadrature(double minDomainValue, double maxDomainValue, const std::vector< double > &positions, const std::vector< double > &weights)
Default constructor.
Definition: 1DQuadrature.C:89
QUESO::Base1DQuadrature::~Base1DQuadrature
virtual ~Base1DQuadrature()
Virtual destructor.
Definition: 1DQuadrature.C:50
QUESO::UniformLegendre1DQuadrature::UniformLegendre1DQuadrature
UniformLegendre1DQuadrature(double minDomainValue, double maxDomainValue, unsigned int order, bool densityIsNormalized)
Default constructor.
Definition: 1DQuadrature.C:118
QUESO::Base1DQuadrature::m_weights
std::vector< double > m_weights
Definition: 1DQuadrature.h:86
QUESO::Base1DQuadrature::maxDomainValue
double maxDomainValue() const
Returns the maximum value of the domain of the (one-dimensional) function.
Definition: 1DQuadrature.C:61
QUESO::Base1DQuadrature::Base1DQuadrature
Base1DQuadrature(double minDomainValue, double maxDomainValue, unsigned int order)
Default constructor.
Definition: 1DQuadrature.C:32
queso_require_equal_to_msg
#define queso_require_equal_to_msg(expr1, expr2, msg)
Definition: asserts.h:73
queso_require_less_msg
#define queso_require_less_msg(expr1, expr2, msg)
Definition: asserts.h:75
QUESO::Base1DQuadrature::weights
const std::vector< double > & weights() const
Array of the weights used in the numerical integration.
Definition: 1DQuadrature.C:80
QUESO::GaussianHermite1DQuadrature::dumbRoutine
void dumbRoutine() const
A bogus method.
Definition: 1DQuadrature.C:649
QUESO::UniformLegendre1DQuadrature::dumbRoutine
void dumbRoutine() const
A bogus method.
Definition: 1DQuadrature.C:420
QUESO::Base1DQuadrature
Base class for one-dimensional quadrature rules (numerical integration of functions).
Definition: 1DQuadrature.h:47
QUESO::WignerInverseChebyshev1st1DQuadrature::dumbRoutine
void dumbRoutine() const
A bogus method.
Definition: 1DQuadrature.C:686
queso_require_not_equal_to_msg
#define queso_require_not_equal_to_msg(expr1, expr2, msg)
Definition: asserts.h:74
QUESO::WignerChebyshev2nd1DQuadrature::WignerChebyshev2nd1DQuadrature
WignerChebyshev2nd1DQuadrature(double minDomainValue, double maxDomainValue, unsigned int order)
Default constructor.
Definition: 1DQuadrature.C:694
QUESO::WignerInverseChebyshev1st1DQuadrature::WignerInverseChebyshev1st1DQuadrature
WignerInverseChebyshev1st1DQuadrature(double minDomainValue, double maxDomainValue, unsigned int order)
TODO: Default constructor.
Definition: 1DQuadrature.C:657
QUESO::Base1DQuadrature::positions
const std::vector< double > & positions() const
Array of the positions for the numerical integration.
Definition: 1DQuadrature.C:73
QUESO::WignerChebyshev2nd1DQuadrature::dumbRoutine
void dumbRoutine() const
Bogus routine.
Definition: 1DQuadrature.C:726
queso_error_msg
#define queso_error_msg(msg)
Definition: asserts.h:47
QUESO::Generic1DQuadrature::dumbRoutine
void dumbRoutine() const
Bogus routine.
Definition: 1DQuadrature.C:110
QUESO::Base1DQuadrature::order
unsigned int order() const
Returns the order of the quadrature rule.
Definition: 1DQuadrature.C:67
QUESO::Base1DQuadrature::minDomainValue
double minDomainValue() const
Returns the minimum value of the domain of the (one-dimensional) function.
Definition: 1DQuadrature.C:55
QUESO::Generic1DQuadrature::~Generic1DQuadrature
~Generic1DQuadrature()
Destructor.
Definition: 1DQuadrature.C:105
Generated on Tue Nov 29 2016 10:53:11 for queso-0.56.0 by  doxygen 1.8.5