Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <array>
7 : #include <cmath>
8 : #include <cstddef>
9 : #include <utility>
10 : #include <vector>
11 :
12 : #include "DataStructures/DataVector.hpp"
13 : #include "DataStructures/DynamicBuffer.hpp"
14 : #include "DataStructures/Tensor/TypeAliases.hpp"
15 : #include "NumericalAlgorithms/SphericalHarmonics/SpherepackHelper.hpp"
16 : #include "Utilities/Blas.hpp"
17 : #include "Utilities/ForceInline.hpp"
18 : #include "Utilities/Gsl.hpp"
19 :
20 : /// Items related to spherical harmonics
21 : namespace ylm {
22 :
23 : /*!
24 : * \ingroup SpectralGroup
25 : *
26 : * \brief Defines the C++ interface to SPHEREPACK.
27 : *
28 : * \details The class `Spherepack` defines the C++ interface to the fortran
29 : * library SPHEREPACK used for computations on the surface of a sphere.
30 : *
31 : * Given a real-valued, scalar function \f$g(\theta, \phi)\f$, SPHEREPACK
32 : * expands it as:
33 : *
34 : * \f{align}
35 : * g(\theta, \phi)
36 : * &=\frac{1}{2}\sum_{l=0}^{l_{\max}}\bar P_l^0(\cos\theta) a_{l0}
37 : * +\sum_{l=1}^{l_{\max}}\sum_{m=1}^{\min(l, m_{\max})}\bar P_l^m(\cos\theta)\{
38 : * a_{lm}\cos m\phi -b_{lm}\sin m\phi\}\label{eq:spherepack_expansion}
39 : * \f}
40 : *
41 : * where \f$a_{lm}\f$ and \f$b_{lm}\f$ are real-valued
42 : * spectral coefficient arrays used by
43 : * SPHEREPACK, \f$P_l^m(x)\f$ are defined as
44 : *
45 : * \f{align}
46 : * \bar P_l^m(x)&=\sqrt{\frac{(2l+1)(l-m)!}{2(l+m)!}}\;P_{lm}(x)
47 : * \f}
48 : *
49 : * and \f$P_{nm}(x)\f$ are the associated Legendre polynomials as defined,
50 : * for example, in Jackson's "Classical Electrodynamics".
51 : *
52 : * #### Relationship to standard spherical harmonics
53 : *
54 : * The standard expansion of \f$g(\theta, \phi)\f$ in terms of scalar
55 : * spherical harmonics is
56 : * \f{align}
57 : * g(\theta, \phi)
58 : * &=
59 : * \sum_{l=0}^{l_{\max}}\sum_{m=-\min(l, m_{\max})}^{\min(l, m_{\max})}
60 : * A_{lm} Y_{lm}(\theta,\phi),
61 : * \f}
62 : * where \f$Y_{lm}(\theta,\phi)\f$ are the usual complex-valued scalar
63 : * spherical harmonics (as defined, for example, in
64 : * Jackson's "Classical Electrodynamics")
65 : * and \f$A_{lm}\f$ are complex coefficients.
66 : *
67 : * The relationship between the complex coefficients \f$A_{lm}\f$ and
68 : * SPHEREPACK's real-valued \f$a_{lm}\f$ and \f$b_{lm}\f$ is
69 : * \f{align}
70 : * a_{l0} & = \sqrt{\frac{2}{\pi}}A_{l0}&\qquad l\geq 0,\\
71 : * a_{lm} & = (-1)^m\sqrt{\frac{2}{\pi}} \mathrm{Re}(A_{lm})
72 : * &\qquad l\geq 1, m\geq 1, \\
73 : * b_{lm} & = (-1)^m\sqrt{\frac{2}{\pi}} \mathrm{Im}(A_{lm})
74 : * &\qquad l\geq 1, m\geq 1.
75 : * \f}
76 : *
77 : * \note If \f$g\f$ is real,
78 : * \f$A_{lm} = (-1)^m A^\star_{l -m}\f$ (where \f${}^\star\f$ means
79 : * a complex conjugate); this is why we don't need to consider \f$m<0\f$
80 : * in the previous formulas or in SPHEREPACK's expansion.
81 : *
82 : * #### Relationship to real-valued spherical harmonics
83 : *
84 : * Sometimes it is useful to expand a real-valued function in the form
85 : * \f{align}
86 : * g(\theta, \phi)
87 : * &= \sum_{l=0}^\infty\sum_{m=0}^l
88 : * \left[
89 : * c_{lm}\mathrm{Re}(Y_{lm}(\theta, \phi))+
90 : * d_{nm}\mathrm{Im}(Y_{lm}(\theta, \phi))
91 : * \right].
92 : * \f}
93 : * The coefficients here are therefore
94 : * \f{align}
95 : * c_{l0} &= A_{l0},\\
96 : * c_{lm} &= 2\mathrm{Re}(A_{lm}) \qquad m\geq 1,\\
97 : * d_{lm} &=-2\mathrm{Im}(A_{lm}).
98 : * \f}
99 : *
100 : * #### Modal and nodal representations
101 : *
102 : * Internally, SPHEREPACK can represent its expansion in two ways which we
103 : * will refer to as modal and nodal representations:
104 : *
105 : * -# modal: The spectral coefficient arrays \f$a_{lm}\f$ and \f$b_{lm}\f$,
106 : * referred to as `spectral_coefs` in the methods below. For this C++ interface,
107 : * they are saved in a single `DataVector`. To help you index the coefficients
108 : * as expected by this interface, use the class `SpherepackIterator`.
109 : *
110 : * -# nodal: The values at certain collocation points, referred to as
111 : * `collocation_values` in the methods below. This is an array of the expanded
112 : * function \f$g(\theta,\phi)\f$ evaluated at collocation values
113 : * \f$(\theta_i,\phi_j)\f$, where \f$\theta_i\f$ are Gauss-Legendre quadrature
114 : * nodes in the interval \f$(0, \pi)\f$ with \f$i = 0, ..., l_{\max}\f$, and
115 : * \f$\phi_j\f$ is distributed uniformly in \f$(0, 2\pi)\f$ with \f$i = 0, ...,
116 : * 2m_{\max}\f$. The angles of the collocation points can be computed with the
117 : * method `theta_phi_points`.
118 : *
119 : * To convert between the two representations the methods `spec_to_phys` and
120 : * `phys_to_spec` can be used. For internal calculations SPHEREPACK will usually
121 : * convert to spectral coefficients first, so it is in general more efficient to
122 : * use these directly.
123 : *
124 : * Most methods of SPHEREPACK will compute the requested values of e.g.
125 : * `gradient` or `scalar_laplacian` at the collocation points, effectively
126 : * returning an expansion in nodal form as defined above. To evaluate the
127 : * function at arbitrary angles \f$\theta\f$, \f$\phi\f$, these values have to
128 : * be "interpolated" (i.e. the new expansion evaluated) using `interpolate`.
129 : *
130 : * Spherepack stores two types of quantities:
131 : * 1. storage_, which is filled in the constructor and is always const.
132 : * 2. memory_pool_, which is dynamic and thread_local, and is overwritten
133 : * by various member functions that need temporary storage.
134 : */
135 1 : class Spherepack {
136 : public:
137 : /// Type returned by gradient function.
138 1 : using FirstDeriv = tnsr::i<DataVector, 2, Frame::ElementLogical>;
139 : /// Type returned by second derivative function.
140 1 : using SecondDeriv = tnsr::ij<DataVector, 2, Frame::ElementLogical>;
141 :
142 : /// Struct to hold cached information at a set of target interpolation
143 : /// points.
144 : template <typename T>
145 1 : struct InterpolationInfo {
146 0 : InterpolationInfo(size_t l_max, size_t m_max, const gsl::span<double> pmm,
147 : const std::array<T, 2>& target_points);
148 0 : T cos_theta;
149 : // cos(m*phi)
150 0 : DynamicBuffer<T> cos_m_phi;
151 : // sin(m*phi)
152 0 : DynamicBuffer<T> sin_m_phi;
153 : // pbar_factor[m] = Pbar(m,m)*sin(theta)^m
154 0 : DynamicBuffer<T> pbar_factor;
155 :
156 0 : size_t size() const { return num_points_; }
157 0 : size_t m_max() const { return m_max_; }
158 0 : size_t l_max() const { return l_max_; }
159 :
160 : private:
161 0 : size_t l_max_;
162 0 : size_t m_max_;
163 0 : size_t num_points_;
164 : };
165 :
166 : /// Here l_max and m_max are the largest fully-represented l and m in
167 : /// the Ylm expansion.
168 1 : Spherepack(size_t l_max, size_t m_max);
169 :
170 : /// @{
171 : /// Static functions to return the correct sizes of vectors of
172 : /// collocation points and spectral coefficients for a given l_max
173 : /// and m_max. Useful for allocating space without having to create
174 : /// a Spherepack.
175 1 : SPECTRE_ALWAYS_INLINE static constexpr size_t physical_size(
176 : const size_t l_max, const size_t m_max) {
177 : return (l_max + 1) * (2 * m_max + 1);
178 : }
179 : /// \note `spectral_size` is the size of the buffer that holds the
180 : /// coefficients; it is not the number of coefficients (which is
181 : /// \f$m_{\rm max}^2+(l_{\rm max}-m_{\rm max})(2m_{\rm max}+1)\f$).
182 : /// To simplify its internal indexing, SPHEREPACK uses a buffer with
183 : /// more space than necessary. See SpherepackIterator for
184 : /// how to index the coefficients in the buffer.
185 1 : SPECTRE_ALWAYS_INLINE static constexpr size_t spectral_size(
186 : const size_t l_max, const size_t m_max) {
187 : return 2 * (l_max + 1) * (m_max + 1);
188 : }
189 : /// @}
190 :
191 : /// @{
192 : /// Sizes in physical and spectral space for this instance.
193 1 : size_t l_max() const { return l_max_; }
194 1 : size_t m_max() const { return m_max_; }
195 1 : size_t physical_size() const { return n_theta_ * n_phi_; }
196 1 : size_t spectral_size() const { return spectral_size_; }
197 : /// @}
198 :
199 0 : std::array<size_t, 2> physical_extents() const {
200 : return {{n_theta_, n_phi_}};
201 : }
202 :
203 : /// @{
204 : /// Collocation points theta and phi.
205 : ///
206 : /// The phi points are uniform in phi, with the first point
207 : /// at phi=0.
208 : ///
209 : /// The theta points are Gauss-Legendre in \f$\cos(\theta)\f$,
210 : /// so there are no points at the poles.
211 1 : SPECTRE_ALWAYS_INLINE const std::vector<double>& theta_points() const {
212 : return storage_.theta;
213 : }
214 1 : SPECTRE_ALWAYS_INLINE const std::vector<double>& phi_points() const {
215 : return storage_.phi;
216 : }
217 1 : std::array<DataVector, 2> theta_phi_points() const;
218 : /// @}
219 :
220 : /// @{
221 : /// Spectral transformations.
222 : /// To act on a slice of the input and output arrays, specify strides
223 : /// and offsets.
224 1 : void phys_to_spec(gsl::not_null<double*> spectral_coefs,
225 : gsl::not_null<const double*> collocation_values,
226 : size_t physical_stride = 1, size_t physical_offset = 0,
227 : size_t spectral_stride = 1,
228 : size_t spectral_offset = 0) const {
229 : phys_to_spec_impl(spectral_coefs, collocation_values, physical_stride,
230 : physical_offset, spectral_stride, spectral_offset, false);
231 : }
232 1 : void spec_to_phys(gsl::not_null<double*> collocation_values,
233 : gsl::not_null<const double*> spectral_coefs,
234 : size_t spectral_stride = 1, size_t spectral_offset = 0,
235 : size_t physical_stride = 1,
236 : size_t physical_offset = 0) const {
237 : spec_to_phys_impl(collocation_values, spectral_coefs, spectral_stride,
238 : spectral_offset, physical_stride, physical_offset, false);
239 : };
240 : /// @}
241 :
242 : /// @{
243 : /// Spectral transformations where `collocation_values` and
244 : /// `spectral_coefs` are assumed to point to 3-dimensional arrays
245 : /// (I1 x S2 topology), and the transformations are done for all
246 : /// 'radial' points at once by internally looping over all values of
247 : /// the offset from zero to `stride`-1 (the physical and spectral
248 : /// strides are equal and are called `stride`).
249 1 : void phys_to_spec_all_offsets(gsl::not_null<double*> spectral_coefs,
250 : gsl::not_null<const double*> collocation_values,
251 : size_t stride) const {
252 : phys_to_spec_impl(spectral_coefs, collocation_values, stride, 0, stride, 0,
253 : true);
254 : }
255 1 : void spec_to_phys_all_offsets(gsl::not_null<double*> collocation_values,
256 : gsl::not_null<const double*> spectral_coefs,
257 : size_t stride) const {
258 : spec_to_phys_impl(collocation_values, spectral_coefs, stride, 0, stride, 0,
259 : true);
260 : };
261 : /// @}
262 :
263 : /// @{
264 : /// Simpler, less general interfaces to `phys_to_spec` and `spec_to_phys`.
265 : /// Acts on a slice of the input and returns a unit-stride result.
266 1 : DataVector phys_to_spec(const DataVector& collocation_values,
267 : size_t physical_stride = 1,
268 : size_t physical_offset = 0) const;
269 1 : DataVector spec_to_phys(const DataVector& spectral_coefs,
270 : size_t spectral_stride = 1,
271 : size_t spectral_offset = 0) const;
272 : /// @}
273 :
274 : /// @{
275 : /// Simpler, less general interfaces to `phys_to_spec_all_offsets`
276 : /// and `spec_to_phys_all_offsets`. Result has the same stride as
277 : /// the input.
278 1 : DataVector phys_to_spec_all_offsets(const DataVector& collocation_values,
279 : size_t stride) const;
280 1 : DataVector spec_to_phys_all_offsets(const DataVector& spectral_coefs,
281 : size_t stride) const;
282 : /// @}
283 :
284 : /// Computes Pfaffian derivative (df/dtheta, csc(theta) df/dphi) at
285 : /// the collocation values.
286 : /// To act on a slice of the input and output arrays, specify stride
287 : /// and offset (assumed to be the same for input and output).
288 1 : void gradient(const std::array<double*, 2>& df,
289 : gsl::not_null<const double*> collocation_values,
290 : size_t physical_stride = 1, size_t physical_offset = 0) const;
291 :
292 : /// Same as `gradient`, but takes the spectral coefficients (rather
293 : /// than collocation values) of the function. This is more
294 : /// efficient if one happens to already have the spectral
295 : /// coefficients.
296 : /// To act on a slice of the input and output arrays, specify strides
297 : /// and offsets.
298 1 : void gradient_from_coefs(const std::array<double*, 2>& df,
299 : gsl::not_null<const double*> spectral_coefs,
300 : size_t spectral_stride = 1,
301 : size_t spectral_offset = 0,
302 : size_t physical_stride = 1,
303 : size_t physical_offset = 0) const {
304 : gradient_from_coefs_impl(df, spectral_coefs, spectral_stride,
305 : spectral_offset, physical_stride, physical_offset,
306 : false);
307 : }
308 :
309 : /// @{
310 : /// Same as `gradient` but pointers are assumed to point to
311 : /// 3-dimensional arrays (I1 x S2 topology), and the gradient is
312 : /// done for all 'radial' points at once by internally looping
313 : /// over all values of the offset from zero to `stride`-1.
314 1 : void gradient_all_offsets(const std::array<double*, 2>& df,
315 : gsl::not_null<const double*> collocation_values,
316 : size_t stride = 1) const;
317 :
318 1 : SPECTRE_ALWAYS_INLINE void gradient_from_coefs_all_offsets(
319 : const std::array<double*, 2>& df,
320 : gsl::not_null<const double*> spectral_coefs, size_t stride = 1) const {
321 : gradient_from_coefs_impl(df, spectral_coefs, stride, 0, stride, 0, true);
322 : }
323 : /// @}
324 :
325 : /// @{
326 : /// Simpler, less general interfaces to `gradient`.
327 : /// Acts on a slice of the input and returns a unit-stride result.
328 1 : FirstDeriv gradient(const DataVector& collocation_values,
329 : size_t physical_stride = 1,
330 : size_t physical_offset = 0) const;
331 1 : FirstDeriv gradient_from_coefs(const DataVector& spectral_coefs,
332 : size_t spectral_stride = 1,
333 : size_t spectral_offset = 0) const;
334 : /// @}
335 :
336 : /// @{
337 : /// Simpler, less general interfaces to `gradient_all_offsets`.
338 : /// Result has the same stride as the input.
339 1 : FirstDeriv gradient_all_offsets(const DataVector& collocation_values,
340 : size_t stride = 1) const;
341 1 : FirstDeriv gradient_from_coefs_all_offsets(const DataVector& spectral_coefs,
342 : size_t stride = 1) const;
343 : /// @}
344 :
345 : /// Computes Laplacian in physical space.
346 : /// To act on a slice of the input and output arrays, specify stride
347 : /// and offset (assumed to be the same for input and output).
348 1 : void scalar_laplacian(gsl::not_null<double*> scalar_laplacian,
349 : gsl::not_null<const double*> collocation_values,
350 : size_t physical_stride = 1,
351 : size_t physical_offset = 0) const;
352 :
353 : /// Same as `scalar_laplacian` above, but the input is the spectral
354 : /// coefficients (rather than collocation values) of the function.
355 : /// This is more efficient if one happens to already have the
356 : /// spectral coefficients.
357 : /// To act on a slice of the input and output arrays, specify strides
358 : /// and offsets.
359 1 : void scalar_laplacian_from_coefs(gsl::not_null<double*> scalar_laplacian,
360 : gsl::not_null<const double*> spectral_coefs,
361 : size_t spectral_stride = 1,
362 : size_t spectral_offset = 0,
363 : size_t physical_stride = 1,
364 : size_t physical_offset = 0) const;
365 :
366 : /// @{
367 : /// Simpler, less general interfaces to `scalar_laplacian`.
368 : /// Acts on a slice of the input and returns a unit-stride result.
369 1 : DataVector scalar_laplacian(const DataVector& collocation_values,
370 : size_t physical_stride = 1,
371 : size_t physical_offset = 0) const;
372 1 : DataVector scalar_laplacian_from_coefs(const DataVector& spectral_coefs,
373 : size_t spectral_stride = 1,
374 : size_t spectral_offset = 0) const;
375 : /// @}
376 :
377 : /// Computes Pfaffian first and second derivative in physical space.
378 : /// The first derivative is \f$df(i) = d_i f\f$, and the
379 : /// second derivative is \f$ddf(i,j) = d_i (d_j f)\f$,
380 : /// where \f$d_0 = d/d\theta\f$ and \f$d_1 = csc(\theta) d/d\phi\f$.
381 : /// ddf is not symmetric.
382 : /// To act on a slice of the input and output arrays, specify stride
383 : /// and offset (assumed to be the same for input and output).
384 1 : void second_derivative(const std::array<double*, 2>& df,
385 : gsl::not_null<SecondDeriv*> ddf,
386 : gsl::not_null<const double*> collocation_values,
387 : size_t physical_stride = 1,
388 : size_t physical_offset = 0) const;
389 :
390 : /// Simpler, less general interface to second_derivative
391 1 : std::pair<FirstDeriv, SecondDeriv> first_and_second_derivative(
392 : const DataVector& collocation_values) const;
393 :
394 : /// Computes the integral over the sphere.
395 1 : SPECTRE_ALWAYS_INLINE double definite_integral(
396 : gsl::not_null<const double*> collocation_values,
397 : size_t physical_stride = 1, size_t physical_offset = 0) const {
398 : // clang-tidy: 'do not use pointer arithmetic'
399 : return ddot_(n_theta_ * n_phi_, storage_.quadrature_weights.data(), 1,
400 : collocation_values.get() + physical_offset, // NOLINT
401 : physical_stride);
402 : }
403 :
404 : /// Returns weights \f$w_i\f$ such that \f$sum_i (c_i w_i)\f$
405 : /// is the definite integral, where \f$c_i\f$ are collocation values
406 : /// at point i.
407 1 : SPECTRE_ALWAYS_INLINE const std::vector<double>& integration_weights() const {
408 : return storage_.quadrature_weights;
409 : }
410 :
411 : /// Adds a constant (i.e. \f$f(\theta,\phi)\f$ += \f$c\f$) to the function
412 : /// given by the spectral coefficients, by modifying the coefficients.
413 1 : SPECTRE_ALWAYS_INLINE static void add_constant(
414 : const gsl::not_null<DataVector*> spectral_coefs, const double c) {
415 : // The factor of sqrt(8) is because of the normalization of
416 : // SPHEREPACK's coefficients.
417 : (*spectral_coefs)[0] += sqrt(8.0) * c;
418 : }
419 :
420 : /// Returns the average of \f$f(\theta,\phi)\f$ over \f$(\theta,\phi)\f$.
421 1 : SPECTRE_ALWAYS_INLINE static double average(
422 : const DataVector& spectral_coefs) {
423 : // The factor of sqrt(8) is because of the normalization of
424 : // SPHEREPACK's coefficients. All other coefficients average to zero.
425 : return spectral_coefs[0] / sqrt(8.0);
426 : }
427 :
428 : /// Sets up the `InterpolationInfo` structure for interpolating onto
429 : /// a set of target \f$(\theta,\phi)\f$ points. Does not depend on
430 : /// the function being interpolated.
431 : template <typename T>
432 1 : InterpolationInfo<T> set_up_interpolation_info(
433 : const std::array<T, 2>& target_points) const;
434 :
435 : /// Interpolates from `collocation_values` onto the points that have
436 : /// been passed into the `set_up_interpolation_info` function.
437 : /// To interpolate a different function on the same spectral grid, there
438 : /// is no need to recompute `interpolation_info`.
439 : /// If you specify stride and offset, acts on a slice of the input values.
440 : /// The output has unit stride.
441 : template <typename T>
442 1 : void interpolate(gsl::not_null<T*> result,
443 : gsl::not_null<const double*> collocation_values,
444 : const InterpolationInfo<T>& interpolation_info,
445 : size_t physical_stride = 1,
446 : size_t physical_offset = 0) const;
447 :
448 : /// Same as `interpolate`, but assumes you have spectral coefficients.
449 : /// This is more efficient if you already have the spectral coefficients
450 : /// available.
451 : /// If you specify stride and offset, acts on a slice of the input coefs.
452 : /// The output has unit stride.
453 : template <typename T, typename R>
454 1 : void interpolate_from_coefs(gsl::not_null<T*> result, const R& spectral_coefs,
455 : const InterpolationInfo<T>& interpolation_info,
456 : size_t spectral_stride = 1,
457 : size_t spectral_offset = 0) const;
458 :
459 : /// Simpler interface to `interpolate`. If you need to call this
460 : /// repeatedly on different `spectral_coefs` or `collocation_values`
461 : /// for the same target points, this is inefficient; instead use
462 : /// `set_up_interpolation_info` and the functions that use
463 : /// `InterpolationInfo`.
464 : template <typename T>
465 1 : T interpolate(const DataVector& collocation_values,
466 : const std::array<T, 2>& target_points) const;
467 : template <typename T>
468 0 : T interpolate_from_coefs(const DataVector& spectral_coefs,
469 : const std::array<T, 2>& target_points) const;
470 :
471 : /// Takes spectral coefficients compatible with `*this`, and either
472 : /// prolongs them or restricts them to be compatible with `target`.
473 : /// This is done by truncation (restriction) or padding with zeros
474 : /// (prolongation).
475 1 : DataVector prolong_or_restrict(const DataVector& spectral_coefs,
476 : const Spherepack& target) const;
477 :
478 : private:
479 : // Spectral transformations and gradient.
480 : // If `loop_over_offset` is true, then `collocation_values` and
481 : // `spectral_coefs` are assumed to point to 3-dimensional
482 : // arrays (I1 x S2 topology), and the transformations are done for
483 : // all 'radial' points at once by looping over all values of the
484 : // offset from zero to stride-1. If `loop_over_offset` is true,
485 : // `physical_stride` must equal `spectral_stride`.
486 0 : void phys_to_spec_impl(gsl::not_null<double*> spectral_coefs,
487 : gsl::not_null<const double*> collocation_values,
488 : size_t physical_stride = 1, size_t physical_offset = 0,
489 : size_t spectral_stride = 1, size_t spectral_offset = 0,
490 : bool loop_over_offset = false) const;
491 0 : void spec_to_phys_impl(gsl::not_null<double*> collocation_values,
492 : gsl::not_null<const double*> spectral_coefs,
493 : size_t spectral_stride = 1, size_t spectral_offset = 0,
494 : size_t physical_stride = 1, size_t physical_offset = 0,
495 : bool loop_over_offset = false) const;
496 0 : void gradient_from_coefs_impl(const std::array<double*, 2>& df,
497 : gsl::not_null<const double*> spectral_coefs,
498 : size_t spectral_stride = 1,
499 : size_t spectral_offset = 0,
500 : size_t physical_stride = 1,
501 : size_t physical_offset = 0,
502 : bool loop_over_offset = false) const;
503 0 : void calculate_collocation_points();
504 0 : void calculate_interpolation_data();
505 0 : void fill_scalar_work_arrays();
506 0 : void fill_vector_work_arrays();
507 0 : size_t l_max_, m_max_, n_theta_, n_phi_;
508 0 : size_t spectral_size_;
509 : // memory_pool_ will be shared by multiple instances of
510 : // Spherepack on the same thread. Because these instances are on
511 : // the same thread, member functions of two or more of these
512 : // instances cannot be called simultaneously. Note that member
513 : // functions do not make any assumptions about the contents of
514 : // memory_pool_ on entry, so between calls to member functions it is
515 : // safe to resize objects in memory_pool_ or to overwrite them with
516 : // arbitrary data.
517 0 : static thread_local Spherepack_detail::MemoryPool memory_pool_;
518 0 : Spherepack_detail::ConstStorage storage_;
519 : }; // class Spherepack
520 :
521 0 : bool operator==(const Spherepack& lhs, const Spherepack& rhs);
522 0 : bool operator!=(const Spherepack& lhs, const Spherepack& rhs);
523 :
524 : } // namespace ylm
|
