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 3018 : makePolynomial(const Distribution * dist)
35 : {
36 3018 : const Uniform * u_dist = dynamic_cast<const Uniform *>(dist);
37 3018 : if (u_dist)
38 9176 : return std::make_unique<const Legendre>(dist->getParam<Real>("lower_bound"),
39 : dist->getParam<Real>("upper_bound"));
40 :
41 724 : const Normal * n_dist = dynamic_cast<const Normal *>(dist);
42 724 : if (n_dist)
43 2896 : 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 392250 : Polynomial::productIntegral(const std::vector<unsigned int> order) const
90 : {
91 392250 : const unsigned int nprod = order.size();
92 :
93 392250 : if (nprod == 1)
94 0 : return (order[0] == 0 ? 1.0 : 0.0);
95 392250 : else if (nprod == 2)
96 0 : return (order[0] == order[1] ? innerProduct(order[0]) : 0.0);
97 :
98 392250 : unsigned int poly_order = std::accumulate(order.begin(), order.end(), 0);
99 392250 : 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 392250 : gaussQuadrature(quad_order, xq, wq);
104 :
105 : Real val = 0.0;
106 1631590 : for (unsigned int q = 0; q < xq.size(); ++q)
107 : {
108 1239340 : Real prod = wq[q];
109 6084690 : for (unsigned int i = 0; i < nprod; ++i)
110 4845350 : prod *= compute(order[i], xq[q], false);
111 1239340 : val += prod;
112 : }
113 :
114 392250 : return val / std::accumulate(wq.begin(), wq.end(), 0.0);
115 : }
116 :
117 2680 : Legendre::Legendre(const Real lower_bound, const Real upper_bound)
118 2680 : : Polynomial(), _lower_bound(lower_bound), _upper_bound(upper_bound)
119 : {
120 2680 : }
121 :
122 : void
123 321 : Legendre::store(std::ostream & stream, void * context) const
124 : {
125 321 : std::string type = "Legendre";
126 321 : dataStore(stream, type, context);
127 321 : dataStore(stream, _lower_bound, context);
128 321 : dataStore(stream, _upper_bound, context);
129 321 : }
130 :
131 : void
132 240 : Legendre::store(nlohmann::json & json) const
133 : {
134 480 : json["type"] = "Legendre";
135 480 : json["lower_bound"] = _lower_bound;
136 240 : json["upper_bound"] = _upper_bound;
137 240 : }
138 :
139 : Real
140 9010810 : Legendre::compute(const unsigned int order, const Real x, const bool normalize) const
141 : {
142 9010810 : if ((x < _lower_bound) || (x > _upper_bound))
143 0 : ::mooseError("The requested polynomial point is outside of distribution bounds.");
144 :
145 9010810 : Real val = legendre(order, x, _lower_bound, _upper_bound);
146 9010810 : if (normalize)
147 4137760 : val /= innerProduct(order);
148 :
149 9010810 : return val;
150 : }
151 :
152 : Real
153 26400 : Legendre::computeDerivative(const unsigned int order, const Real x, const unsigned int m) const
154 : {
155 26400 : if ((x < _lower_bound) || (x > _upper_bound))
156 0 : ::mooseError("The requested polynomial point is outside of distribution bounds.");
157 :
158 26400 : 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 26400 : return Jac * computeDerivativeRef(order, xref, m);
161 : }
162 :
163 : Real
164 61152 : Legendre::computeDerivativeRef(const unsigned int order,
165 : const Real xref,
166 : const unsigned int m) const
167 : {
168 61152 : if (m == 0)
169 31856 : return legendre(order, xref);
170 29296 : else if (m > order)
171 : return 0.0;
172 27040 : else if (order == 1)
173 : return 1.0;
174 : else
175 17376 : return static_cast<Real>(order + m - 1) * computeDerivativeRef(order - 1, xref, m - 1) +
176 17376 : xref * computeDerivativeRef(order - 1, xref, m);
177 : }
178 :
179 : void
180 395710 : Legendre::gaussQuadrature(const unsigned int order,
181 : std::vector<Real> & points,
182 : std::vector<Real> & weights) const
183 : {
184 395710 : gauss_legendre(order, points, weights, _lower_bound, _upper_bound);
185 395710 : }
186 :
187 : void
188 11088 : Legendre::clenshawQuadrature(const unsigned int order,
189 : std::vector<Real> & points,
190 : std::vector<Real> & weights) const
191 : {
192 11088 : clenshaw_curtis(order, points, weights);
193 11088 : Real dx = (_upper_bound - _lower_bound) / 2.0;
194 11088 : Real xav = (_upper_bound + _lower_bound) / 2.0;
195 29832 : for (unsigned int i = 0; i < points.size(); ++i)
196 18744 : points[i] = points[i] * dx + xav;
197 11088 : }
198 :
199 : Real
200 4571174 : Legendre::innerProduct(const unsigned int order) const
201 : {
202 4571174 : return 1. / (2. * static_cast<Real>(order) + 1.);
203 : }
204 :
205 : Real
206 9042667 : legendre(const unsigned int order, const Real x, const Real lower_bound, const Real upper_bound)
207 : {
208 9042667 : 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 9042667 : if (order < 16)
218 : {
219 : Real val = 0;
220 22643234 : for (unsigned int k = 0; k <= (order % 2 == 0 ? order / 2 : (order - 1) / 2); ++k)
221 : {
222 : Real coeff =
223 13600567 : Real(Utility::binomial(order, k)) * Real(Utility::binomial(2 * order - 2 * k, order));
224 13600567 : Real sgn = (k % 2 == 0 ? 1.0 : -1.0);
225 13600567 : unsigned int ord = order - 2 * k;
226 13600567 : val += sgn * coeff * pow(xref, ord);
227 : }
228 9042667 : 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 918 : Hermite::Hermite(const Real mu, const Real sig) : Polynomial(), _mu(mu), _sig(sig) {}
241 :
242 : Real
243 421510 : Hermite::innerProduct(const unsigned int order) const
244 : {
245 421510 : return (Real)Utility::factorial(order);
246 : }
247 :
248 : void
249 201 : Hermite::store(std::ostream & stream, void * context) const
250 : {
251 201 : std::string type = "Hermite";
252 201 : dataStore(stream, type, context);
253 201 : dataStore(stream, _mu, context);
254 201 : dataStore(stream, _sig, context);
255 201 : }
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 162180 : Hermite::compute(const unsigned int order, const Real x, const bool normalize) const
267 : {
268 162180 : Real val = hermite(order, x, _mu, _sig);
269 162180 : if (normalize)
270 84240 : val /= innerProduct(order);
271 162180 : return val;
272 : }
273 :
274 : Real
275 26400 : Hermite::computeDerivative(const unsigned int order, const Real x, const unsigned int m) const
276 : {
277 26400 : if (m > order)
278 : return 0.0;
279 :
280 24144 : Real xref = (x - _mu) / _sig;
281 24144 : Real val = hermite(order - m, xref);
282 33808 : for (unsigned int i = 0; i < m; ++i)
283 9664 : val *= static_cast<Real>(order - i) / _sig;
284 :
285 : return val;
286 : }
287 :
288 : void
289 8658 : Hermite::gaussQuadrature(const unsigned int order,
290 : std::vector<Real> & points,
291 : std::vector<Real> & weights) const
292 : {
293 8658 : gauss_hermite(order, points, weights, _mu, _sig);
294 8658 : }
295 :
296 : Real
297 186325 : hermite(const unsigned int order, const Real x, const Real mu, const Real sig)
298 : {
299 186325 : 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 186325 : if (order < 13)
314 : {
315 : Real val = 0;
316 583440 : for (unsigned int m = 0; m <= (order % 2 == 0 ? order / 2 : (order - 1) / 2); ++m)
317 : {
318 397115 : Real sgn = (m % 2 == 0 ? 1.0 : -1.0);
319 : Real coeff =
320 693985 : 1.0 / Real(Utility::factorial(m) * Utility::factorial(order - 2 * m)) / pow(2.0, m);
321 : unsigned int ord = order - 2 * m;
322 397115 : val += sgn * coeff * pow(xref, ord);
323 : }
324 186325 : 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 395711 : 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 395711 : unsigned int n = order + 1;
340 395711 : points.resize(n);
341 395711 : weights.resize(n);
342 :
343 395711 : DenseMatrix<Real> mat(n, n);
344 395711 : DenseVector<Real> lambda(n);
345 395711 : DenseVector<Real> lambdai(n);
346 395711 : DenseMatrix<Real> vec(n, n);
347 1241685 : for (unsigned int i = 1; i < n; ++i)
348 : {
349 845974 : Real ri = i;
350 845974 : mat(i, i - 1) = ri / std::sqrt(((2. * ri - 1.) * (2. * ri + 1.)));
351 845974 : mat(i - 1, i) = mat(i, i - 1);
352 : }
353 : mat.evd_right(lambda, lambdai, vec);
354 :
355 395711 : Real dx = (upper_bound - lower_bound) / 2.0;
356 395711 : Real xav = (upper_bound + lower_bound) / 2.0;
357 1637396 : for (unsigned int i = 0; i < n; ++i)
358 : {
359 1241685 : points[i] = lambda(i) * dx + xav;
360 1241685 : weights[i] = vec(0, i) * vec(0, i);
361 : }
362 791422 : }
363 :
364 : void
365 8659 : 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 8659 : unsigned int n = order + 1;
373 8659 : points.resize(n);
374 8659 : weights.resize(n);
375 :
376 8659 : DenseMatrix<Real> mat(n, n);
377 8659 : DenseVector<Real> lambda(n);
378 8659 : DenseVector<Real> lambdai(n);
379 8659 : DenseMatrix<Real> vec(n, n);
380 18056 : for (unsigned int i = 1; i < n; ++i)
381 : {
382 9397 : mat(i, i - 1) = std::sqrt(static_cast<Real>(i));
383 9397 : mat(i - 1, i) = mat(i, i - 1);
384 : }
385 : mat.evd_right(lambda, lambdai, vec);
386 :
387 26715 : for (unsigned int i = 0; i < n; ++i)
388 : {
389 18056 : points[i] = mu + lambda(i) * sig;
390 18056 : weights[i] = vec(0, i) * vec(0, i);
391 : }
392 17318 : }
393 :
394 : void
395 11089 : clenshaw_curtis(const unsigned int order, std::vector<Real> & points, std::vector<Real> & weights)
396 : {
397 : // Number of points needed
398 11089 : unsigned int N = order + (order % 2);
399 11089 : points.resize(N + 1);
400 11089 : weights.resize(N + 1);
401 :
402 11089 : if (N == 0)
403 : {
404 7392 : points[0] = 0;
405 7392 : weights[0] = 1;
406 7392 : return;
407 : }
408 :
409 3697 : std::vector<Real> dk(N / 2 + 1);
410 11227 : for (unsigned int k = 0; k <= (N / 2); ++k)
411 14924 : dk[k] = ((k == 0 || k == (N / 2)) ? 1.0 : 2.0) / (1.0 - 4.0 * (Real)k * (Real)k);
412 :
413 11227 : for (unsigned int n = 0; n <= (N / 2); ++n)
414 : {
415 7530 : Real theta = (Real)n * M_PI / ((Real)N);
416 7530 : points[n] = -std::cos(theta);
417 23010 : for (unsigned int k = 0; k <= (N / 2); ++k)
418 : {
419 15480 : Real Dnk =
420 15480 : ((n == 0 || n == (N / 2)) ? 0.5 : 1.0) * std::cos((Real)k * theta * 2.0) / ((Real)N);
421 15480 : weights[n] += Dnk * dk[k];
422 : }
423 : }
424 :
425 7530 : for (unsigned int n = 0; n < (N / 2); ++n)
426 : {
427 3833 : points[N - n] = -points[n];
428 3833 : weights[N - n] = weights[n];
429 : }
430 3697 : weights[N / 2] *= 2.0;
431 : }
432 :
433 372 : TensorGrid::TensorGrid(const std::vector<unsigned int> & npoints,
434 372 : std::vector<std::unique_ptr<const Polynomial>> & poly)
435 : {
436 372 : 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 372 : std::vector<std::vector<Real>> qpoints_1D(poly.size());
441 372 : std::vector<std::vector<Real>> qweights_1D(poly.size());
442 1402 : for (unsigned int d = 0; d < poly.size(); ++d)
443 1030 : poly[d]->gaussQuadrature(npoints[d] - 1, qpoints_1D[d], qweights_1D[d]);
444 :
445 744 : _quad = std::make_unique<const StochasticTools::WeightedCartesianProduct<Real, Real>>(
446 : qpoints_1D, qweights_1D);
447 372 : }
448 :
449 22 : SmolyakGrid::SmolyakGrid(const unsigned int max_order,
450 22 : std::vector<std::unique_ptr<const Polynomial>> & poly)
451 22 : : _ndim(poly.size())
452 : {
453 :
454 : // Compute full tensor tuple
455 22 : std::vector<std::vector<unsigned int>> tuple_1d(_ndim);
456 154 : for (unsigned int d = 0; d < _ndim; ++d)
457 : {
458 132 : tuple_1d[d].resize(max_order);
459 660 : for (unsigned int i = 0; i < max_order; ++i)
460 528 : tuple_1d[d][i] = i;
461 : }
462 22 : StochasticTools::CartesianProduct<unsigned int> tensor_tuple(tuple_1d);
463 :
464 22 : _npoints.push_back(0);
465 22 : unsigned int maxq = max_order - 1;
466 22 : unsigned int minq = (max_order > _ndim ? max_order - _ndim : 0);
467 90134 : for (std::size_t p = 0; p < tensor_tuple.numRows(); ++p)
468 : {
469 90112 : std::vector<unsigned int> dorder = tensor_tuple.computeRow(p);
470 90112 : unsigned int q = std::accumulate(dorder.begin(), dorder.end(), 0);
471 90112 : if (q <= maxq && q >= minq)
472 : {
473 1848 : int sgn = ((max_order - q - 1) % 2 == 0 ? 1 : -1);
474 3696 : _coeff.push_back(sgn * Utility::binomial(_ndim - 1, _ndim + q - max_order));
475 :
476 1848 : std::vector<std::vector<Real>> qpoints_1D(_ndim);
477 1848 : std::vector<std::vector<Real>> qweights_1D(_ndim);
478 12936 : for (unsigned int d = 0; d < poly.size(); ++d)
479 11088 : poly[d]->gaussQuadrature(dorder[d], qpoints_1D[d], qweights_1D[d]);
480 :
481 3696 : _quad.push_back(std::make_unique<const StochasticTools::WeightedCartesianProduct<Real, Real>>(
482 : qpoints_1D, qweights_1D));
483 1848 : _npoints.push_back(_npoints.back() + _quad.back()->numRows());
484 1848 : }
485 : }
486 22 : }
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 54600 : SmolyakGrid::quadraturePoint(const unsigned int n, const unsigned int dim) const
497 : {
498 54600 : unsigned int ind = gridIndex(n);
499 54600 : return _quad[ind]->computeValue(n - _npoints[ind], dim);
500 : }
501 :
502 : Real
503 4550 : SmolyakGrid::quadratureWeight(const unsigned int n) const
504 : {
505 4550 : unsigned int ind = gridIndex(n);
506 4550 : return static_cast<Real>(_coeff[ind]) * _quad[ind]->computeWeight(n - _npoints[ind]);
507 : }
508 :
509 : unsigned int
510 59150 : SmolyakGrid::gridIndex(const unsigned int n) const
511 : {
512 2710890 : for (unsigned int i = 0; i < _npoints.size() - 1; ++i)
513 2710890 : 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 22 : ClenshawCurtisGrid::ClenshawCurtisGrid(const unsigned int max_order,
522 22 : std::vector<std::unique_ptr<const Polynomial>> & poly)
523 22 : : _ndim(poly.size())
524 : {
525 : // Compute full tensor tuple
526 22 : std::vector<std::vector<unsigned int>> tuple_1d(_ndim);
527 154 : for (unsigned int d = 0; d < _ndim; ++d)
528 : {
529 132 : tuple_1d[d].resize(max_order);
530 660 : for (unsigned int i = 0; i < max_order; ++i)
531 528 : tuple_1d[d][i] = i;
532 : }
533 22 : 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 22 : unsigned int maxq = max_order - 1;
539 22 : unsigned int minq = (max_order > _ndim ? max_order - _ndim : 0);
540 90134 : for (std::size_t p = 0; p < tensor_tuple.numRows(); ++p)
541 : {
542 90112 : std::vector<unsigned int> dorder = tensor_tuple.computeRow(p);
543 90112 : unsigned int q = std::accumulate(dorder.begin(), dorder.end(), 0);
544 90112 : if (q <= maxq && q >= minq)
545 : {
546 1848 : int sgn = ((max_order - q - 1) % 2 == 0 ? 1 : -1);
547 1848 : Real coeff = static_cast<Real>(sgn * Utility::binomial(_ndim - 1, _ndim + q - max_order));
548 :
549 1848 : std::vector<std::vector<Real>> qpoints_1D(_ndim);
550 1848 : std::vector<std::vector<Real>> qweights_1D(_ndim);
551 12936 : for (unsigned int d = 0; d < poly.size(); ++d)
552 11088 : poly[d]->clenshawQuadrature(dorder[d], qpoints_1D[d], qweights_1D[d]);
553 :
554 1848 : StochasticTools::WeightedCartesianProduct<Real, Real> quad(qpoints_1D, qweights_1D);
555 :
556 24112 : for (unsigned int i = 0; i < quad.numRows(); ++i)
557 22264 : quad_map[quad.computeRow(i)] += coeff * quad.computeWeight(i);
558 1848 : }
559 : }
560 :
561 22 : _quadrature.reserve(quad_map.size());
562 15554 : for (const auto & it : quad_map)
563 15532 : _quadrature.push_back(it);
564 22 : }
565 : } // namespace PolynomialQuadrature
566 :
567 : template <>
568 : void
569 522 : dataStore(std::ostream & stream,
570 : std::unique_ptr<const PolynomialQuadrature::Polynomial> & ptr,
571 : void * context)
572 : {
573 522 : ptr->store(stream, context);
574 522 : }
575 :
576 : template <>
577 : void
578 578 : dataLoad(std::istream & stream,
579 : std::unique_ptr<const PolynomialQuadrature::Polynomial> & ptr,
580 : void * context)
581 : {
582 : std::string poly_type;
583 578 : dataLoad(stream, poly_type, context);
584 578 : if (poly_type == "Legendre")
585 : {
586 : Real lower_bound, upper_bound;
587 385 : dataLoad(stream, lower_bound, context);
588 385 : dataLoad(stream, upper_bound, context);
589 385 : ptr = std::make_unique<const PolynomialQuadrature::Legendre>(lower_bound, upper_bound);
590 : }
591 193 : else if (poly_type == "Hermite")
592 : {
593 : Real mean, stddev;
594 193 : dataLoad(stream, mean, context);
595 193 : dataLoad(stream, stddev, context);
596 193 : ptr = std::make_unique<const PolynomialQuadrature::Hermite>(mean, stddev);
597 : }
598 : else
599 0 : ::mooseError("Unknown Polynomaial type: ", poly_type);
600 578 : }
|
