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