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 "Legendre.h"
11 : #include <functional>
12 :
13 : /**
14 : * The highest order of Legendre polynomials calculated directly instead of via the recurrence
15 : * relation
16 : */
17 : #define MAX_DIRECT_CALCULATION_LEGENDRE 12
18 :
19 0 : Legendre::Legendre() : SingleSeriesBasisInterface() {}
20 :
21 405 : Legendre::Legendre(const std::vector<MooseEnum> & domain,
22 : const std::vector<std::size_t> & order,
23 : MooseEnum expansion_type,
24 405 : MooseEnum generation_type)
25 405 : : SingleSeriesBasisInterface(domain, order, calculatedNumberOfTermsBasedOnOrder(order))
26 : {
27 405 : if (expansion_type == "orthonormal")
28 1 : _evaluateExpansionWrapper = [this]() { this->evaluateOrthonormal(); };
29 404 : else if (expansion_type == "sqrt_mu")
30 33516 : _evaluateExpansionWrapper = [this]() { this->evaluateSqrtMu(); };
31 360 : else if (expansion_type == "standard")
32 127032 : _evaluateExpansionWrapper = [this]() { this->evaluateStandard(); };
33 : else
34 0 : mooseError("The specified type of normalization for expansion does not exist");
35 :
36 405 : if (generation_type == "orthonormal")
37 218619 : _evaluateGenerationWrapper = [this]() { this->evaluateOrthonormal(); };
38 47 : else if (generation_type == "sqrt_mu")
39 62806 : _evaluateGenerationWrapper = [this]() { this->evaluateSqrtMu(); };
40 2 : else if (generation_type == "standard")
41 3 : _evaluateGenerationWrapper = [this]() { this->evaluateStandard(); };
42 : else
43 0 : mooseError("The specified type of normalization for generation does not exist");
44 405 : }
45 :
46 : std::size_t
47 405 : Legendre::calculatedNumberOfTermsBasedOnOrder(const std::vector<std::size_t> & order) const
48 : {
49 405 : return order[0] + 1;
50 : }
51 :
52 : void
53 380 : Legendre::checkPhysicalBounds(const std::vector<Real> & bounds) const
54 : {
55 : // Each Legendre series should have a min and max bound
56 380 : if (bounds.size() != 2)
57 0 : mooseError("Legend: Invalid number of bounds specified for single series!");
58 380 : }
59 :
60 : void
61 218261 : Legendre::evaluateOrthonormal()
62 : {
63 : std::size_t k;
64 : const Real & x = _standardized_location[0];
65 218261 : const Real x2 = x * x;
66 :
67 : /*
68 : * Use direct formula to efficiently evaluate the polynomials for n <= 12
69 : *
70 : * The performance benefit diminishes for higher n. It is expected that the cost of the direct
71 : * calculation nears that of the recurrence relation in the neighborhood of n == 15, although this
72 : * theory is untested due to only implementing the direct calculations up to n == 12.
73 : *
74 : * If you want to calculate the higher-order Legendre Coefficients and code them in then be my
75 : * guest.
76 : */
77 : // clang-format off
78 218261 : switch (_orders[0])
79 : {
80 13 : default:
81 : case MAX_DIRECT_CALCULATION_LEGENDRE: /* 12 */
82 13 : save(12, ((((((676039 * x2 - 1939938) * x2 + 2078505) * x2 - 1021020) * x2 + 225225) * x2 - 18018) * x2 + 231) / 1024
83 : * 12.5);
84 : libmesh_fallthrough();
85 :
86 13 : case 11:
87 13 : save(11, (((((88179 * x2 - 230945) * x2 + 218790) * x2 - 90090) * x2 + 15015) * x2 - 693) * x / 256
88 : * 11.5);
89 : libmesh_fallthrough();
90 :
91 13 : case 10:
92 13 : save(10, (((((46189 * x2 - 109395) * x2 + 90090) * x2 - 30030) * x2 + 3465) * x2 - 63) / 256
93 : * 10.5);
94 : libmesh_fallthrough();
95 :
96 13 : case 9:
97 13 : save(9, ((((12155 * x2 - 25740) * x2 + 18018) * x2 - 4620) * x2 + 315) * x / 128
98 : * 9.5);
99 : libmesh_fallthrough();
100 :
101 13 : case 8:
102 13 : save(8, ((((6435 * x2 - 12012) * x2 + 6930) * x2 - 1260) * x2 + 35) / 128
103 : * 8.5);
104 : libmesh_fallthrough();
105 :
106 13 : case 7:
107 13 : save(7, (((429 * x2 - 693) * x2 + 315) * x2 - 35) * x / 16
108 : * 7.5);
109 : libmesh_fallthrough();
110 :
111 13 : case 6:
112 13 : save(6, (((231 * x2 - 315) * x2 + 105) * x2 - 5) / 16
113 : * 6.5);
114 : libmesh_fallthrough();
115 :
116 21773 : case 5:
117 21773 : save(5, ((63 * x2 - 70) * x2 + 15) * x / 8
118 : * 5.5);
119 : libmesh_fallthrough();
120 :
121 32653 : case 4:
122 32653 : save(4, ((35 * x2 - 30) * x2 + 3) / 8
123 : * 4.5);
124 : libmesh_fallthrough();
125 :
126 218261 : case 3:
127 218261 : save(3, (5 * x2 - 3) * x / 2
128 : * 3.5);
129 : libmesh_fallthrough();
130 :
131 218261 : case 2:
132 218261 : save(2, (3 * x2 - 1) / 2
133 : * 2.5);
134 : libmesh_fallthrough();
135 :
136 218261 : case 1:
137 218261 : save(1, x
138 : * 1.5);
139 : libmesh_fallthrough();
140 :
141 218261 : case 0:
142 218261 : save(0, 1
143 : * 0.5);
144 : }
145 : // clang-format on
146 :
147 : /*
148 : * Evaluate any remaining polynomials.
149 : *
150 : * The original recurrence relation is:
151 : * (2 * k - 1) * x * L_(k-1) - (k - 1) * L_(k-2)
152 : * L_k = ---------------------------------------------
153 : * k
154 : *
155 : * However, for FXs we are using a the orthonormalized version of the polynomials, so each
156 : * polynomial L_k is multiplied by:
157 : * (2 * k + 1)
158 : * ----------- essentially: k + 0.5
159 : * 2
160 : * Reversing this in the previous polynomials and implementing for the current polynomial results
161 : * in the orthonormalized recurrence:
162 : * (2 * k + 1) / (k - 1) \
163 : * L_k = ----------- * | x * L_(k-1) - ----------- * L_(k-2) |
164 : * k \ (2 * k - 3) /
165 : *
166 : * The options are 1) to use this form, or 2) to not apply the orthonormalization at first, and
167 : * then loop through all the values in a second loop and then apply the orthonormalization.
168 : */
169 218336 : for (k = MAX_DIRECT_CALCULATION_LEGENDRE + 1; k <= _orders[0]; ++k)
170 75 : save(k, ((k + k + 1) / Real(k)) * (x * load(k - 1) - ((k - 1) / (k + k - 3.0)) * load(k - 2)));
171 218261 : }
172 :
173 : void
174 222906 : Legendre::evaluateStandard()
175 : {
176 : std::size_t k;
177 222906 : const Real x = _standardized_location[0];
178 222906 : const Real x2 = x * x;
179 :
180 : /*
181 : * Use direct formula to efficiently evaluate the polynomials for n <= 12
182 : *
183 : * The performance benefit diminishes for higher n. It is expected that the cost of the direct
184 : * calculation nears that of the recurrence relation in the neighborhood of n == 15, although this
185 : * theory is untested due to only implementing the direct calculations up to n == 12.
186 : *
187 : * If you want to calculate the higher-order Legendre Coefficients and
188 : * code them in then be my guest.
189 : */
190 : // clang-format off
191 222906 : switch (_orders[0])
192 : {
193 14 : default:
194 : case MAX_DIRECT_CALCULATION_LEGENDRE: /* 12 */
195 14 : save(12, ((((((676039 * x2 - 1939938) * x2 + 2078505) * x2 - 1021020) * x2 + 225225) * x2 - 18018) * x2 + 231) / 1024);
196 : libmesh_fallthrough();
197 :
198 14 : case 11:
199 14 : save(11, (((((88179 * x2 - 230945) * x2 + 218790) * x2 - 90090) * x2 + 15015) * x2 - 693) * x / 256);
200 : libmesh_fallthrough();
201 :
202 14 : case 10:
203 14 : save(10, (((((46189 * x2 - 109395) * x2 + 90090) * x2 - 30030) * x2 + 3465) * x2 - 63) / 256);
204 : libmesh_fallthrough();
205 :
206 14 : case 9:
207 14 : save(9, ((((12155 * x2 - 25740) * x2 + 18018) * x2 - 4620) * x2 + 315) * x / 128);
208 : libmesh_fallthrough();
209 :
210 14 : case 8:
211 14 : save(8, ((((6435 * x2 - 12012) * x2 + 6930) * x2 - 1260) * x2 + 35) / 128);
212 : libmesh_fallthrough();
213 :
214 14 : case 7:
215 14 : save(7, (((429 * x2 - 693) * x2 + 315) * x2 - 35) * x / 16);
216 : libmesh_fallthrough();
217 :
218 14 : case 6:
219 14 : save(6, (((231 * x2 - 315) * x2 + 105) * x2 - 5) / 16);
220 : libmesh_fallthrough();
221 :
222 23215 : case 5:
223 23215 : save(5, ((63 * x2 - 70) * x2 + 15) * x / 8);
224 : libmesh_fallthrough();
225 :
226 29685 : case 4:
227 29685 : save(4, ((35 * x2 - 30) * x2 + 3) / 8);
228 : libmesh_fallthrough();
229 :
230 222906 : case 3:
231 222906 : save(3, (5 * x2 - 3) * x / 2);
232 : libmesh_fallthrough();
233 :
234 222906 : case 2:
235 222906 : save(2, (3 * x2 - 1) / 2);
236 : libmesh_fallthrough();
237 :
238 222906 : case 1:
239 222906 : save(1, x);
240 : libmesh_fallthrough();
241 :
242 222906 : case 0:
243 222906 : save(0, 1);
244 : }
245 : // clang-format on
246 :
247 : /*
248 : * Evaluate any remaining polynomials.
249 : * The recurrence relation is:
250 : * (2 * k - 1) * x * L_(k-1) - (k - 1) * L_(k-2)
251 : * L_k = ---------------------------------------------
252 : * k
253 : */
254 222984 : for (k = MAX_DIRECT_CALCULATION_LEGENDRE + 1; k <= _orders[0]; ++k)
255 78 : save(k, (((2 * k - 1) * x * load(k - 1)) - ((k - 1) * load(k - 2))) / Real(k));
256 222906 : }
257 :
258 : void
259 96233 : Legendre::evaluateSqrtMu()
260 : {
261 96233 : evaluateStandard();
262 481177 : for (size_t i = 0; i < getNumberOfTerms(); ++i)
263 384944 : save(i, load(i) * std::sqrt(i + 0.5));
264 96233 : }
265 :
266 : const std::vector<Real> &
267 0 : Legendre::getStandardizedFunctionLimits() const
268 : {
269 : // Lazily instantiate the function limits array
270 0 : static const std::vector<Real> standardizedFunctionLimits = {-1, 1};
271 :
272 0 : return standardizedFunctionLimits;
273 : }
274 :
275 : Real
276 320 : Legendre::getStandardizedFunctionVolume() const
277 : {
278 320 : return 2.0; // Span of [-1, 1]
279 : }
280 :
281 : std::vector<Real>
282 441052 : Legendre::getStandardizedLocation(const std::vector<Real> & location) const
283 : {
284 441052 : const Real difference = location[0] - _physical_bounds[0];
285 441052 : const Real span = _physical_bounds[1] - _physical_bounds[0];
286 :
287 : // Convert to [0, 1] (assuming that location[0] is within _physical_bounds)
288 441052 : const Real ratio = difference / span;
289 :
290 : // Legendre space is [-1, 1]
291 441052 : return {ratio * 2 - 1};
292 : }
293 :
294 : bool
295 307070 : Legendre::isInPhysicalBounds(const Point & point) const
296 : {
297 307070 : std::vector<Real> location = extractLocationFromPoint(point);
298 :
299 307070 : if (location[0] < _physical_bounds[0] || _physical_bounds[1] < location[0])
300 : return false;
301 : else
302 300632 : return true;
303 307070 : }
|
