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 "MooseUtils.h"
11 : #include "Zernike.h"
12 : #include <functional>
13 :
14 : /**
15 : * The higherst order of Zernike polynomials calculated directly instead of via the recurrence
16 : * relation
17 : */
18 : #define MAX_DIRECT_CALCULATION_ZERNIKE 10
19 :
20 0 : Zernike::Zernike() : SingleSeriesBasisInterface() {}
21 :
22 11 : Zernike::Zernike(const std::vector<MooseEnum> & domain,
23 : const std::vector<std::size_t> & order,
24 : MooseEnum expansion_type,
25 11 : MooseEnum generation_type)
26 11 : : SingleSeriesBasisInterface(domain, order, calculatedNumberOfTermsBasedOnOrder(order))
27 : {
28 11 : if (expansion_type == "orthonormal")
29 1 : _evaluateExpansionWrapper = [this]() { this->evaluateOrthonormal(); };
30 10 : else if (expansion_type == "sqrt_mu")
31 0 : _evaluateExpansionWrapper = [this]() { this->evaluateSqrtMu(); };
32 10 : else if (expansion_type == "standard")
33 14 : _evaluateExpansionWrapper = [this]() { this->evaluateStandard(); };
34 : else
35 0 : mooseError("The specified type of normalization for expansion does not exist");
36 :
37 11 : if (generation_type == "orthonormal")
38 14 : _evaluateGenerationWrapper = [this]() { this->evaluateOrthonormal(); };
39 3 : else if (generation_type == "sqrt_mu")
40 2 : _evaluateGenerationWrapper = [this]() { this->evaluateSqrtMu(); };
41 2 : else if (generation_type == "standard")
42 3 : _evaluateGenerationWrapper = [this]() { this->evaluateStandard(); };
43 : else
44 0 : mooseError("The specified type of normalization for generation does not exist");
45 11 : }
46 :
47 : std::size_t
48 11 : Zernike::calculatedNumberOfTermsBasedOnOrder(const std::vector<std::size_t> & order) const
49 : {
50 11 : return ((order[0] + 1) * (order[0] + 2)) / 2;
51 : }
52 :
53 : void
54 0 : Zernike::checkPhysicalBounds(const std::vector<Real> & bounds) const
55 : {
56 : /*
57 : * Each single series is assumed to be a function of a single coordinate, which should only have
58 : * two bounds.
59 : */
60 0 : if (bounds.size() != 3)
61 0 : mooseError("Zernike: Invalid number of bounds specified for single series!");
62 0 : }
63 :
64 : // clang-format off
65 : void
66 6 : Zernike::evaluateOrthonormal()
67 : {
68 : std::size_t n;
69 : long j, q;
70 : Real H1, H2, H3;
71 : const Real & rho = _standardized_location[0];
72 6 : const Real rho2 = rho * rho;
73 6 : const Real rho4 = rho2 * rho2;
74 :
75 6 : if (MooseUtils::absoluteFuzzyEqual(rho, 0.0))
76 : {
77 0 : for (n = 0; n <= _orders[0]; n += 2)
78 : {
79 0 : j = simpleDoubleToSingle(n, 0);
80 :
81 0 : if ((n / 2) % 2 != 0)
82 0 : save(j, -1 * (n + 1) / M_PI);
83 : else
84 0 : save(j, 1 * (n + 1) / M_PI);
85 : }
86 :
87 : return;
88 : }
89 :
90 6 : switch (_orders[0])
91 : {
92 3 : default:
93 : case MAX_DIRECT_CALCULATION_ZERNIKE: /* 10 */
94 3 : save(65, rho4 * rho4 * rho2
95 3 : * 22 / M_PI);
96 3 : save(64, (10 * rho2 - 9) * rho4 * rho4
97 3 : * 22 / M_PI);
98 3 : save(63, ((45 * rho2 - 72) * rho2 + 28) * rho4 * rho2
99 3 : * 22 / M_PI);
100 3 : save(62, (((120 * rho2 - 252) * rho2 + 168) * rho2- 35) * rho4
101 3 : * 22 / M_PI);
102 3 : save(61, ((((210 * rho2 - 504) * rho2 + 420) * rho2 - 140) * rho2 + 15) * rho2
103 3 : * 22 / M_PI);
104 3 : save(60, (((((252 * rho2 - 630) * rho2 + 560) * rho2 - 210) * rho2 + 30) * rho2 - 1)
105 3 : * 11 / M_PI);
106 : libmesh_fallthrough();
107 :
108 3 : case 9:
109 3 : save(54, rho4 * rho4 * rho
110 3 : * 20 / M_PI);
111 3 : save(53, (9 * rho2 - 8) * rho4 * rho2 * rho
112 3 : * 20 / M_PI);
113 3 : save(52, ((36 * rho2 - 56) * rho2 + 21) * rho4 * rho
114 3 : * 20 / M_PI);
115 3 : save(51, (((84 * rho2 - 168) * rho2 + 105) * rho2 - 20) * rho2 * rho
116 3 : * 20 / M_PI);
117 3 : save(50, ((((126 * rho2 - 280) * rho2 + 210) * rho2 - 60) * rho2 + 5) * rho
118 3 : * 20 / M_PI);
119 : libmesh_fallthrough();
120 :
121 3 : case 8:
122 3 : save(44, rho4 * rho4
123 3 : * 18 / M_PI);
124 3 : save(43, (8 * rho2 - 7) * rho4 * rho2
125 3 : * 18 / M_PI);
126 3 : save(42, ((28 * rho2 - 42) * rho2 + 15) * rho4
127 3 : * 18 / M_PI);
128 3 : save(41, (((56 * rho2 - 105) * rho2 + 60) * rho2 - 10) * rho2
129 3 : * 18 / M_PI);
130 3 : save(40, ((((70 * rho2 - 140) * rho2 + 90) * rho2 - 20) * rho2 + 1)
131 3 : * 9 / M_PI);
132 : libmesh_fallthrough();
133 :
134 3 : case 7:
135 3 : save(35, rho4 * rho2 * rho
136 3 : * 16 / M_PI);
137 3 : save(34, (7 * rho2 - 6) * rho4 * rho
138 3 : * 16 / M_PI);
139 3 : save(33, ((21 * rho2 - 30) * rho2 + 10) * rho2 * rho
140 3 : * 16 / M_PI);
141 3 : save(32, (((35 * rho2 - 60) * rho2 + 30) * rho2 - 4) * rho
142 3 : * 16 / M_PI);
143 : libmesh_fallthrough();
144 :
145 3 : case 6:
146 3 : save(27, rho4 * rho2
147 3 : * 14 / M_PI);
148 3 : save(26, (6 * rho2 - 5) * rho4
149 3 : * 14 / M_PI);
150 3 : save(25, ((15 * rho2 - 20) * rho2 + 6) * rho2
151 3 : * 14 / M_PI);
152 3 : save(24, (((20 * rho2 - 30) * rho2 + 12) * rho2 - 1)
153 3 : * 7 / M_PI);
154 : libmesh_fallthrough();
155 :
156 3 : case 5:
157 3 : save(20, rho4 * rho
158 3 : * 12 / M_PI);
159 3 : save(19, (5 * rho2 - 4) * rho2 * rho
160 3 : * 12 / M_PI);
161 3 : save(18, ((10 * rho2 - 12) * rho2 + 3) * rho
162 3 : * 12 / M_PI);
163 : libmesh_fallthrough();
164 :
165 6 : case 4:
166 6 : save(14, rho4
167 6 : * 10 / M_PI);
168 6 : save(13, (4 * rho2 - 3) * rho2
169 6 : * 10 / M_PI);
170 6 : save(12, ((6 * rho2 - 6) * rho2 + 1)
171 6 : * 5 / M_PI);
172 : libmesh_fallthrough();
173 :
174 6 : case 3:
175 6 : save(9, rho2 * rho
176 6 : * 8 / M_PI);
177 6 : save(8, (3 * rho2 - 2) * rho
178 6 : * 8 / M_PI);
179 : libmesh_fallthrough();
180 :
181 6 : case 2:
182 6 : save(5, rho2
183 6 : * 6 / M_PI);
184 6 : save(4, (2 * rho2 - 1)
185 6 : * 3 / M_PI);
186 : libmesh_fallthrough();
187 :
188 6 : case 1:
189 6 : save(2, rho
190 6 : * 4 / M_PI);
191 : libmesh_fallthrough();
192 :
193 6 : case 0:
194 6 : save(0, 1
195 : * 1 / M_PI);
196 : }
197 :
198 27 : for (n = MAX_DIRECT_CALCULATION_ZERNIKE + 1; n <= _orders[0]; ++n)
199 : {
200 21 : j = simpleDoubleToSingle(n, n);
201 21 : save(j, pow(rho, n)
202 21 : * (n + n + 2) / M_PI);
203 :
204 21 : j--;
205 21 : save(j, n * load(j + 1) - (n + 1) * load(j - (n + n)));
206 :
207 141 : for (q = n; q >= 4; q -= 2)
208 : {
209 120 : H3 = (-4 * (q - 2) * (q - 3)) / ((n + q - 2) * (n - q + 4.0));
210 120 : H2 = (H3 * (n + q) * (n - q + 2)) / (4.0 * (q - 1)) + (q - 2);
211 120 : H1 = q * (q - 1) / 2 - q * H2 + (H3 * (n + q + 2) * (n - q)) / 8.0;
212 120 : j--;
213 120 : if (q == 4)
214 9 : save(j, (H1 * load(j + 2) + (H2 + H3 / rho2) * load(j + 1))
215 : * 0.5);
216 : else
217 111 : save(j, H1 * load(j + 2) + (H2 + H3 / rho2) * load(j + 1));
218 : }
219 : }
220 :
221 6 : fillOutNegativeRankAndApplyAzimuthalComponent();
222 : }
223 : // clang-format on
224 :
225 : void
226 1 : Zernike::evaluateSqrtMu()
227 : {
228 1 : evaluateStandard();
229 1 : save(0, load(0) / std::sqrt(M_PI));
230 : size_t j = 1;
231 5 : for (size_t n = 1; n < _orders[0] + 1; ++n)
232 : {
233 18 : for (size_t m = 0; m < n + 1; ++m)
234 : {
235 14 : if (m != 0 && n / m == 2 && n % m == 0)
236 2 : save(j, load(j) * std::sqrt((n + 1) / M_PI));
237 : else
238 12 : save(j, load(j) * std::sqrt((2 * n + 2) / M_PI));
239 14 : ++j;
240 : }
241 : }
242 1 : }
243 :
244 : void
245 6 : Zernike::evaluateStandard()
246 : {
247 : std::size_t n;
248 : long j, q;
249 : Real H1, H2, H3;
250 : const Real & rho = _standardized_location[0];
251 6 : const Real rho2 = rho * rho;
252 6 : const Real rho4 = rho2 * rho2;
253 :
254 6 : if (MooseUtils::absoluteFuzzyLessEqual(rho, 0))
255 : {
256 10 : for (n = 0; n <= _orders[0]; n += 2)
257 : {
258 9 : j = simpleDoubleToSingle(n, 0);
259 :
260 9 : if ((n / 2) % 2 != 0)
261 4 : save(j, -1);
262 : else
263 5 : save(j, 1);
264 : }
265 :
266 : return;
267 : }
268 :
269 5 : switch (_orders[0])
270 : {
271 2 : default:
272 : case MAX_DIRECT_CALCULATION_ZERNIKE: /* 10 */
273 2 : save(65, rho4 * rho4 * rho2);
274 2 : save(64, (10 * rho2 - 9) * rho4 * rho4);
275 2 : save(63, ((45 * rho2 - 72) * rho2 + 28) * rho4 * rho2);
276 2 : save(62, (((120 * rho2 - 252) * rho2 + 168) * rho2 - 35) * rho4);
277 2 : save(61, ((((210 * rho2 - 504) * rho2 + 420) * rho2 - 140) * rho2 + 15) * rho2);
278 2 : save(60, ((((252 * rho2 - 630) * rho2 + 560) * rho2 - 210) * rho2 + 30) * rho2 - 1);
279 : libmesh_fallthrough();
280 :
281 2 : case 9:
282 2 : save(54, rho4 * rho4 * rho);
283 2 : save(53, (9 * rho2 - 8) * rho4 * rho2 * rho);
284 2 : save(52, ((36 * rho2 - 56) * rho2 + 21) * rho4 * rho);
285 2 : save(51, (((84 * rho2 - 168) * rho2 + 105) * rho2 - 20) * rho2 * rho);
286 2 : save(50, ((((126 * rho2 - 280) * rho2 + 210) * rho2 - 60) * rho2 + 5) * rho);
287 : libmesh_fallthrough();
288 :
289 2 : case 8:
290 2 : save(44, rho4 * rho4);
291 2 : save(43, (8 * rho2 - 7) * rho4 * rho2);
292 2 : save(42, ((28 * rho2 - 42) * rho2 + 15) * rho4);
293 2 : save(41, (((56 * rho2 - 105) * rho2 + 60) * rho2 - 10) * rho2);
294 2 : save(40, (((70 * rho2 - 140) * rho2 + 90) * rho2 - 20) * rho2 + 1);
295 : libmesh_fallthrough();
296 :
297 2 : case 7:
298 2 : save(35, rho4 * rho2 * rho);
299 2 : save(34, (7 * rho2 - 6) * rho4 * rho);
300 2 : save(33, ((21 * rho2 - 30) * rho2 + 10) * rho2 * rho);
301 2 : save(32, (((35 * rho2 - 60) * rho2 + 30) * rho2 - 4) * rho);
302 : libmesh_fallthrough();
303 :
304 2 : case 6:
305 2 : save(27, rho4 * rho2);
306 2 : save(26, (6 * rho2 - 5) * rho4);
307 2 : save(25, ((15 * rho2 - 20) * rho2 + 6) * rho2);
308 2 : save(24, ((20 * rho2 - 30) * rho2 + 12) * rho2 - 1);
309 : libmesh_fallthrough();
310 :
311 3 : case 5:
312 3 : save(20, rho4 * rho);
313 3 : save(19, (5 * rho2 - 4) * rho2 * rho);
314 3 : save(18, ((10 * rho2 - 12) * rho2 + 3) * rho);
315 : libmesh_fallthrough();
316 :
317 5 : case 4:
318 5 : save(14, rho4);
319 5 : save(13, (4 * rho2 - 3) * rho2);
320 5 : save(12, (6 * rho2 - 6) * rho2 + 1);
321 : libmesh_fallthrough();
322 :
323 5 : case 3:
324 5 : save(9, rho2 * rho);
325 5 : save(8, (3 * rho2 - 2) * rho);
326 : libmesh_fallthrough();
327 :
328 5 : case 2:
329 5 : save(5, rho2);
330 5 : save(4, 2 * rho2 - 1);
331 : libmesh_fallthrough();
332 :
333 5 : case 1:
334 5 : save(2, rho);
335 : libmesh_fallthrough();
336 :
337 5 : case 0:
338 5 : save(0, 1);
339 : }
340 :
341 19 : for (n = MAX_DIRECT_CALCULATION_ZERNIKE + 1; n <= _orders[0]; ++n)
342 : {
343 14 : j = simpleDoubleToSingle(n, n);
344 14 : save(j, pow(rho, n));
345 :
346 14 : j--;
347 14 : save(j, n * load(j + 1) - (n - 1) * load(j - (n + n)));
348 :
349 94 : for (q = n; q >= 4; q -= 2)
350 : {
351 80 : H3 = (-4 * (q - 2) * (q - 3)) / ((n + q - 2) * (n - q + 4.0));
352 80 : H2 = (H3 * (n + q) * (n - q + 2)) / (4.0 * (q - 1)) + (q - 2);
353 80 : H1 = q * (q - 1) / 2 - q * H2 + (H3 * (n + q + 2) * (n - q)) / 8.0;
354 80 : j--;
355 80 : save(j, H1 * load(j + 2) + (H2 + H3 / rho2) * load(j + 1));
356 : }
357 : }
358 :
359 5 : fillOutNegativeRankAndApplyAzimuthalComponent();
360 : }
361 :
362 : void
363 11 : Zernike::fillOutNegativeRankAndApplyAzimuthalComponent()
364 : {
365 : std::size_t n;
366 : long j, m, q, a;
367 : const Real & phi = _standardized_location[1];
368 :
369 : j = 0;
370 121 : for (n = 1; n <= _orders[0]; ++n)
371 : {
372 110 : j += n;
373 554 : for (m = 0, q = a = n; m < q; ++m, --q, a -= 2)
374 : {
375 444 : save(j + m, load(j + q) * sin(a * phi));
376 444 : save(j + q, load(j + q) * cos(a * phi));
377 : }
378 : }
379 11 : }
380 :
381 : const std::vector<Real> &
382 0 : Zernike::getStandardizedFunctionLimits() const
383 : {
384 : // Lazily instantiate the function limits array
385 0 : static const std::vector<Real> standardizedFunctionLimits = {0, 1, -M_PI, M_PI};
386 :
387 0 : return standardizedFunctionLimits;
388 : }
389 :
390 : Real
391 0 : Zernike::getStandardizedFunctionVolume() const
392 : {
393 0 : return M_PI; // The area of a unit disc is pi
394 : }
395 :
396 : std::vector<Real>
397 0 : Zernike::getStandardizedLocation(const std::vector<Real> & location) const
398 : {
399 : // Get the offset corresponding to the 'x' direction
400 0 : const Real offset1 = location[0] - _physical_bounds[0];
401 : // Get the offset corresponding to the 'y' direction
402 0 : const Real offset2 = location[1] - _physical_bounds[1];
403 : // Get the user-provided radius bound
404 : const Real & radius = _physical_bounds[2];
405 : // Covert to a radis and normalize
406 0 : const Real standardizedRadius = sqrt(offset1 * offset1 + offset2 * offset2) / radius;
407 : // Get the angle
408 0 : const Real theta = atan2(offset2, offset1);
409 :
410 0 : return {standardizedRadius, theta};
411 : }
412 :
413 : bool
414 0 : Zernike::isInPhysicalBounds(const Point & point) const
415 : {
416 : /*
417 : * Because Zernike polynomials live in RZ space, the easiest approach to check
418 : * this is to convert the physical location into a standardized location, then
419 : * check against the radius and theta bounds.
420 : */
421 0 : const std::vector<Real> location = extractLocationFromPoint(point);
422 0 : const std::vector<Real> standardized_location = getStandardizedLocation(location);
423 :
424 : /*
425 : * The radius (standardized_location[0]) is always positive, so only check
426 : * against the maximum radius (1). The theta components should always be in
427 : * bounds.
428 : */
429 0 : if (standardized_location[0] > 1.0)
430 : return false;
431 : else
432 0 : return true;
433 0 : }
434 :
435 : std::size_t
436 44 : Zernike::simpleDoubleToSingle(std::size_t n, long m) const
437 : {
438 44 : return (n * (n + 2) + m) / 2;
439 : }
|
