Projection.hpp
1 // 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 #include <ostream>
10 
11 /// \cond
12 class Matrix;
13 template <size_t Dim>
14 class Mesh;
15 /// \endcond
16 
17 namespace Spectral {
18 
19 /// The portion of a mesh covered by a child mesh.
20 enum class ChildSize { Full, UpperHalf, LowerHalf };
21 
22 /// The portion of an element covered by a mortar.
24 
25 std::ostream& operator<<(std::ostream& os, ChildSize mortar_size) noexcept;
26 
27 /// Determine whether data needs to be projected between a child mesh and its
28 /// parent mesh. If no projection is necessary the data may be used as-is.
29 /// Projection is necessary if the child is either p-refined or h-refined
30 /// relative to its parent, or both. This operation is symmetric, i.e. it is
31 /// irrelevant in which order the child and the parent mesh are passed in.
32 template <size_t Dim>
33 bool needs_projection(const Mesh<Dim>& mesh1, const Mesh<Dim>& mesh2,
34  const std::array<ChildSize, Dim>& child_sizes) noexcept;
35 
36 /*!
37  * \brief The projection matrix from a child mesh to its parent.
38  *
39  * The projection matrices returned by this function (and by
40  * projection_matrix_parent_to_child()) define orthogonal projection operators
41  * between the spaces of functions on a parent mesh and its children. These
42  * projections are usually the correct way to transfer data between meshes in
43  * a mesh-refinement hierarchy, as well as between an element face and its
44  * adjacent mortars.
45  *
46  * These functions assume that the `child_mesh` is at least as fine as the
47  * `parent_mesh`, i.e. functions on the `parent_mesh` can be represented exactly
48  * on the `child_mesh`. In practice this means that functions can be projected
49  * to a mortar (the `child_mesh`) from both adjacent element faces (the
50  * `parent_mesh`) without losing accuracy. Similarly, functions in a
51  * mesh-refinement hierarchy don't lose accuracy when an element is split
52  * (h-refined). For this reason, the `projection_matrix_child_to_parent` is
53  * sometimes referred to as a "restriction operator" and the
54  * `projection_matrix_parent_to_child` as a "prolongation operator".
55  *
56  * \par Massive quantities
57  * If the quantity that should be projected is not a function over the
58  * computational grid but a "massive" residual, i.e. a quantity
59  * \f$\int_{\Omega_k} f(x) \psi_p(x) \mathrm{d}V\f$ where \f$\psi_p\f$ are the
60  * basis functions on the mesh, then pass `true` for the parameter
61  * `operand_is_massive` (default is `false`). The restriction operator for this
62  * case is just the transpose of the prolongation operator, i.e. just an
63  * interpolation matrix transpose. Note that the "massive" residual already
64  * takes the difference in element size between parent and children into account
65  * by including a Jacobian in the volume element of the integral.
66  *
67  * \par Implementation details
68  * The half-interval projections are based on an equation derived by
69  * Saul. This shows that the projection from the spectral basis for
70  * the entire interval to the spectral basis for the upper half
71  * interval is
72  * \f{equation*}
73  * T_{jk} = \frac{2 j + 1}{2} 2^j \sum_{n=0}^{j-k} \binom{j}{k+n}
74  * \binom{(j + k + n - 1)/2}{j} \frac{(k + n)!^2}{(2 k + n + 1)! n!}
75  * \f}
76  */
78  const Mesh<1>& child_mesh, const Mesh<1>& parent_mesh, ChildSize size,
79  bool operand_is_massive = false) noexcept;
80 
81 /// The projection matrix from a child mesh to its parent, in `Dim` dimensions.
82 template <size_t Dim>
83 std::array<std::reference_wrapper<const Matrix>, Dim>
84 projection_matrix_child_to_parent(const Mesh<Dim>& child_mesh,
85  const Mesh<Dim>& parent_mesh,
86  const std::array<ChildSize, Dim>& child_sizes,
87  bool operand_is_massive = false) noexcept;
88 
89 /// The projection matrix from a parent mesh to one of its children.
90 ///
91 /// \see projection_matrix_child_to_parent()
92 const Matrix& projection_matrix_parent_to_child(const Mesh<1>& parent_mesh,
93  const Mesh<1>& child_mesh,
94  ChildSize size) noexcept;
95 
96 /// The projection matrix from a parent mesh to one of its children, in `Dim`
97 /// dimensions
98 template <size_t Dim>
99 std::array<std::reference_wrapper<const Matrix>, Dim>
101  const Mesh<Dim>& parent_mesh, const Mesh<Dim>& child_mesh,
102  const std::array<ChildSize, Dim>& child_sizes) noexcept;
103 
104 } // namespace Spectral
Matrix
A dynamically sized matrix of doubles with column-major storage.
Definition: Matrix.hpp:19
Spectral::needs_projection
bool needs_projection(const Mesh< Dim > &mesh1, const Mesh< Dim > &mesh2, const std::array< ChildSize, Dim > &child_sizes) noexcept
Determine whether data needs to be projected between a child mesh and its parent mesh....
Spectral::projection_matrix_child_to_parent
const Matrix & projection_matrix_child_to_parent(const Mesh< 1 > &child_mesh, const Mesh< 1 > &parent_mesh, ChildSize size, bool operand_is_massive=false) noexcept
The projection matrix from a child mesh to its parent.
functional
std::ostream
cstddef
Spectral::projection_matrix_parent_to_child
const Matrix & projection_matrix_parent_to_child(const Mesh< 1 > &parent_mesh, const Mesh< 1 > &child_mesh, ChildSize size) noexcept
The projection matrix from a parent mesh to one of its children.
array
Spectral::ChildSize
ChildSize
The portion of a mesh covered by a child mesh.
Definition: Projection.hpp:20
Mesh
Holds the number of grid points, basis, and quadrature in each direction of the computational grid.
Definition: Mesh.hpp:49
ostream
dg::mortar_size
MortarSize< Dim - 1 > mortar_size(const ElementId< Dim > &self, const ElementId< Dim > &neighbor, size_t dimension, const OrientationMap< Dim > &orientation) noexcept
Spectral
Functionality associated with a particular choice of basis functions and quadrature for spectral oper...
Definition: ComplexDataView.hpp:13