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 <cstddef>
8 : #include <functional>
9 :
10 : #include "DataStructures/DataVector.hpp"
11 : #include "DataStructures/Matrix.hpp"
12 : #include "DataStructures/Variables.hpp"
13 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
14 : #include "NumericalAlgorithms/Spectral/SegmentSize.hpp"
15 : #include "Utilities/ConstantExpressions.hpp"
16 : #include "Utilities/Gsl.hpp"
17 :
18 : namespace Spectral {
19 : /// Determine whether data needs to be projected between a child mesh and its
20 : /// parent mesh. If no projection is necessary the data may be used as-is.
21 : /// Projection is necessary if the child is either p-refined or h-refined
22 : /// relative to its parent, or both. This operation is symmetric, i.e. it is
23 : /// irrelevant in which order the child and the parent mesh are passed in.
24 : template <size_t Dim>
25 1 : bool needs_projection(const Mesh<Dim>& mesh1, const Mesh<Dim>& mesh2,
26 : const std::array<SegmentSize, Dim>& child_sizes);
27 :
28 : /*!
29 : * \brief The projection matrix from a child mesh to its parent.
30 : *
31 : * The projection matrices returned by this function (and by
32 : * projection_matrix_parent_to_child()) define orthogonal projection operators
33 : * between the spaces of functions on a parent mesh and its children. These
34 : * projections are usually the correct way to transfer data between meshes in
35 : * a mesh-refinement hierarchy, as well as between an element face and its
36 : * adjacent mortars.
37 : *
38 : * These functions assume that the `child_mesh` is at least as fine as the
39 : * `parent_mesh`, i.e. functions on the `parent_mesh` can be represented exactly
40 : * on the `child_mesh`. In practice this means that functions can be projected
41 : * to a mortar (the `child_mesh`) from both adjacent element faces (the
42 : * `parent_mesh`) without losing accuracy. Similarly, functions in a
43 : * mesh-refinement hierarchy don't lose accuracy when an element is split
44 : * (h-refined). For this reason, the `projection_matrix_child_to_parent` is
45 : * sometimes referred to as a "restriction operator" and the
46 : * `projection_matrix_parent_to_child` as a "prolongation operator".
47 : *
48 : * \par Massive quantities
49 : * If the quantity that should be projected is not a function over the
50 : * computational grid but a "massive" residual, i.e. a quantity
51 : * \f$\int_{\Omega_k} f(x) \psi_p(x) \mathrm{d}V\f$ where \f$\psi_p\f$ are the
52 : * basis functions on the mesh, then pass `true` for the parameter
53 : * `operand_is_massive` (default is `false`). The restriction operator for this
54 : * case is just the transpose of the prolongation operator, i.e. just an
55 : * interpolation matrix transpose. Note that the "massive" residual already
56 : * takes the difference in element size between parent and children into account
57 : * by including a Jacobian in the volume element of the integral.
58 : *
59 : * \par Implementation details
60 : * The half-interval projections are based on an equation derived by
61 : * Saul. This shows that the projection from the spectral basis for
62 : * the entire interval to the spectral basis for the upper half
63 : * interval is
64 : * \f{equation*}
65 : * T_{jk} = \frac{2 j + 1}{2} 2^j \sum_{n=0}^{j-k} \binom{j}{k+n}
66 : * \binom{(j + k + n - 1)/2}{j} \frac{(k + n)!^2}{(2 k + n + 1)! n!}
67 : * \f}
68 : *
69 : * \note This and the other matrix-returning projection functions in this file
70 : * (projection_matrix_parent_to_child(), projection_matrices()) only support
71 : * tensor-product bases (Legendre, Fourier, ZernikeB2). They cannot be used for
72 : * spherical-shell meshes that use the `SphericalHarmonic` basis, because the
73 : * angular projection couples the \f$\theta\f$ and \f$\phi\f$ directions and so
74 : * cannot be expressed as the per-dimension matrices these functions return. Use
75 : * Spectral::project() instead whenever spherical-shell meshes must be
76 : * supported; it dispatches to these matrices for tensor-product meshes and to a
77 : * Spherepack-based angular projection for spherical-shell meshes.
78 : */
79 1 : const Matrix& projection_matrix_child_to_parent(
80 : const Mesh<1>& child_mesh, const Mesh<1>& parent_mesh, SegmentSize size,
81 : bool operand_is_massive = false);
82 :
83 : /// The projection matrix from a child mesh to its parent, in `Dim` dimensions.
84 : template <size_t Dim>
85 : std::array<std::reference_wrapper<const Matrix>, Dim>
86 1 : projection_matrix_child_to_parent(
87 : const Mesh<Dim>& child_mesh, const Mesh<Dim>& parent_mesh,
88 : const std::array<SegmentSize, Dim>& child_sizes,
89 : bool operand_is_massive = false);
90 :
91 : /// The projection matrix from a parent mesh to one of its children.
92 : ///
93 : /// \see projection_matrix_child_to_parent()
94 1 : const Matrix& projection_matrix_parent_to_child(const Mesh<1>& parent_mesh,
95 : const Mesh<1>& child_mesh,
96 : SegmentSize size);
97 :
98 : /// The projection matrix from a parent mesh to one of its children, in `Dim`
99 : /// dimensions
100 : template <size_t Dim>
101 : std::array<std::reference_wrapper<const Matrix>, Dim>
102 1 : projection_matrix_parent_to_child(
103 : const Mesh<Dim>& parent_mesh, const Mesh<Dim>& child_mesh,
104 : const std::array<SegmentSize, Dim>& child_sizes);
105 :
106 : /// The projection matrices from a source mesh to a target mesh
107 : /// covering given portions of an element
108 : template <size_t Dim>
109 1 : std::array<std::reference_wrapper<const Matrix>, Dim> projection_matrices(
110 : const Mesh<Dim>& source_mesh, const Mesh<Dim>& target_mesh,
111 : const std::array<SegmentSize, Dim>& source_sizes,
112 : const std::array<SegmentSize, Dim>& target_sizes,
113 : bool operand_is_massive = false);
114 :
115 : /// Change the angular resolution (`l_max`, with `m_max == l_max`) of volume
116 : /// data on a spherical shell from `source_data` to `result_data`, for
117 : /// `num_components` components each laid out with the radial dimension varying
118 : /// fastest. Uses Spherepack prolong/restrict. \see Spectral::project() for the
119 : /// higher-level interface to call this.
120 : ///
121 : /// For non-massive operands this is the L2 (Galerkin) projection,
122 : /// $P_\mathrm{L2}$. For massive operands (`operand_is_massive == true`) the
123 : /// restriction is the transpose of the prolongation (interpolation) operator,
124 : /// $I^T$, which is the L2 projection conjugated by the diagonal matrices $W$ of
125 : /// angular quadrature weights:
126 : /// $I^T = W_\mathrm{target} P_\mathrm{L2} W_\mathrm{source}^{-1}$.
127 : /// This means we divide by the source weights before projecting and multiply by
128 : /// the target weights afterwards. See projection_matrix_child_to_parent() for
129 : /// details on massive operands.
130 1 : void project_spherical_harmonics(gsl::not_null<double*> result_data,
131 : const double* source_data,
132 : size_t num_components,
133 : size_t num_radial_points, size_t l_max_source,
134 : size_t l_max_target, bool operand_is_massive);
135 :
136 : /// @{
137 : /*!
138 : * \brief Project volume data from `source_mesh` to `target_mesh`, writing the
139 : * (resized) result into `*result`.
140 : *
141 : * This is the unified entry point for projection regardless of basis. It
142 : * handles both tensor-product meshes and spherical-shell meshes (Legendre
143 : * radial dimension plus `SphericalHarmonic` angular dimensions with
144 : * `m_max == l_max`):
145 : * - For tensor-product meshes it delegates to the per-dimension projection
146 : * matrices above.
147 : * - For spherical-shell meshes it projects the radial dimension with the 1D
148 : * projection matrix and changes the angular `l_max` with a Spherepack
149 : * prolong/restrict. The angular dimensions cannot be h-refined, so their
150 : * segment sizes must be `Full`.
151 : *
152 : * `source_sizes` and `target_sizes` carry h-refinement information (which
153 : * portion of an element the source/target cover); pass `Full` in every
154 : * dimension for pure p-refinement.
155 : *
156 : * \warning `*result` must not point to `source` (the result is resized to the
157 : * target number of grid points).
158 : */
159 : template <typename VectorType, size_t Dim>
160 1 : void project(gsl::not_null<VectorType*> result, const VectorType& source,
161 : const Mesh<Dim>& source_mesh, const Mesh<Dim>& target_mesh,
162 : const std::array<SegmentSize, Dim>& source_sizes,
163 : const std::array<SegmentSize, Dim>& target_sizes,
164 : bool operand_is_massive = false);
165 :
166 : template <size_t Dim, typename TagList>
167 1 : void project(const gsl::not_null<Variables<TagList>*> result,
168 : const Variables<TagList>& source, const Mesh<Dim>& source_mesh,
169 : const Mesh<Dim>& target_mesh,
170 : const std::array<SegmentSize, Dim>& source_sizes,
171 : const std::array<SegmentSize, Dim>& target_sizes,
172 : const bool operand_is_massive = false) {
173 : // Type-erase to the vector implementation with multiple components
174 : using VectorType = typename Variables<TagList>::vector_type;
175 : using ValueType = typename Variables<TagList>::value_type;
176 : result->initialize(target_mesh.number_of_grid_points());
177 : VectorType result_view(result->data(), result->size());
178 : // NOLINTNEXTLINE(cppcoreguidelines-pro-type-const-cast)
179 : const VectorType source_view(const_cast<ValueType*>(source.data()),
180 : source.size());
181 : project(make_not_null(&result_view), source_view, source_mesh, target_mesh,
182 : source_sizes, target_sizes, operand_is_massive);
183 : }
184 :
185 : template <size_t Dim, typename T>
186 1 : T project(const T& source, const Mesh<Dim>& source_mesh,
187 : const Mesh<Dim>& target_mesh,
188 : const std::array<SegmentSize, Dim>& source_sizes,
189 : const std::array<SegmentSize, Dim>& target_sizes,
190 : const bool operand_is_massive = false) {
191 : T result{};
192 : project(make_not_null(&result), source, source_mesh, target_mesh,
193 : source_sizes, target_sizes, operand_is_massive);
194 : return result;
195 : }
196 : /// @}
197 :
198 : /// @{
199 : /// \brief Performs a perfect hash of the mortars into $2^{d-1}$ slots on the
200 : /// range $[0, 2^{d-1})$.
201 : ///
202 : /// This is particularly useful when hashing into statically-sized maps based
203 : /// on the number of dimensions.
204 : template <size_t DimMinusOne>
205 1 : size_t hash(const std::array<Spectral::SegmentSize, DimMinusOne>& mortar_size);
206 :
207 : template <size_t Dim>
208 0 : struct MortarSizeHash {
209 : template <size_t MaxSize>
210 0 : static constexpr bool is_perfect = MaxSize == two_to_the(Dim);
211 :
212 0 : size_t operator()(
213 : const std::array<Spectral::SegmentSize, Dim - 1>& mortar_size);
214 : };
215 : /// @}
216 : } // namespace Spectral
|
