Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <type_traits>
8 :
9 : #include "NumericalAlgorithms/Spectral/Basis.hpp"
10 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
11 : #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
12 : #include "Utilities/ErrorHandling/Error.hpp"
13 :
14 : namespace Spectral::detail {
15 : template <typename F>
16 : decltype(auto) get_spectral_quantity_for_mesh(F&& f, const Mesh<1>& mesh) {
17 : const auto num_points = mesh.extents(0);
18 : // Switch on runtime values of basis and quadrature to select
19 : // corresponding template specialization. For basis functions spanning
20 : // multiple dimensions we can generalize this function to take a
21 : // higher-dimensional Mesh.
22 : switch (mesh.basis(0)) {
23 : case Basis::Legendre:
24 : switch (mesh.quadrature(0)) {
25 : case Quadrature::Gauss:
26 : return f(std::integral_constant<Basis, Basis::Legendre>{},
27 : std::integral_constant<Quadrature, Quadrature::Gauss>{},
28 : num_points);
29 : case Quadrature::GaussLobatto:
30 : return f(
31 : std::integral_constant<Basis, Basis::Legendre>{},
32 : std::integral_constant<Quadrature, Quadrature::GaussLobatto>{},
33 : num_points);
34 : default:
35 : ERROR("Missing quadrature case for spectral quantity");
36 : }
37 : case Basis::Chebyshev:
38 : switch (mesh.quadrature(0)) {
39 : case Quadrature::Gauss:
40 : return f(std::integral_constant<Basis, Basis::Chebyshev>{},
41 : std::integral_constant<Quadrature, Quadrature::Gauss>{},
42 : num_points);
43 : case Quadrature::GaussLobatto:
44 : return f(
45 : std::integral_constant<Basis, Basis::Chebyshev>{},
46 : std::integral_constant<Quadrature, Quadrature::GaussLobatto>{},
47 : num_points);
48 : default:
49 : ERROR("Missing quadrature case for spectral quantity");
50 : }
51 : case Basis::Cartoon:
52 : switch (mesh.quadrature(0)) {
53 : case Quadrature::AxialSymmetry:
54 : return f(
55 : std::integral_constant<Basis, Basis::Cartoon>{},
56 : std::integral_constant<Quadrature, Quadrature::AxialSymmetry>{},
57 : num_points);
58 : case Quadrature::SphericalSymmetry:
59 : return f(std::integral_constant<Basis, Basis::Cartoon>{},
60 : std::integral_constant<Quadrature,
61 : Quadrature::SphericalSymmetry>{},
62 : num_points);
63 : default:
64 : ERROR(
65 : "Only Axial and Spherical Symmetry quadratures are allowed for "
66 : "a Cartoon basis.");
67 : }
68 : case Basis::Fourier:
69 : switch (mesh.quadrature(0)) {
70 : case Quadrature::Equiangular:
71 : return f(
72 : std::integral_constant<Basis, Basis::Fourier>{},
73 : std::integral_constant<Quadrature, Quadrature::Equiangular>{},
74 : num_points);
75 : default:
76 : ERROR("Missing quadrature case for spectral quantity");
77 : }
78 : case Basis::ZernikeB1:
79 : switch (mesh.quadrature(0)) {
80 : case Quadrature::GaussRadauUpper:
81 : return f(
82 : std::integral_constant<Basis, Basis::ZernikeB1>{},
83 : std::integral_constant<Quadrature, Quadrature::GaussRadauUpper>{},
84 : num_points);
85 : default:
86 : ERROR("Missing quadrature case for spectral quantity");
87 : }
88 : case Basis::ZernikeB2:
89 : switch (mesh.quadrature(0)) {
90 : case Quadrature::GaussRadauUpper:
91 : return f(
92 : std::integral_constant<Basis, Basis::ZernikeB2>{},
93 : std::integral_constant<Quadrature, Quadrature::GaussRadauUpper>{},
94 : num_points);
95 : case Quadrature::Equiangular:
96 : // While ZernikeB2 with Equiangular is treated differently in
97 : // operators (e.g. partial derivatives, interpolation), in terms of
98 : // its spectral quantities it is Fourier with Equiangular
99 : return f(
100 : std::integral_constant<Basis, Basis::Fourier>{},
101 : std::integral_constant<Quadrature, Quadrature::Equiangular>{},
102 : num_points);
103 : default:
104 : ERROR("Missing quadrature case for spectral quantity");
105 : }
106 : case Basis::ZernikeB3:
107 : switch (mesh.quadrature(0)) {
108 : case Quadrature::GaussRadauUpper:
109 : return f(
110 : std::integral_constant<Basis, Basis::ZernikeB3>{},
111 : std::integral_constant<Quadrature, Quadrature::GaussRadauUpper>{},
112 : num_points);
113 : default:
114 : ERROR("Missing quadrature case for spectral quantity");
115 : }
116 : case Basis::FiniteDifference:
117 : switch (mesh.quadrature(0)) {
118 : case Quadrature::CellCentered:
119 : return f(
120 : std::integral_constant<Basis, Basis::FiniteDifference>{},
121 : std::integral_constant<Quadrature, Quadrature::CellCentered>{},
122 : num_points);
123 : case Quadrature::FaceCentered:
124 : return f(
125 : std::integral_constant<Basis, Basis::FiniteDifference>{},
126 : std::integral_constant<Quadrature, Quadrature::FaceCentered>{},
127 : num_points);
128 : default:
129 : ERROR(
130 : "Only CellCentered and FaceCentered are supported for finite "
131 : "difference quadrature.");
132 : }
133 : case Basis::SphericalHarmonic:
134 : ERROR(
135 : "Basis::SphericalHarmonic is a two-dimensional basis and is not "
136 : "supported for this function. If you want the collocation points, "
137 : "use the function logical_coordinates.");
138 : default:
139 : ERROR("Missing basis case for spectral quantity. The missing basis is: "
140 : << mesh.basis(0));
141 : }
142 : }
143 :
144 : template <typename F>
145 : decltype(auto) get_two_indexed_spectral_quantity_for_mesh(F&& f,
146 : const Mesh<1>& mesh,
147 : const size_t m,
148 : const size_t N) {
149 : const auto num_points = mesh.extents(0);
150 : switch (mesh.basis(0)) {
151 : case Basis::ZernikeB1:
152 : switch (mesh.quadrature(0)) {
153 : case Quadrature::GaussRadauUpper:
154 : return f(
155 : std::integral_constant<Basis, Basis::ZernikeB1>{},
156 : std::integral_constant<Quadrature, Quadrature::GaussRadauUpper>{},
157 : num_points, m, N);
158 : default:
159 : ERROR("Missing quadrature case for two-indexed spectral quantity");
160 : }
161 : case Basis::ZernikeB2:
162 : switch (mesh.quadrature(0)) {
163 : case Quadrature::GaussRadauUpper:
164 : return f(
165 : std::integral_constant<Basis, Basis::ZernikeB2>{},
166 : std::integral_constant<Quadrature, Quadrature::GaussRadauUpper>{},
167 : num_points, m, N);
168 : default:
169 : ERROR("Missing quadrature case for two-indexed spectral quantity");
170 : }
171 : case Basis::ZernikeB3:
172 : switch (mesh.quadrature(0)) {
173 : case Quadrature::GaussRadauUpper:
174 : return f(
175 : std::integral_constant<Basis, Basis::ZernikeB3>{},
176 : std::integral_constant<Quadrature, Quadrature::GaussRadauUpper>{},
177 : num_points, m, N);
178 : default:
179 : ERROR("Missing quadrature case for two-indexed spectral quantity");
180 : }
181 : default:
182 : ERROR(
183 : "Missing basis case for two-indexed spectral quantity. The "
184 : "missing basis is: "
185 : << mesh.basis(0));
186 : }
187 : }
188 :
189 : } // namespace Spectral::detail
|
