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