www.mooseframework.org
Zernike.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
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 
18 #define MAX_DIRECT_CALCULATION_ZERNIKE 10
19 
20 Zernike::Zernike() : SingleSeriesBasisInterface() {}
21 
22 Zernike::Zernike(const std::vector<MooseEnum> & domain,
23  const std::vector<std::size_t> & order,
24  MooseEnum expansion_type,
25  MooseEnum generation_type)
26  : SingleSeriesBasisInterface(domain, order, calculatedNumberOfTermsBasedOnOrder(order))
27 {
28  if (expansion_type == "orthonormal")
29  _evaluateExpansionWrapper = [this]() { this->evaluateOrthonormal(); };
30  else if (expansion_type == "sqrt_mu")
31  _evaluateExpansionWrapper = [this]() { this->evaluateSqrtMu(); };
32  else if (expansion_type == "standard")
33  _evaluateExpansionWrapper = [this]() { this->evaluateStandard(); };
34  else
35  mooseError("The specified type of normalization for expansion does not exist");
36 
37  if (generation_type == "orthonormal")
38  _evaluateGenerationWrapper = [this]() { this->evaluateOrthonormal(); };
39  else if (generation_type == "sqrt_mu")
40  _evaluateGenerationWrapper = [this]() { this->evaluateSqrtMu(); };
41  else if (generation_type == "standard")
42  _evaluateGenerationWrapper = [this]() { this->evaluateStandard(); };
43  else
44  mooseError("The specified type of normalization for generation does not exist");
45 }
46 
47 std::size_t
48 Zernike::calculatedNumberOfTermsBasedOnOrder(const std::vector<std::size_t> & order) const
49 {
50  return ((order[0] + 1) * (order[0] + 2)) / 2;
51 }
52 
53 void
54 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  if (bounds.size() != 3)
61  mooseError("Zernike: Invalid number of bounds specified for single series!");
62 }
63 
64 // clang-format off
65 void
66 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  const Real rho2 = rho * rho;
73  const Real rho4 = rho2 * rho2;
74 
75  if (MooseUtils::absoluteFuzzyEqual(rho, 0.0))
76  {
77  for (n = 0; n <= _orders[0]; n += 2)
78  {
79  j = simpleDoubleToSingle(n, 0);
80 
81  if ((n / 2) % 2 != 0)
82  save(j, -1 * (n + 1) / M_PI);
83  else
84  save(j, 1 * (n + 1) / M_PI);
85  }
86 
87  return;
88  }
89 
90  switch (_orders[0])
91  {
92  default:
93  case MAX_DIRECT_CALCULATION_ZERNIKE: /* 10 */
94  save(65, rho4 * rho4 * rho2
95  * 22 / M_PI);
96  save(64, (10 * rho2 - 9) * rho4 * rho4
97  * 22 / M_PI);
98  save(63, ((45 * rho2 - 72) * rho2 + 28) * rho4 * rho2
99  * 22 / M_PI);
100  save(62, (((120 * rho2 - 252) * rho2 + 168) * rho2- 35) * rho4
101  * 22 / M_PI);
102  save(61, ((((210 * rho2 - 504) * rho2 + 420) * rho2 - 140) * rho2 + 15) * rho2
103  * 22 / M_PI);
104  save(60, (((((252 * rho2 - 630) * rho2 + 560) * rho2 - 210) * rho2 + 30) * rho2 - 1)
105  * 11 / M_PI);
106  libmesh_fallthrough();
107 
108  case 9:
109  save(54, rho4 * rho4 * rho
110  * 20 / M_PI);
111  save(53, (9 * rho2 - 8) * rho4 * rho2 * rho
112  * 20 / M_PI);
113  save(52, ((36 * rho2 - 56) * rho2 + 21) * rho4 * rho
114  * 20 / M_PI);
115  save(51, (((84 * rho2 - 168) * rho2 + 105) * rho2 - 20) * rho2 * rho
116  * 20 / M_PI);
117  save(50, ((((126 * rho2 - 280) * rho2 + 210) * rho2 - 60) * rho2 + 5) * rho
118  * 20 / M_PI);
119  libmesh_fallthrough();
120 
121  case 8:
122  save(44, rho4 * rho4
123  * 18 / M_PI);
124  save(43, (8 * rho2 - 7) * rho4 * rho2
125  * 18 / M_PI);
126  save(42, ((28 * rho2 - 42) * rho2 + 15) * rho4
127  * 18 / M_PI);
128  save(41, (((56 * rho2 - 105) * rho2 + 60) * rho2 - 10) * rho2
129  * 18 / M_PI);
130  save(40, ((((70 * rho2 - 140) * rho2 + 90) * rho2 - 20) * rho2 + 1)
131  * 9 / M_PI);
132  libmesh_fallthrough();
133 
134  case 7:
135  save(35, rho4 * rho2 * rho
136  * 16 / M_PI);
137  save(34, (7 * rho2 - 6) * rho4 * rho
138  * 16 / M_PI);
139  save(33, ((21 * rho2 - 30) * rho2 + 10) * rho2 * rho
140  * 16 / M_PI);
141  save(32, (((35 * rho2 - 60) * rho2 + 30) * rho2 - 4) * rho
142  * 16 / M_PI);
143  libmesh_fallthrough();
144 
145  case 6:
146  save(27, rho4 * rho2
147  * 14 / M_PI);
148  save(26, (6 * rho2 - 5) * rho4
149  * 14 / M_PI);
150  save(25, ((15 * rho2 - 20) * rho2 + 6) * rho2
151  * 14 / M_PI);
152  save(24, (((20 * rho2 - 30) * rho2 + 12) * rho2 - 1)
153  * 7 / M_PI);
154  libmesh_fallthrough();
155 
156  case 5:
157  save(20, rho4 * rho
158  * 12 / M_PI);
159  save(19, (5 * rho2 - 4) * rho2 * rho
160  * 12 / M_PI);
161  save(18, ((10 * rho2 - 12) * rho2 + 3) * rho
162  * 12 / M_PI);
163  libmesh_fallthrough();
164 
165  case 4:
166  save(14, rho4
167  * 10 / M_PI);
168  save(13, (4 * rho2 - 3) * rho2
169  * 10 / M_PI);
170  save(12, ((6 * rho2 - 6) * rho2 + 1)
171  * 5 / M_PI);
172  libmesh_fallthrough();
173 
174  case 3:
175  save(9, rho2 * rho
176  * 8 / M_PI);
177  save(8, (3 * rho2 - 2) * rho
178  * 8 / M_PI);
179  libmesh_fallthrough();
180 
181  case 2:
182  save(5, rho2
183  * 6 / M_PI);
184  save(4, (2 * rho2 - 1)
185  * 3 / M_PI);
186  libmesh_fallthrough();
187 
188  case 1:
189  save(2, rho
190  * 4 / M_PI);
191  libmesh_fallthrough();
192 
193  case 0:
194  save(0, 1
195  * 1 / M_PI);
196  }
197 
198  for (n = MAX_DIRECT_CALCULATION_ZERNIKE + 1; n <= _orders[0]; ++n)
199  {
200  j = simpleDoubleToSingle(n, n);
201  save(j, pow(rho, n)
202  * (n + n + 2) / M_PI);
203 
204  j--;
205  save(j, n * load(j + 1) - (n + 1) * load(j - (n + n)));
206 
207  for (q = n; q >= 4; q -= 2)
208  {
209  H3 = (-4 * (q - 2) * (q - 3)) / ((n + q - 2) * (n - q + 4.0));
210  H2 = (H3 * (n + q) * (n - q + 2)) / (4.0 * (q - 1)) + (q - 2);
211  H1 = q * (q - 1) / 2 - q * H2 + (H3 * (n + q + 2) * (n - q)) / 8.0;
212  j--;
213  if (q == 4)
214  save(j, (H1 * load(j + 2) + (H2 + H3 / rho2) * load(j + 1))
215  * 0.5);
216  else
217  save(j, H1 * load(j + 2) + (H2 + H3 / rho2) * load(j + 1));
218  }
219  }
220 
221  fillOutNegativeRankAndApplyAzimuthalComponent();
222 }
223 // clang-format on
224 
225 void
226 Zernike::evaluateSqrtMu()
227 {
228  evaluateStandard();
229  save(0, load(0) / std::sqrt(M_PI));
230  size_t j = 1;
231  for (size_t n = 1; n < _orders[0] + 1; ++n)
232  {
233  for (size_t m = 0; m < n + 1; ++m)
234  {
235  if (m != 0 && n / m == 2 && n % m == 0)
236  save(j, load(j) * std::sqrt((n + 1) / M_PI));
237  else
238  save(j, load(j) * std::sqrt((2 * n + 2) / M_PI));
239  ++j;
240  }
241  }
242 }
243 
244 void
245 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  const Real rho2 = rho * rho;
252  const Real rho4 = rho2 * rho2;
253 
254  if (MooseUtils::absoluteFuzzyLessEqual(rho, 0))
255  {
256  for (n = 0; n <= _orders[0]; n += 2)
257  {
258  j = simpleDoubleToSingle(n, 0);
259 
260  if ((n / 2) % 2 != 0)
261  save(j, -1);
262  else
263  save(j, 1);
264  }
265 
266  return;
267  }
268 
269  switch (_orders[0])
270  {
271  default:
272  case MAX_DIRECT_CALCULATION_ZERNIKE: /* 10 */
273  save(65, rho4 * rho4 * rho2);
274  save(64, (10 * rho2 - 9) * rho4 * rho4);
275  save(63, ((45 * rho2 - 72) * rho2 + 28) * rho4 * rho2);
276  save(62, (((120 * rho2 - 252) * rho2 + 168) * rho2 - 35) * rho4);
277  save(61, ((((210 * rho2 - 504) * rho2 + 420) * rho2 - 140) * rho2 + 15) * rho2);
278  save(60, ((((252 * rho2 - 630) * rho2 + 560) * rho2 - 210) * rho2 + 30) * rho2 - 1);
279  libmesh_fallthrough();
280 
281  case 9:
282  save(54, rho4 * rho4 * rho);
283  save(53, (9 * rho2 - 8) * rho4 * rho2 * rho);
284  save(52, ((36 * rho2 - 56) * rho2 + 21) * rho4 * rho);
285  save(51, (((84 * rho2 - 168) * rho2 + 105) * rho2 - 20) * rho2 * rho);
286  save(50, ((((126 * rho2 - 280) * rho2 + 210) * rho2 - 60) * rho2 + 5) * rho);
287  libmesh_fallthrough();
288 
289  case 8:
290  save(44, rho4 * rho4);
291  save(43, (8 * rho2 - 7) * rho4 * rho2);
292  save(42, ((28 * rho2 - 42) * rho2 + 15) * rho4);
293  save(41, (((56 * rho2 - 105) * rho2 + 60) * rho2 - 10) * rho2);
294  save(40, (((70 * rho2 - 140) * rho2 + 90) * rho2 - 20) * rho2 + 1);
295  libmesh_fallthrough();
296 
297  case 7:
298  save(35, rho4 * rho2 * rho);
299  save(34, (7 * rho2 - 6) * rho4 * rho);
300  save(33, ((21 * rho2 - 30) * rho2 + 10) * rho2 * rho);
301  save(32, (((35 * rho2 - 60) * rho2 + 30) * rho2 - 4) * rho);
302  libmesh_fallthrough();
303 
304  case 6:
305  save(27, rho4 * rho2);
306  save(26, (6 * rho2 - 5) * rho4);
307  save(25, ((15 * rho2 - 20) * rho2 + 6) * rho2);
308  save(24, ((20 * rho2 - 30) * rho2 + 12) * rho2 - 1);
309  libmesh_fallthrough();
310 
311  case 5:
312  save(20, rho4 * rho);
313  save(19, (5 * rho2 - 4) * rho2 * rho);
314  save(18, ((10 * rho2 - 12) * rho2 + 3) * rho);
315  libmesh_fallthrough();
316 
317  case 4:
318  save(14, rho4);
319  save(13, (4 * rho2 - 3) * rho2);
320  save(12, (6 * rho2 - 6) * rho2 + 1);
321  libmesh_fallthrough();
322 
323  case 3:
324  save(9, rho2 * rho);
325  save(8, (3 * rho2 - 2) * rho);
326  libmesh_fallthrough();
327 
328  case 2:
329  save(5, rho2);
330  save(4, 2 * rho2 - 1);
331  libmesh_fallthrough();
332 
333  case 1:
334  save(2, rho);
335  libmesh_fallthrough();
336 
337  case 0:
338  save(0, 1);
339  }
340 
341  for (n = MAX_DIRECT_CALCULATION_ZERNIKE + 1; n <= _orders[0]; ++n)
342  {
343  j = simpleDoubleToSingle(n, n);
344  save(j, pow(rho, n));
345 
346  j--;
347  save(j, n * load(j + 1) - (n - 1) * load(j - (n + n)));
348 
349  for (q = n; q >= 4; q -= 2)
350  {
351  H3 = (-4 * (q - 2) * (q - 3)) / ((n + q - 2) * (n - q + 4.0));
352  H2 = (H3 * (n + q) * (n - q + 2)) / (4.0 * (q - 1)) + (q - 2);
353  H1 = q * (q - 1) / 2 - q * H2 + (H3 * (n + q + 2) * (n - q)) / 8.0;
354  j--;
355  save(j, H1 * load(j + 2) + (H2 + H3 / rho2) * load(j + 1));
356  }
357  }
358 
359  fillOutNegativeRankAndApplyAzimuthalComponent();
360 }
361 
362 void
363 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  for (n = 1; n <= _orders[0]; ++n)
371  {
372  j += n;
373  for (m = 0, q = a = n; m < q; ++m, --q, a -= 2)
374  {
375  save(j + m, load(j + q) * sin(a * phi));
376  save(j + q, load(j + q) * cos(a * phi));
377  }
378  }
379 }
380 
381 const std::vector<Real> &
382 Zernike::getStandardizedFunctionLimits() const
383 {
384  // Lazily instantiate the function limits array
385  static const std::vector<Real> standardizedFunctionLimits = {0, 1, -M_PI, M_PI};
386 
387  return standardizedFunctionLimits;
388 }
389 
390 Real
391 Zernike::getStandardizedFunctionVolume() const
392 {
393  return M_PI; // The area of a unit disc is pi
394 }
395 
396 std::vector<Real>
397 Zernike::getStandardizedLocation(const std::vector<Real> & location) const
398 {
399  // Get the offset corresponding to the 'x' direction
400  const Real offset1 = location[0] - _physical_bounds[0];
401  // Get the offset corresponding to the 'y' direction
402  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  const Real standardizedRadius = sqrt(offset1 * offset1 + offset2 * offset2) / radius;
407  // Get the angle
408  const Real theta = atan2(offset2, offset1);
409 
410  return {standardizedRadius, theta};
411 }
412 
413 bool
414 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  const std::vector<Real> location = extractLocationFromPoint(point);
422  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  if (standardized_location[0] > 1.0)
430  return false;
431  else
432  return true;
433 }
434 
435 std::size_t
436 Zernike::simpleDoubleToSingle(std::size_t n, long m) const
437 {
438  return (n * (n + 2) + m) / 2;
439 }
Zernike::getStandardizedLocation
virtual std::vector< Real > getStandardizedLocation(const std::vector< Real > &location) const override
Standardize the location according to the requirements of the underlying basis, which may actually co...
Definition: Zernike.C:397
SingleSeriesBasisInterface::_evaluateGenerationWrapper
std::function< void()> _evaluateGenerationWrapper
The generation evaluation wrapper.
Definition: SingleSeriesBasisInterface.h:91
Zernike::simpleDoubleToSingle
std::size_t simpleDoubleToSingle(std::size_t n, long m) const
Maps the double order/rank idices to a single linear index.
Definition: Zernike.C:436
SingleSeriesBasisInterface::_evaluateExpansionWrapper
std::function< void()> _evaluateExpansionWrapper
The expansion evaluation wrapper.
Definition: SingleSeriesBasisInterface.h:88
pow
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
Definition: ExpressionBuilder.h:673
Zernike::calculatedNumberOfTermsBasedOnOrder
virtual std::size_t calculatedNumberOfTermsBasedOnOrder(const std::vector< std::size_t > &order) const override
Returns the number of terms in the single series given a value for the order.
Definition: Zernike.C:48
FunctionalBasisInterface::load
Real load(std::size_t index) const
Helper function to load a value from #_series.
Definition: FunctionalBasisInterface.C:132
Zernike::checkPhysicalBounds
virtual void checkPhysicalBounds(const std::vector< Real > &bounds) const override
Checks the physical bounds according to the actual implementation.
Definition: Zernike.C:54
expansion_type
MooseEnum expansion_type
Zernike::evaluateStandard
virtual void evaluateStandard()
Evaluates the standard form of the basis functions.
Definition: Zernike.C:245
SingleSeriesBasisInterface::_standardized_location
std::vector< Real > _standardized_location
The standardized location of evaluation.
Definition: SingleSeriesBasisInterface.h:85
SingleSeriesBasisInterface
This class is a simple wrapper around FunctionalBasisInterface, and intended for use by any single fu...
Definition: SingleSeriesBasisInterface.h:19
Zernike::evaluateSqrtMu
virtual void evaluateSqrtMu()
Evaluates the 1/sqrt(mu) normalized form of the basis functions.
Definition: Zernike.C:226
generation_type
MooseEnum generation_type
SingleSeriesBasisInterface::_physical_bounds
std::vector< Real > _physical_bounds
The physical bounds of the series.
Definition: SingleSeriesBasisInterface.h:82
Zernike::Zernike
Zernike()
Definition: Zernike.C:20
Zernike::getStandardizedFunctionVolume
virtual Real getStandardizedFunctionVolume() const override
Returns the volume within the standardized function local_limits.
Definition: Zernike.C:391
Zernike.h
Zernike::evaluateOrthonormal
virtual void evaluateOrthonormal()
Evaluates the orthonormal form of the basis functions.
Definition: Zernike.C:66
FunctionalBasisInterface::save
void save(std::size_t index, Real value)
Helper function to store a value in #_series.
Definition: FunctionalBasisInterface.C:138
Zernike::isInPhysicalBounds
virtual bool isInPhysicalBounds(const Point &point) const override
Determines if the point provided is in within the physical bounds.
Definition: Zernike.C:414
SingleSeriesBasisInterface::_orders
std::vector< std::size_t > _orders
The order of the series.
Definition: SingleSeriesBasisInterface.h:79
SingleSeriesBasisInterface::extractLocationFromPoint
std::vector< Real > extractLocationFromPoint(const Point &point) const
Convert a spatial point to a location that the series will use to determine the value at which to eva...
Definition: SingleSeriesBasisInterface.C:122
Zernike::getStandardizedFunctionLimits
virtual const std::vector< Real > & getStandardizedFunctionLimits() const override
Returns a vector of the lower and upper bounds of the standard functional space.
Definition: Zernike.C:382
Zernike::fillOutNegativeRankAndApplyAzimuthalComponent
void fillOutNegativeRankAndApplyAzimuthalComponent()
Helper function used by evaluateGeneration() and evaluateExpansion().
Definition: Zernike.C:363
Generated on Sat Jan 25 2020 12:04:37 for www.mooseframework.org by   doxygen 1.8.16