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