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