Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "PolynomialQuadrature.h"
11 :
12 : #include "MooseError.h"
13 : #include "DataIO.h"
14 :
15 : // For computing legendre quadrature
16 : #include "libmesh/dense_matrix_impl.h"
17 :
18 : // For quickly computing polynomials
19 : #ifdef LIBMESH_HAVE_EXTERNAL_BOOST
20 : #pragma GCC diagnostic push
21 : #pragma GCC diagnostic ignored "-Wparentheses"
22 : #include <boost/math/special_functions/legendre.hpp>
23 : #include <boost/math/special_functions/hermite.hpp>
24 : #pragma GCC diagnostic pop
25 : #endif
26 :
27 : #include <cmath>
28 : #include <memory>
29 :
30 : namespace PolynomialQuadrature
31 : {
32 :
33 : std::unique_ptr<const Polynomial>
34 1266 : makePolynomial(const Distribution * dist)
35 : {
36 1266 : const Uniform * u_dist = dynamic_cast<const Uniform *>(dist);
37 1266 : if (u_dist)
38 3888 : return std::make_unique<const Legendre>(dist->getParam<Real>("lower_bound"),
39 : dist->getParam<Real>("upper_bound"));
40 :
41 294 : const Normal * n_dist = dynamic_cast<const Normal *>(dist);
42 294 : if (n_dist)
43 1176 : return std::make_unique<const Hermite>(dist->getParam<Real>("mean"),
44 : dist->getParam<Real>("standard_deviation"));
45 :
46 0 : ::mooseError("Polynomials for '", dist->type(), "' distributions have not been implemented.");
47 : return nullptr;
48 : }
49 :
50 : void
51 0 : Polynomial::store(std::ostream & /*stream*/, void * /*context*/) const
52 : {
53 : // Cannot be pure virtual because for dataLoad operations the base class must be constructed
54 0 : ::mooseError("Polynomial child class must override this method.");
55 : }
56 :
57 : void
58 0 : Polynomial::store(nlohmann::json & /*json*/) const
59 : {
60 : // Cannot be pure virtual because for dataLoad operations the base class must be constructed
61 0 : ::mooseError("Polynomial child class must override this method.");
62 : }
63 :
64 : Real
65 0 : Polynomial::compute(const unsigned int /*order*/, const Real /*x*/, const bool /*normalize*/) const
66 : {
67 0 : ::mooseError("Polynomial type has not been implemented.");
68 : return 0;
69 : }
70 :
71 : Real
72 0 : Polynomial::computeDerivative(const unsigned int /*order*/,
73 : const Real /*x*/,
74 : const unsigned int /*n*/) const
75 : {
76 0 : ::mooseError("Polynomial type has not been implemented.");
77 : return 0;
78 : }
79 :
80 : void
81 0 : Polynomial::clenshawQuadrature(const unsigned int /*order*/,
82 : std::vector<Real> & /*points*/,
83 : std::vector<Real> & /*weights*/) const
84 : {
85 0 : ::mooseError("Clenshaw quadrature has not been implemented for this polynomial type.");
86 : }
87 :
88 : Real
89 196125 : Polynomial::productIntegral(const std::vector<unsigned int> order) const
90 : {
91 196125 : const unsigned int nprod = order.size();
92 :
93 196125 : if (nprod == 1)
94 0 : return (order[0] == 0 ? 1.0 : 0.0);
95 196125 : else if (nprod == 2)
96 0 : return (order[0] == order[1] ? innerProduct(order[0]) : 0.0);
97 :
98 196125 : unsigned int poly_order = std::accumulate(order.begin(), order.end(), 0);
99 196125 : unsigned int quad_order = (poly_order % 2 == 0 ? poly_order : poly_order - 1) / 2;
100 :
101 : std::vector<Real> xq;
102 : std::vector<Real> wq;
103 196125 : gaussQuadrature(quad_order, xq, wq);
104 :
105 : Real val = 0.0;
106 815795 : for (unsigned int q = 0; q < xq.size(); ++q)
107 : {
108 619670 : Real prod = wq[q];
109 3042345 : for (unsigned int i = 0; i < nprod; ++i)
110 2422675 : prod *= compute(order[i], xq[q], false);
111 619670 : val += prod;
112 : }
113 :
114 196125 : return val / std::accumulate(wq.begin(), wq.end(), 0.0);
115 196125 : }
116 :
117 1142 : Legendre::Legendre(const Real lower_bound, const Real upper_bound)
118 1142 : : Polynomial(), _lower_bound(lower_bound), _upper_bound(upper_bound)
119 : {
120 1142 : }
121 :
122 : void
123 251 : Legendre::store(std::ostream & stream, void * context) const
124 : {
125 251 : std::string type = "Legendre";
126 251 : dataStore(stream, type, context);
127 251 : dataStore(stream, _lower_bound, context);
128 251 : dataStore(stream, _upper_bound, context);
129 251 : }
130 :
131 : void
132 120 : Legendre::store(nlohmann::json & json) const
133 : {
134 240 : json["type"] = "Legendre";
135 240 : json["lower_bound"] = _lower_bound;
136 120 : json["upper_bound"] = _upper_bound;
137 120 : }
138 :
139 : Real
140 4506685 : Legendre::compute(const unsigned int order, const Real x, const bool normalize) const
141 : {
142 4506685 : if ((x < _lower_bound) || (x > _upper_bound))
143 0 : ::mooseError("The requested polynomial point is outside of distribution bounds.");
144 :
145 4506685 : Real val = legendre(order, x, _lower_bound, _upper_bound);
146 4506685 : if (normalize)
147 2068880 : val /= innerProduct(order);
148 :
149 4506685 : return val;
150 : }
151 :
152 : Real
153 12800 : Legendre::computeDerivative(const unsigned int order, const Real x, const unsigned int m) const
154 : {
155 12800 : if ((x < _lower_bound) || (x > _upper_bound))
156 0 : ::mooseError("The requested polynomial point is outside of distribution bounds.");
157 :
158 12800 : Real xref = 2.0 / (_upper_bound - _lower_bound) * (x - (_upper_bound + _lower_bound) / 2.0);
159 : Real Jac = pow(2.0 / (_upper_bound - _lower_bound), m);
160 12800 : return Jac * computeDerivativeRef(order, xref, m);
161 : }
162 :
163 : Real
164 29504 : Legendre::computeDerivativeRef(const unsigned int order,
165 : const Real xref,
166 : const unsigned int m) const
167 : {
168 29504 : if (m == 0)
169 15312 : return legendre(order, xref);
170 14192 : else if (m > order)
171 : return 0.0;
172 13080 : else if (order == 1)
173 : return 1.0;
174 : else
175 8352 : return static_cast<Real>(order + m - 1) * computeDerivativeRef(order - 1, xref, m - 1) +
176 8352 : xref * computeDerivativeRef(order - 1, xref, m);
177 : }
178 :
179 : void
180 196250 : Legendre::gaussQuadrature(const unsigned int order,
181 : std::vector<Real> & points,
182 : std::vector<Real> & weights) const
183 : {
184 196250 : gauss_legendre(order, points, weights, _lower_bound, _upper_bound);
185 196250 : }
186 :
187 : void
188 4032 : Legendre::clenshawQuadrature(const unsigned int order,
189 : std::vector<Real> & points,
190 : std::vector<Real> & weights) const
191 : {
192 4032 : clenshaw_curtis(order, points, weights);
193 4032 : Real dx = (_upper_bound - _lower_bound) / 2.0;
194 4032 : Real xav = (_upper_bound + _lower_bound) / 2.0;
195 10848 : for (unsigned int i = 0; i < points.size(); ++i)
196 6816 : points[i] = points[i] * dx + xav;
197 4032 : }
198 :
199 : Real
200 2272586 : Legendre::innerProduct(const unsigned int order) const
201 : {
202 2272586 : return 1. / (2. * static_cast<Real>(order) + 1.);
203 : }
204 :
205 : Real
206 4521998 : legendre(const unsigned int order, const Real x, const Real lower_bound, const Real upper_bound)
207 : {
208 4521998 : Real xref = 2 / (upper_bound - lower_bound) * (x - (upper_bound + lower_bound) / 2);
209 : #ifdef LIBMESH_HAVE_EXTERNAL_BOOST
210 : // Using Legendre polynomials from boost library (if available)
211 : // https://www.boost.org/doc/libs/1_46_1/libs/math/doc/sf_and_dist/html/math_toolkit/special/sf_poly/legendre.html
212 : return boost::math::legendre_p(order, xref);
213 : #else
214 : // Using explicit expression of polynomial coefficients:
215 : // P_n(x) = 1/2^n\sum_{k=1}^{floor(n/2)} (-1)^k * nchoosek(n, k) * nchoosek(2n-2k, n) * x^(n-2k)
216 : // https://en.wikipedia.org/wiki/Legendre_polynomials
217 4521998 : if (order < 16)
218 : {
219 : Real val = 0;
220 11321819 : for (unsigned int k = 0; k <= (order % 2 == 0 ? order / 2 : (order - 1) / 2); ++k)
221 : {
222 : Real coeff =
223 6799821 : Real(Utility::binomial(order, k)) * Real(Utility::binomial(2 * order - 2 * k, order));
224 6799821 : Real sgn = (k % 2 == 0 ? 1.0 : -1.0);
225 6799821 : unsigned int ord = order - 2 * k;
226 6799821 : val += sgn * coeff * pow(xref, ord);
227 : }
228 4521998 : return val / pow(2.0, order);
229 : }
230 : else
231 : {
232 0 : Real ord = order;
233 0 : return ((2.0 * ord - 1.0) * xref * legendre(order - 1, xref) -
234 0 : (ord - 1.0) * legendre(order - 2, xref)) /
235 0 : ord;
236 : }
237 : #endif
238 : }
239 :
240 380 : Hermite::Hermite(const Real mu, const Real sig) : Polynomial(), _mu(mu), _sig(sig) {}
241 :
242 : Real
243 199248 : Hermite::innerProduct(const unsigned int order) const
244 : {
245 199248 : return (Real)Utility::factorial(order);
246 : }
247 :
248 : void
249 101 : Hermite::store(std::ostream & stream, void * context) const
250 : {
251 101 : std::string type = "Hermite";
252 101 : dataStore(stream, type, context);
253 101 : dataStore(stream, _mu, context);
254 101 : dataStore(stream, _sig, context);
255 101 : }
256 :
257 : void
258 0 : Hermite::store(nlohmann::json & json) const
259 : {
260 0 : json["type"] = "Hermite";
261 0 : json["mu"] = _mu;
262 0 : json["sig"] = _sig;
263 0 : }
264 :
265 : Real
266 80970 : Hermite::compute(const unsigned int order, const Real x, const bool normalize) const
267 : {
268 80970 : Real val = hermite(order, x, _mu, _sig);
269 80970 : if (normalize)
270 42120 : val /= innerProduct(order);
271 80970 : return val;
272 : }
273 :
274 : Real
275 12800 : Hermite::computeDerivative(const unsigned int order, const Real x, const unsigned int m) const
276 : {
277 12800 : if (m > order)
278 : return 0.0;
279 :
280 11688 : Real xref = (x - _mu) / _sig;
281 11688 : Real val = hermite(order - m, xref);
282 16416 : for (unsigned int i = 0; i < m; ++i)
283 4728 : val *= static_cast<Real>(order - i) / _sig;
284 :
285 : return val;
286 : }
287 :
288 : void
289 4287 : Hermite::gaussQuadrature(const unsigned int order,
290 : std::vector<Real> & points,
291 : std::vector<Real> & weights) const
292 : {
293 4287 : gauss_hermite(order, points, weights, _mu, _sig);
294 4287 : }
295 :
296 : Real
297 92659 : hermite(const unsigned int order, const Real x, const Real mu, const Real sig)
298 : {
299 92659 : Real xref = (x - mu) / sig;
300 : #ifdef LIBMESH_HAVE_EXTERNAL_BOOST
301 : // Using Hermite polynomials from boost library (if available)
302 : // https://www.boost.org/doc/libs/1_46_1/libs/math/doc/sf_and_dist/html/math_toolkit/special/sf_poly/hermite.html
303 : // Need to do some modification since boost does physicists hermite polynomials:
304 : // H_n^prob(x) = 2^(-n/2)H_n^phys(x / sqrt(2))
305 : xref /= M_SQRT2; // 1 / sqrt(2)
306 : Real val = boost::math::hermite(order, xref);
307 : val /= pow(M_SQRT2, order); // 2^(-order / 2)
308 : return val;
309 : #else
310 : // Using explicit expression of polynomial coefficients:
311 : // H_n(x) = n!\sum_{m=1}^{floor(n/2)} (-1)^m / (m!(n-2m)!) * x^(n-2m) / 2^m
312 : // https://en.wikipedia.org/wiki/Hermite_polynomials
313 92659 : if (order < 13)
314 : {
315 : Real val = 0;
316 289866 : for (unsigned int m = 0; m <= (order % 2 == 0 ? order / 2 : (order - 1) / 2); ++m)
317 : {
318 197207 : Real sgn = (m % 2 == 0 ? 1.0 : -1.0);
319 : Real coeff =
320 344555 : 1.0 / Real(Utility::factorial(m) * Utility::factorial(order - 2 * m)) / pow(2.0, m);
321 : unsigned int ord = order - 2 * m;
322 197207 : val += sgn * coeff * pow(xref, ord);
323 : }
324 92659 : return val * Utility::factorial(order);
325 : }
326 : else
327 0 : return xref * hermite(order - 1, xref) -
328 0 : (static_cast<Real>(order) - 1.0) * hermite(order - 2, xref);
329 : #endif
330 : }
331 :
332 : void
333 196251 : gauss_legendre(const unsigned int order,
334 : std::vector<Real> & points,
335 : std::vector<Real> & weights,
336 : const Real lower_bound,
337 : const Real upper_bound)
338 : {
339 196251 : unsigned int n = order + 1;
340 196251 : points.resize(n);
341 196251 : weights.resize(n);
342 :
343 196251 : DenseMatrix<Real> mat(n, n);
344 196251 : DenseVector<Real> lambda(n);
345 196251 : DenseVector<Real> lambdai(n);
346 196251 : DenseMatrix<Real> vec(n, n);
347 618265 : for (unsigned int i = 1; i < n; ++i)
348 : {
349 422014 : Real ri = i;
350 422014 : mat(i, i - 1) = ri / std::sqrt(((2. * ri - 1.) * (2. * ri + 1.)));
351 422014 : mat(i - 1, i) = mat(i, i - 1);
352 : }
353 : mat.evd_right(lambda, lambdai, vec);
354 :
355 196251 : Real dx = (upper_bound - lower_bound) / 2.0;
356 196251 : Real xav = (upper_bound + lower_bound) / 2.0;
357 814516 : for (unsigned int i = 0; i < n; ++i)
358 : {
359 618265 : points[i] = lambda(i) * dx + xav;
360 618265 : weights[i] = vec(0, i) * vec(0, i);
361 : }
362 392502 : }
363 :
364 : void
365 4288 : gauss_hermite(const unsigned int order,
366 : std::vector<Real> & points,
367 : std::vector<Real> & weights,
368 : const Real mu,
369 : const Real sig)
370 : {
371 : // Number of points needed
372 4288 : unsigned int n = order + 1;
373 4288 : points.resize(n);
374 4288 : weights.resize(n);
375 :
376 4288 : DenseMatrix<Real> mat(n, n);
377 4288 : DenseVector<Real> lambda(n);
378 4288 : DenseVector<Real> lambdai(n);
379 4288 : DenseMatrix<Real> vec(n, n);
380 8862 : for (unsigned int i = 1; i < n; ++i)
381 : {
382 4574 : mat(i, i - 1) = std::sqrt(static_cast<Real>(i));
383 4574 : mat(i - 1, i) = mat(i, i - 1);
384 : }
385 : mat.evd_right(lambda, lambdai, vec);
386 :
387 13150 : for (unsigned int i = 0; i < n; ++i)
388 : {
389 8862 : points[i] = mu + lambda(i) * sig;
390 8862 : weights[i] = vec(0, i) * vec(0, i);
391 : }
392 8576 : }
393 :
394 : void
395 4033 : clenshaw_curtis(const unsigned int order, std::vector<Real> & points, std::vector<Real> & weights)
396 : {
397 : // Number of points needed
398 4033 : unsigned int N = order + (order % 2);
399 4033 : points.resize(N + 1);
400 4033 : weights.resize(N + 1);
401 :
402 4033 : if (N == 0)
403 : {
404 2688 : points[0] = 0;
405 2688 : weights[0] = 1;
406 2688 : return;
407 : }
408 :
409 1345 : std::vector<Real> dk(N / 2 + 1);
410 4087 : for (unsigned int k = 0; k <= (N / 2); ++k)
411 5432 : dk[k] = ((k == 0 || k == (N / 2)) ? 1.0 : 2.0) / (1.0 - 4.0 * (Real)k * (Real)k);
412 :
413 4087 : for (unsigned int n = 0; n <= (N / 2); ++n)
414 : {
415 2742 : Real theta = (Real)n * M_PI / ((Real)N);
416 2742 : points[n] = -std::cos(theta);
417 8394 : for (unsigned int k = 0; k <= (N / 2); ++k)
418 : {
419 5652 : Real Dnk =
420 5652 : ((n == 0 || n == (N / 2)) ? 0.5 : 1.0) * std::cos((Real)k * theta * 2.0) / ((Real)N);
421 5652 : weights[n] += Dnk * dk[k];
422 : }
423 : }
424 :
425 2742 : for (unsigned int n = 0; n < (N / 2); ++n)
426 : {
427 1397 : points[N - n] = -points[n];
428 1397 : weights[N - n] = weights[n];
429 : }
430 1345 : weights[N / 2] *= 2.0;
431 1345 : }
432 :
433 138 : TensorGrid::TensorGrid(const std::vector<unsigned int> & npoints,
434 138 : std::vector<std::unique_ptr<const Polynomial>> & poly)
435 : {
436 138 : if (npoints.size() != poly.size())
437 0 : ::mooseError("List of number of 1D quadrature points must be same size as number of Polynomial "
438 : "objects.");
439 :
440 138 : std::vector<std::vector<Real>> qpoints_1D(poly.size());
441 138 : std::vector<std::vector<Real>> qweights_1D(poly.size());
442 518 : for (unsigned int d = 0; d < poly.size(); ++d)
443 380 : poly[d]->gaussQuadrature(npoints[d] - 1, qpoints_1D[d], qweights_1D[d]);
444 :
445 276 : _quad = std::make_unique<const StochasticTools::WeightedCartesianProduct<Real, Real>>(
446 : qpoints_1D, qweights_1D);
447 138 : }
448 :
449 8 : SmolyakGrid::SmolyakGrid(const unsigned int max_order,
450 8 : std::vector<std::unique_ptr<const Polynomial>> & poly)
451 8 : : _ndim(poly.size())
452 : {
453 :
454 : // Compute full tensor tuple
455 8 : std::vector<std::vector<unsigned int>> tuple_1d(_ndim);
456 56 : for (unsigned int d = 0; d < _ndim; ++d)
457 : {
458 48 : tuple_1d[d].resize(max_order);
459 240 : for (unsigned int i = 0; i < max_order; ++i)
460 192 : tuple_1d[d][i] = i;
461 : }
462 8 : StochasticTools::CartesianProduct<unsigned int> tensor_tuple(tuple_1d);
463 :
464 8 : _npoints.push_back(0);
465 8 : unsigned int maxq = max_order - 1;
466 8 : unsigned int minq = (max_order > _ndim ? max_order - _ndim : 0);
467 32776 : for (std::size_t p = 0; p < tensor_tuple.numRows(); ++p)
468 : {
469 32768 : std::vector<unsigned int> dorder = tensor_tuple.computeRow(p);
470 32768 : unsigned int q = std::accumulate(dorder.begin(), dorder.end(), 0);
471 32768 : if (q <= maxq && q >= minq)
472 : {
473 672 : int sgn = ((max_order - q - 1) % 2 == 0 ? 1 : -1);
474 1344 : _coeff.push_back(sgn * Utility::binomial(_ndim - 1, _ndim + q - max_order));
475 :
476 672 : std::vector<std::vector<Real>> qpoints_1D(_ndim);
477 672 : std::vector<std::vector<Real>> qweights_1D(_ndim);
478 4704 : for (unsigned int d = 0; d < poly.size(); ++d)
479 4032 : poly[d]->gaussQuadrature(dorder[d], qpoints_1D[d], qweights_1D[d]);
480 :
481 1344 : _quad.push_back(std::make_unique<const StochasticTools::WeightedCartesianProduct<Real, Real>>(
482 : qpoints_1D, qweights_1D));
483 672 : _npoints.push_back(_npoints.back() + _quad.back()->numRows());
484 672 : }
485 32768 : }
486 8 : }
487 :
488 : std::vector<Real>
489 0 : SmolyakGrid::quadraturePoint(const unsigned int n) const
490 : {
491 0 : unsigned int ind = gridIndex(n);
492 0 : return _quad[ind]->computeRow(n - _npoints[ind]);
493 : }
494 :
495 : Real
496 27300 : SmolyakGrid::quadraturePoint(const unsigned int n, const unsigned int dim) const
497 : {
498 27300 : unsigned int ind = gridIndex(n);
499 27300 : return _quad[ind]->computeValue(n - _npoints[ind], dim);
500 : }
501 :
502 : Real
503 2275 : SmolyakGrid::quadratureWeight(const unsigned int n) const
504 : {
505 2275 : unsigned int ind = gridIndex(n);
506 2275 : return static_cast<Real>(_coeff[ind]) * _quad[ind]->computeWeight(n - _npoints[ind]);
507 : }
508 :
509 : unsigned int
510 29575 : SmolyakGrid::gridIndex(const unsigned int n) const
511 : {
512 1355445 : for (unsigned int i = 0; i < _npoints.size() - 1; ++i)
513 1355445 : if (_npoints[i + 1] > n)
514 : return i;
515 :
516 0 : ::mooseError("Point index requested is greater than number of points.");
517 :
518 : return 0;
519 : }
520 :
521 8 : ClenshawCurtisGrid::ClenshawCurtisGrid(const unsigned int max_order,
522 8 : std::vector<std::unique_ptr<const Polynomial>> & poly)
523 8 : : _ndim(poly.size())
524 : {
525 : // Compute full tensor tuple
526 8 : std::vector<std::vector<unsigned int>> tuple_1d(_ndim);
527 56 : for (unsigned int d = 0; d < _ndim; ++d)
528 : {
529 48 : tuple_1d[d].resize(max_order);
530 240 : for (unsigned int i = 0; i < max_order; ++i)
531 192 : tuple_1d[d][i] = i;
532 : }
533 8 : StochasticTools::CartesianProduct<unsigned int> tensor_tuple(tuple_1d);
534 :
535 : // Curtis clenshaw has a lot of nested points,
536 : // so it behooves us to avoid duplicate points
537 : std::map<std::vector<Real>, Real> quad_map;
538 8 : unsigned int maxq = max_order - 1;
539 8 : unsigned int minq = (max_order > _ndim ? max_order - _ndim : 0);
540 32776 : for (std::size_t p = 0; p < tensor_tuple.numRows(); ++p)
541 : {
542 32768 : std::vector<unsigned int> dorder = tensor_tuple.computeRow(p);
543 32768 : unsigned int q = std::accumulate(dorder.begin(), dorder.end(), 0);
544 32768 : if (q <= maxq && q >= minq)
545 : {
546 672 : int sgn = ((max_order - q - 1) % 2 == 0 ? 1 : -1);
547 672 : Real coeff = static_cast<Real>(sgn * Utility::binomial(_ndim - 1, _ndim + q - max_order));
548 :
549 672 : std::vector<std::vector<Real>> qpoints_1D(_ndim);
550 672 : std::vector<std::vector<Real>> qweights_1D(_ndim);
551 4704 : for (unsigned int d = 0; d < poly.size(); ++d)
552 4032 : poly[d]->clenshawQuadrature(dorder[d], qpoints_1D[d], qweights_1D[d]);
553 :
554 672 : StochasticTools::WeightedCartesianProduct<Real, Real> quad(qpoints_1D, qweights_1D);
555 :
556 8768 : for (unsigned int i = 0; i < quad.numRows(); ++i)
557 8096 : quad_map[quad.computeRow(i)] += coeff * quad.computeWeight(i);
558 672 : }
559 32768 : }
560 :
561 8 : _quadrature.reserve(quad_map.size());
562 5656 : for (const auto & it : quad_map)
563 5648 : _quadrature.push_back(it);
564 8 : }
565 : } // namespace PolynomialQuadrature
566 :
567 : template <>
568 : void
569 352 : dataStore(std::ostream & stream,
570 : std::unique_ptr<const PolynomialQuadrature::Polynomial> & ptr,
571 : void * context)
572 : {
573 352 : ptr->store(stream, context);
574 352 : }
575 :
576 : template <>
577 : void
578 254 : dataLoad(std::istream & stream,
579 : std::unique_ptr<const PolynomialQuadrature::Polynomial> & ptr,
580 : void * context)
581 : {
582 : std::string poly_type;
583 254 : dataLoad(stream, poly_type, context);
584 254 : if (poly_type == "Legendre")
585 : {
586 : Real lower_bound, upper_bound;
587 169 : dataLoad(stream, lower_bound, context);
588 169 : dataLoad(stream, upper_bound, context);
589 169 : ptr = std::make_unique<const PolynomialQuadrature::Legendre>(lower_bound, upper_bound);
590 : }
591 85 : else if (poly_type == "Hermite")
592 : {
593 : Real mean, stddev;
594 85 : dataLoad(stream, mean, context);
595 85 : dataLoad(stream, stddev, context);
596 85 : ptr = std::make_unique<const PolynomialQuadrature::Hermite>(mean, stddev);
597 : }
598 : else
599 0 : ::mooseError("Unknown Polynomaial type: ", poly_type);
600 254 : }
|
