PRISMS-PF Manual
Loading...
Searching...
No Matches
symmetry.h
Go to the documentation of this file.
1// SPDX-FileCopyrightText: © 2026 PRISMS Center at the University of Michigan
2// SPDX-License-Identifier: GNU Lesser General Public Version 2.1
3
4#pragma once
5
6#include <deal.II/base/vectorization.h>
7
8#include <prismspf/utilities/vectorized_operations.h>
9
10#include <prismspf/config.h>
11
12PRISMS_PF_BEGIN_NAMESPACE
13
14namespace Symmetry
15{
28 template <unsigned int N, typename T>
29 [[nodiscard]] inline DEAL_II_ALWAYS_INLINE T
30 cos_arctan(const T &x) noexcept
31 {
32 // NOLINTBEGIN(readability-magic-numbers,cppcoreguidelines-avoid-magic-numbers)
33
34 using std::atan;
35 using std::cos;
36 using std::sqrt;
37
38 if constexpr (N <= 3)
39 {
40 const T t1 = T(1) / sqrt(T(1) + x * x); // cos θ
41 const T t2 = t1 * t1; // (cos θ)^2
42
43 if constexpr (N == 0)
44 {
45 return T(1);
46 }
47 else if constexpr (N == 1)
48 {
49 return t1;
50 }
51 else if constexpr (N == 2)
52 {
53 return T(2) * t2 - T(1);
54 }
55 else if constexpr (N == 3)
56 {
57 return t1 * (T(4) * t2 - T(3));
58 }
59 else if constexpr (N == 4)
60 {
61 return T(8) * t2 * t2 - T(8) * t2 + T(1);
62 }
63 }
64 else
65 {
66 return cos(T(N) * atan(x));
67 }
68
69 // NOLINTEND(readability-magic-numbers,cppcoreguidelines-avoid-magic-numbers)
70 }
71
84 template <unsigned int N, typename T>
85 [[nodiscard]] inline DEAL_II_ALWAYS_INLINE T
86 sin_arctan(const T &x) noexcept
87 {
88 // NOLINTBEGIN(readability-magic-numbers,cppcoreguidelines-avoid-magic-numbers)
89
90 using std::atan;
91 using std::sin;
92 using std::sqrt;
93
94 if constexpr (N <= 3)
95 {
96 const T t1 = T(1) / sqrt(T(1) + x * x); // cos θ
97 const T s1 = x * t1; // sin θ
98
99 if constexpr (N == 0)
100 {
101 return T(0);
102 }
103 else if constexpr (N == 1)
104 {
105 return s1;
106 }
107 else if constexpr (N == 2)
108 {
109 return T(2) * s1 * t1;
110 }
111 else if constexpr (N == 3)
112 {
113 return s1 * (T(3) - T(4) * s1 * s1);
114 }
115 else if constexpr (N == 4)
116 {
117 return T(4) * s1 * t1 * (T(1) - T(2) * s1 * s1);
118 }
119 }
120 else
121 {
122 return sin(T(N) * atan(x));
123 }
124
125 // NOLINTEND(readability-magic-numbers,cppcoreguidelines-avoid-magic-numbers)
126 }
127
140 template <unsigned int N, typename T>
141 [[nodiscard]] inline DEAL_II_ALWAYS_INLINE T
142 cos_theta(const T &nx, const T &ny) noexcept
143 {
144 // NOLINTBEGIN(readability-magic-numbers,cppcoreguidelines-avoid-magic-numbers)
145
146 using std::atan2;
147 using std::cos;
148
149 if constexpr (N <= 5)
150 {
151 const T t1 = nx; // cos θ
152 const T t2 = t1 * t1; // (cos θ)^2
153 const T t3 = t2 * t1; // (cos θ)^3
154 const T t4 = t2 * t2; // (cos θ)^4
155 const T t5 = t3 * t2; // (cos θ)^5
156 const T t6 = t3 * t3; // (cos θ)^6
157 const T t7 = t4 * t3; // (cos θ)^7
158 const T t8 = t4 * t4; // (cos θ)^8
159 const T t9 = t5 * t4; // (cos θ)^9
160 const T t10 = t5 * t5; // (cos θ)^10
161 const T t11 = t6 * t5; // (cos θ)^11
162 const T t12 = t6 * t6; // (cos θ)^12
163
164 if constexpr (N == 0)
165 {
166 return T(1);
167 }
168 else if constexpr (N == 1)
169 {
170 return t1;
171 }
172 else if constexpr (N == 2)
173 {
174 return T(2) * t2 - T(1);
175 }
176 else if constexpr (N == 3)
177 {
178 return T(4) * t3 - T(3) * t1;
179 }
180 else if constexpr (N == 4)
181 {
182 return T(8) * t4 - T(8) * t2 + T(1);
183 }
184 else if constexpr (N == 5)
185 {
186 return T(16) * t5 - T(20) * t3 + T(5) * t1;
187 }
188 else if constexpr (N == 6)
189 {
190 return T(32) * t6 - T(48) * t4 + T(18) * t2 - T(1);
191 }
192 else if constexpr (N == 7)
193 {
194 return T(64) * t7 - T(112) * t5 + T(56) * t3 - T(7) * t1;
195 }
196 else if constexpr (N == 8)
197 {
198 return T(128) * t8 - T(256) * t6 + T(160) * t4 - T(32) * t2 + T(1);
199 }
200 else if constexpr (N == 9)
201 {
202 return T(256) * t9 - T(576) * t7 + T(432) * t5 - T(120) * t3 + T(9) * t1;
203 }
204 else if constexpr (N == 10)
205 {
206 return T(512) * t10 - T(1280) * t8 + T(1120) * t6 - T(400) * t4 + T(50) * t2 -
207 T(1);
208 }
209 else if constexpr (N == 11)
210 {
211 return T(1024) * t11 - T(2816) * t9 + T(2816) * t7 - T(1232) * t5 +
212 T(220) * t3 - T(11) * t1;
213 }
214 else if constexpr (N == 12)
215 {
216 return T(2048) * t12 - T(6144) * t10 + T(6912) * t8 - T(3584) * t6 +
217 T(840) * t4 - T(72) * t2 + T(1);
218 }
219 }
220 else
221 {
222 return cos(T(N) * atan2(ny, nx));
223 }
224
225 // NOLINTEND(readability-magic-numbers,cppcoreguidelines-avoid-magic-numbers)
226 }
227
240 template <unsigned int N, typename T>
241 [[nodiscard]] inline DEAL_II_ALWAYS_INLINE T
242 sin_theta(const T &nx, const T &ny) noexcept
243 {
244 // NOLINTBEGIN(readability-magic-numbers,cppcoreguidelines-avoid-magic-numbers)
245
246 using std::atan2;
247 using std::sin;
248
249 if constexpr (N <= 5)
250 {
251 const T t1 = nx; // cos θ
252 const T s1 = ny; // sin θ
253 const T t2 = t1 * t1; // (cos θ)^2
254 const T t3 = t2 * t1; // (cos θ)^3
255 const T t4 = t2 * t2; // (cos θ)^4
256 const T t5 = t3 * t2; // (cos θ)^5
257 const T t6 = t3 * t3; // (cos θ)^6
258 const T t7 = t4 * t3; // (cos θ)^7
259 const T t8 = t4 * t4; // (cos θ)^8
260 const T t9 = t5 * t4; // (cos θ)^9
261 const T t10 = t5 * t5; // (cos θ)^10
262 const T t11 = t6 * t5; // (cos θ)^11
263
264 if constexpr (N == 0)
265 {
266 return T(0);
267 }
268 else if constexpr (N == 1)
269 {
270 return s1;
271 }
272 else if constexpr (N == 2)
273 {
274 return T(2) * s1 * t1;
275 }
276 else if constexpr (N == 3)
277 {
278 return s1 * (T(4) * t2 - T(1));
279 }
280 else if constexpr (N == 4)
281 {
282 return s1 * (T(8) * t3 - T(4) * t1);
283 }
284 else if constexpr (N == 5)
285 {
286 return s1 * (T(16) * t4 - T(12) * t2 + T(1));
287 }
288 else if constexpr (N == 6)
289 {
290 return s1 * (T(32) * t5 - T(32) * t3 + T(6) * t1);
291 }
292 else if constexpr (N == 7)
293 {
294 return s1 * (T(64) * t6 - T(80) * t4 + T(24) * t2 - T(1));
295 }
296 else if constexpr (N == 8)
297 {
298 return s1 * (T(128) * t7 - T(192) * t5 + T(80) * t3 - T(8) * t1);
299 }
300 else if constexpr (N == 9)
301 {
302 return s1 * (T(256) * t8 - T(448) * t6 + T(240) * t4 - T(40) * t2 + T(1));
303 }
304 else if constexpr (N == 10)
305 {
306 return s1 *
307 (T(512) * t9 - T(1024) * t7 + T(672) * t5 - T(160) * t3 + T(10) * t1);
308 }
309 else if constexpr (N == 11)
310 {
311 return s1 * (T(1024) * t10 - T(2304) * t8 + T(1792) * t6 - T(560) * t4 +
312 T(60) * t2 - T(1));
313 }
314 else if constexpr (N == 12)
315 {
316 return s1 * (T(2048) * t11 - T(5120) * t9 + T(4608) * t7 - T(1792) * t5 +
317 T(280) * t3 - T(12) * t1);
318 }
319 }
320 else
321 {
322 return sin(T(N) * atan2(ny, nx));
323 }
324
325 // NOLINTEND(readability-magic-numbers,cppcoreguidelines-avoid-magic-numbers)
326 }
327
340 template <unsigned int N, typename T>
341 [[nodiscard]] inline DEAL_II_ALWAYS_INLINE T
342 cos_psi(const T &nx, const T &ny, const T &nz) noexcept
343 {
344 // NOLINTBEGIN(readability-magic-numbers,cppcoreguidelines-avoid-magic-numbers)
345
346 using std::atan2;
347 using std::cos;
348 using std::sqrt;
349
350 if constexpr (N <= 5)
351 {
352 const T t1 = nz; // cos θ
353 const T t2 = t1 * t1; // (cos θ)^2
354 const T t3 = t2 * t1; // (cos θ)^3
355 const T t4 = t2 * t2; // (cos θ)^4
356 const T t5 = t3 * t2; // (cos θ)^5
357 const T t6 = t3 * t3; // (cos θ)^6
358 const T t7 = t4 * t3; // (cos θ)^7
359 const T t8 = t4 * t4; // (cos θ)^8
360 const T t9 = t5 * t4; // (cos θ)^9
361 const T t10 = t5 * t5; // (cos θ)^10
362 const T t11 = t6 * t5; // (cos θ)^11
363 const T t12 = t6 * t6; // (cos θ)^12
364
365 if constexpr (N == 0)
366 {
367 return T(1);
368 }
369 else if constexpr (N == 1)
370 {
371 return t1;
372 }
373 else if constexpr (N == 2)
374 {
375 return T(2) * t2 - T(1);
376 }
377 else if constexpr (N == 3)
378 {
379 return T(4) * t3 - T(3) * t1;
380 }
381 else if constexpr (N == 4)
382 {
383 return T(8) * t4 - T(8) * t2 + T(1);
384 }
385 else if constexpr (N == 5)
386 {
387 return T(16) * t5 - T(20) * t3 + T(5) * t1;
388 }
389 else if constexpr (N == 6)
390 {
391 return T(32) * t6 - T(48) * t4 + T(18) * t2 - T(1);
392 }
393 else if constexpr (N == 7)
394 {
395 return T(64) * t7 - T(112) * t5 + T(56) * t3 - T(7) * t1;
396 }
397 else if constexpr (N == 8)
398 {
399 return T(128) * t8 - T(256) * t6 + T(160) * t4 - T(32) * t2 + T(1);
400 }
401 else if constexpr (N == 9)
402 {
403 return T(256) * t9 - T(576) * t7 + T(432) * t5 - T(120) * t3 + T(9) * t1;
404 }
405 else if constexpr (N == 10)
406 {
407 return T(512) * t10 - T(1280) * t8 + T(1120) * t6 - T(400) * t4 + T(50) * t2 -
408 T(1);
409 }
410 else if constexpr (N == 11)
411 {
412 return T(1024) * t11 - T(2816) * t9 + T(2816) * t7 - T(1232) * t5 +
413 T(220) * t3 - T(11) * t1;
414 }
415 else if constexpr (N == 12)
416 {
417 return T(2048) * t12 - T(6144) * t10 + T(6912) * t8 - T(3584) * t6 +
418 T(840) * t4 - T(72) * t2 + T(1);
419 }
420 }
421 else
422 {
423 return cos(T(N) * atan2(sqrt(nx * nx + ny * ny), nz));
424 }
425
426 // NOLINTEND(readability-magic-numbers,cppcoreguidelines-avoid-magic-numbers)
427 }
428
441 template <unsigned int N, typename T>
442 [[nodiscard]] inline DEAL_II_ALWAYS_INLINE T
443 sin_psi(const T &nx, const T &ny, const T &nz) noexcept
444 {
445 // NOLINTBEGIN(readability-magic-numbers,cppcoreguidelines-avoid-magic-numbers)
446
447 using std::atan2;
448 using std::sin;
449 using std::sqrt;
450
451 if constexpr (N <= 5)
452 {
453 const T s1 = sqrt(nx * nx + ny * ny); // sin θ
454 const T t1 = nz; // cos θ
455 const T t2 = t1 * t1; // (cos θ)^2
456 const T t3 = t2 * t1; // (cos θ)^3
457 const T t4 = t2 * t2; // (cos θ)^4
458 const T t5 = t3 * t2; // (cos θ)^5
459 const T t6 = t3 * t3; // (cos θ)^6
460 const T t7 = t4 * t3; // (cos θ)^7
461 const T t8 = t4 * t4; // (cos θ)^8
462 const T t9 = t5 * t4; // (cos θ)^9
463 const T t10 = t5 * t5; // (cos θ)^10
464 const T t11 = t6 * t5; // (cos θ)^11
465
466 if constexpr (N == 0)
467 {
468 return T(0);
469 }
470 else if constexpr (N == 1)
471 {
472 return s1;
473 }
474 else if constexpr (N == 2)
475 {
476 return T(2) * s1 * t1;
477 }
478 else if constexpr (N == 3)
479 {
480 return s1 * (T(4) * t2 - T(1));
481 }
482 else if constexpr (N == 4)
483 {
484 return s1 * (T(8) * t3 - T(4) * t1);
485 }
486 else if constexpr (N == 5)
487 {
488 return s1 * (T(16) * t4 - T(12) * t2 + T(1));
489 }
490 else if constexpr (N == 6)
491 {
492 return s1 * (T(32) * t5 - T(32) * t3 + T(6) * t1);
493 }
494 else if constexpr (N == 7)
495 {
496 return s1 * (T(64) * t6 - T(80) * t4 + T(24) * t2 - T(1));
497 }
498 else if constexpr (N == 8)
499 {
500 return s1 * (T(128) * t7 - T(192) * t5 + T(80) * t3 - T(8) * t1);
501 }
502 else if constexpr (N == 9)
503 {
504 return s1 * (T(256) * t8 - T(448) * t6 + T(240) * t4 - T(40) * t2 + T(1));
505 }
506 else if constexpr (N == 10)
507 {
508 return s1 *
509 (T(512) * t9 - T(1024) * t7 + T(672) * t5 - T(160) * t3 + T(10) * t1);
510 }
511 else if constexpr (N == 11)
512 {
513 return s1 * (T(1024) * t10 - T(2304) * t8 + T(1792) * t6 - T(560) * t4 +
514 T(60) * t2 - T(1));
515 }
516 else if constexpr (N == 12)
517 {
518 return s1 * (T(2048) * t11 - T(5120) * t9 + T(4608) * t7 - T(1792) * t5 +
519 T(280) * t3 - T(12) * t1);
520 }
521 }
522 else
523 {
524 return sin(T(N) * atan2(sqrt(nx * nx + ny * ny), nz));
525 }
526
527 // NOLINTEND(readability-magic-numbers,cppcoreguidelines-avoid-magic-numbers)
528 }
529
530} // namespace Symmetry
531
532PRISMS_PF_END_NAMESPACE
PRISMS_PF_BEGIN_NAMESPACE
Definition conditional_ostreams.cc:20
Symmetry
Definition symmetry.h:15
Symmetry::sin_theta
DEAL_II_ALWAYS_INLINE T sin_theta(const T &nx, const T &ny) noexcept
Compute sin(N * theta) where theta = arctan(ny/nx).
Definition symmetry.h:242
Symmetry::cos_arctan
DEAL_II_ALWAYS_INLINE T cos_arctan(const T &x) noexcept
Compute cos(N * arctan(x)).
Definition symmetry.h:30
Symmetry::sin_arctan
DEAL_II_ALWAYS_INLINE T sin_arctan(const T &x) noexcept
Compute sin(N * arctan(x)).
Definition symmetry.h:86
Symmetry::sin_psi
DEAL_II_ALWAYS_INLINE T sin_psi(const T &nx, const T &ny, const T &nz) noexcept
Compute sin(N * psi) where theta = arctan(sqrt(nx^2+ny^2)/ny).
Definition symmetry.h:443
Symmetry::cos_psi
DEAL_II_ALWAYS_INLINE T cos_psi(const T &nx, const T &ny, const T &nz) noexcept
Compute cos(N * psi) where theta = arctan(sqrt(nx^2+ny^2)/ny).
Definition symmetry.h:342
Symmetry::cos_theta
DEAL_II_ALWAYS_INLINE T cos_theta(const T &nx, const T &ny) noexcept
Compute cos(N * theta) where theta = arctan(ny/nx).
Definition symmetry.h:142
std::atan2
inline ::dealii::VectorizedArray< Number, width > atan2(const ::dealii::VectorizedArray< Number, width > &y, const ::dealii::VectorizedArray< Number, width > &x)
Definition vectorized_operations.h:43
vectorized_operations.h