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