MortarHelpers.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <algorithm>
7 #include <array>
8 #include <cstddef>
9 #include <functional>
10 #include <type_traits>
11 #include <utility> // IWYU pragma: keep // for std::forward
12 
17 #include "Domain/Mesh.hpp"
18 #include "ErrorHandling/Assert.hpp"
20 #include "NumericalAlgorithms/LinearOperators/ApplyMatrices.hpp"
21 #include "NumericalAlgorithms/Spectral/Projection.hpp"
22 #include "Utilities/Algorithm.hpp"
23 #include "Utilities/Gsl.hpp"
24 #include "Utilities/MakeArray.hpp"
25 #include "Utilities/TMPL.hpp"
26 /// \cond
27 template <size_t VolumeDim>
28 class ElementId;
29 template <size_t VolumeDim>
30 class OrientationMap;
31 // IWYU pragma: no_forward_declare Variables
32 /// \endcond
33 
34 namespace dg {
35 /// \ingroup DiscontinuousGalerkinGroup
36 /// Find a mesh for a mortar capable of representing data from either
37 /// of two faces.
38 template <size_t Dim>
39 Mesh<Dim> mortar_mesh(const Mesh<Dim>& face_mesh1,
40  const Mesh<Dim>& face_mesh2) noexcept;
41 
42 /// \ingroup DiscontinuousGalerkinGroup
43 /// Determine the size of the mortar (i.e., the part of the face it
44 /// covers) for communicating with a neighbor. This is the size
45 /// relative to the size of \p self, and will not generally agree with
46 /// that determined by \p neighbor.
47 template <size_t Dim>
49  const ElementId<Dim>& self, const ElementId<Dim>& neighbor,
50  size_t dimension, const OrientationMap<Dim>& orientation) noexcept;
51 
52 /// \ingroup DiscontinuousGalerkinGroup
53 /// Project variables from a face to a mortar.
54 template <typename Tags, size_t Dim>
55 Variables<Tags> project_to_mortar(
56  const Variables<Tags>& vars,
57  const Mesh<Dim>& face_mesh,
58  const Mesh<Dim>& mortar_mesh,
60  const Matrix identity{};
61  auto projection_matrices = make_array<Dim>(std::cref(identity));
62 
63  const auto face_slice_meshes = face_mesh.slices();
64  const auto mortar_slice_meshes = mortar_mesh.slices();
65  for (size_t i = 0; i < Dim; ++i) {
66  const auto& face_slice_mesh = gsl::at(face_slice_meshes, i);
67  const auto& mortar_slice_mesh = gsl::at(mortar_slice_meshes, i);
68  const auto& slice_size = gsl::at(mortar_size, i);
69  if (slice_size != Spectral::MortarSize::Full or
70  face_slice_mesh != mortar_slice_mesh) {
71  gsl::at(projection_matrices, i) = projection_matrix_element_to_mortar(
72  slice_size, mortar_slice_mesh, face_slice_mesh);
73  }
74  }
75  return apply_matrices(projection_matrices, vars, face_mesh.extents());
76 }
77 
78 /// \ingroup DiscontinuousGalerkinGroup
79 /// Project variables from a mortar to a face.
80 template <typename Tags, size_t Dim>
81 Variables<Tags> project_from_mortar(
82  const Variables<Tags>& vars,
83  const Mesh<Dim>& face_mesh,
84  const Mesh<Dim>& mortar_mesh,
86  ASSERT(face_mesh != mortar_mesh or
88  [](const Spectral::MortarSize& size) noexcept {
89  return size != Spectral::MortarSize::Full;
90  }),
91  "project_from_mortar should not be called if the interface mesh and "
92  "mortar mesh are identical. Please elide the copy instead.");
93 
94  const Matrix identity{};
95  auto projection_matrices = make_array<Dim>(std::cref(identity));
96 
97  const auto face_slice_meshes = face_mesh.slices();
98  const auto mortar_slice_meshes = mortar_mesh.slices();
99  for (size_t i = 0; i < Dim; ++i) {
100  const auto& face_slice_mesh = gsl::at(face_slice_meshes, i);
101  const auto& mortar_slice_mesh = gsl::at(mortar_slice_meshes, i);
102  const auto& slice_size = gsl::at(mortar_size, i);
103  if (slice_size != Spectral::MortarSize::Full or
104  face_slice_mesh != mortar_slice_mesh) {
105  gsl::at(projection_matrices, i) = projection_matrix_mortar_to_element(
106  slice_size, face_slice_mesh, mortar_slice_mesh);
107  }
108  }
109  return apply_matrices(projection_matrices, vars, mortar_mesh.extents());
110 }
111 
112 namespace MortarHelpers_detail {
113 template <typename NormalDotNumericalFluxComputer,
114  typename... NumericalFluxTags, typename... SelfTags,
115  typename... PackagedTags>
116 void apply_normal_dot_numerical_flux(
117  const gsl::not_null<Variables<tmpl::list<NumericalFluxTags...>>*>
118  numerical_fluxes,
119  const NormalDotNumericalFluxComputer& normal_dot_numerical_flux_computer,
120  const Variables<tmpl::list<SelfTags...>>& self_packaged_data,
121  const Variables<tmpl::list<PackagedTags...>>&
122  neighbor_packaged_data) noexcept {
123  normal_dot_numerical_flux_computer(
124  make_not_null(&get<NumericalFluxTags>(*numerical_fluxes))...,
125  get<PackagedTags>(self_packaged_data)...,
126  get<PackagedTags>(neighbor_packaged_data)...);
127 }
128 } // namespace MortarHelpers_detail
129 
130 /// \ingroup DiscontinuousGalerkinGroup
131 /// Compute the lifted data resulting from computing the numerical flux.
132 ///
133 /// \details
134 /// This applies the numerical flux, projects the result to the face
135 /// mesh if necessary, and then lifts it to the volume (still
136 /// presented only on the face mesh as all other points are zero).
137 ///
138 /// Projection must happen after the numerical flux calculation so
139 /// that the elements on either side of the mortar calculate the same
140 /// result. Projection must happen before flux lifting because we
141 /// want the factor of the magnitude of the unit normal added during
142 /// the lift to cancel the Jacobian factor in integrals to preserve
143 /// conservation; this only happens if the two operations are done on
144 /// the same grid.
145 template <typename FluxCommTypes, typename NormalDotNumericalFluxComputer,
146  size_t Dim, typename LocalData>
148  const NormalDotNumericalFluxComputer& normal_dot_numerical_flux_computer,
149  LocalData&& local_data,
150  const typename FluxCommTypes::PackagedData& remote_data,
151  const Mesh<Dim>& face_mesh, const Mesh<Dim>& mortar_mesh,
152  const size_t extent_perpendicular_to_boundary,
154  -> db::item_type<
157  typename FluxCommTypes::LocalData>,
158  "Second argument must be a FluxCommTypes::LocalData");
159  using variables_tag =
162  normal_dot_numerical_fluxes(mortar_mesh.number_of_grid_points(), 0.0);
163  MortarHelpers_detail::apply_normal_dot_numerical_flux(
164  make_not_null(&normal_dot_numerical_fluxes),
165  normal_dot_numerical_flux_computer, local_data.mortar_data, remote_data);
166 
167  tmpl::for_each<db::get_variables_tags_list<variables_tag>>(
168  [&normal_dot_numerical_fluxes, &local_data](const auto tag) noexcept {
169  using Tag = tmpl::type_from<decltype(tag)>;
170  auto& numerical_flux =
171  get<Tags::NormalDotNumericalFlux<Tag>>(normal_dot_numerical_fluxes);
172  const auto& local_flux =
173  get<Tags::NormalDotFlux<Tag>>(local_data.mortar_data);
174  for (size_t i = 0; i < numerical_flux.size(); ++i) {
175  numerical_flux[i] -= local_flux[i];
176  }
177  });
178 
179  const bool refining =
180  face_mesh != mortar_mesh or
181  std::any_of(mortar_size.begin(), mortar_size.end(),
182  [](const Spectral::MortarSize s) noexcept {
183  return s != Spectral::MortarSize::Full;
184  });
185 
186  return dg::lift_flux(
187  refining ? project_from_mortar(normal_dot_numerical_fluxes, face_mesh,
189  : std::move(normal_dot_numerical_fluxes),
190  extent_perpendicular_to_boundary,
191  std::forward<LocalData>(local_data).magnitude_of_face_normal);
192 }
193 } // namespace dg
Defines class Matrix.
decltype(auto) any_of(const Container &c, UnaryPredicate &&unary_predicate)
Convenience wrapper around std::any_of.
Definition: Algorithm.hpp:148
Definition: InitializeElement.hpp:63
Defines function make_array.
An ElementId uniquely labels an Element. It is constructed from the BlockId of the Block to which the...
Definition: ElementId.hpp:36
Holds the number of grid points, basis, and quadrature in each direction of the computational grid...
Definition: Mesh.hpp:49
Mesh< Dim > mortar_mesh(const Mesh< Dim > &face_mesh1, const Mesh< Dim > &face_mesh2) noexcept
Find a mesh for a mortar capable of representing data from either of two faces.
Definition: MortarHelpers.cpp:19
void apply_matrices(const gsl::not_null< Variables< VariableTags > *> result, const std::array< MatrixType, Dim > &matrices, const Variables< VariableTags > &u, const Index< Dim > &extents) noexcept
Multiply by matrices in each dimension.
Definition: ApplyMatrices.hpp:62
Define prefixes for DataBox tags.
Variables< Tags > project_to_mortar(const Variables< Tags > &vars, const Mesh< Dim > &face_mesh, const Mesh< Dim > &mortar_mesh, const std::array< Spectral::MortarSize, Dim > &mortar_size) noexcept
Project variables from a face to a mortar.
Definition: MortarHelpers.hpp:55
auto compute_boundary_flux_contribution(const NormalDotNumericalFluxComputer &normal_dot_numerical_flux_computer, LocalData &&local_data, const typename FluxCommTypes::PackagedData &remote_data, const Mesh< Dim > &face_mesh, const Mesh< Dim > &mortar_mesh, const size_t extent_perpendicular_to_boundary, const std::array< Spectral::MortarSize, Dim > &mortar_size) noexcept -> db::item_type< db::remove_tag_prefix< typename FluxCommTypes::normal_dot_fluxes_tag >>
Compute the lifted data resulting from computing the numerical flux.
Definition: MortarHelpers.hpp:147
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
MortarSize
The portion of an element covered by a mortar.
Definition: Projection.hpp:18
typename DataBox_detail::remove_tag_prefix_impl< Tag >::type remove_tag_prefix
Remove a prefix from Tag, also removing it from the variables tags if the unwrapped tag is a Tags::Va...
Definition: DataBoxTag.hpp:540
constexpr bool is_same_v
Variable template for is_same.
Definition: TypeTraits.hpp:221
Defines function lift_flux.
A dynamically sized matrix of doubles with column-major storage.
Definition: Matrix.hpp:19
Defines class Variables.
Defines the class template Mesh.
std::array< Spectral::MortarSize, Dim - 1 > mortar_size(const ElementId< Dim > &self, const ElementId< Dim > &neighbor, const size_t dimension, const OrientationMap< Dim > &orientation) noexcept
Determine the size of the mortar (i.e., the part of the face it covers) for communicating with a neig...
Definition: MortarHelpers.cpp:36
Variables< Tags > project_from_mortar(const Variables< Tags > &vars, const Mesh< Dim > &face_mesh, const Mesh< Dim > &mortar_mesh, const std::array< Spectral::MortarSize, Dim > &mortar_size) noexcept
Project variables from a mortar to a face.
Definition: MortarHelpers.hpp:81
Defines macro ASSERT.
Wraps the template metaprogramming library used (brigand)
typename DataBox_detail::item_type_impl< TagList, Tag >::type item_type
Get the type that is returned by the Tag. If it is a base tag then a TagList must be passed as a seco...
Definition: DataBoxTag.hpp:410
Defines functions and classes from the GSL.
gsl::not_null< T * > make_not_null(T *ptr) noexcept
Construct a not_null from a pointer. Often this will be done as an implicit conversion, but it may be necessary to perform the conversion explicitly when type deduction is desired.
Definition: Gsl.hpp:863
tnsr::Ij< DataType, Dim, Frame::NoFrame > identity(const DataType &used_for_type) noexcept
returns the Identity matrix
Definition: Identity.hpp:14
Defines classes SimpleTag, PrefixTag, ComputeTag and several functions for retrieving tag info...
auto lift_flux(Variables< tmpl::list< FluxTags... >> flux, const size_t extent_perpendicular_to_boundary, Scalar< DataVector > magnitude_of_face_normal) noexcept -> Variables< tmpl::list< db::remove_tag_prefix< FluxTags >... >>
Lifts the flux contribution from an interface to the volume.
Definition: LiftFlux.hpp:40
Require a pointer to not be a nullptr
Definition: ConservativeFromPrimitive.hpp:12
constexpr T & at(std::array< T, N > &arr, Size index)
Retrieve a entry from a container, with checks in Debug mode that the index being retrieved is valid...
Definition: Gsl.hpp:124