VariablesHelpers.hpp
Go to the documentation of this file.
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 /// \file
5 /// Defines helper functions for use with Variables class
6 
7 #pragma once
8 
9 #include <boost/range/combine.hpp>
10 #include <boost/tuple/tuple.hpp>
11 #include <cstddef>
12 #include <ostream>
13 #include <vector>
14 
15 #include "DataStructures/DataVector.hpp"
16 #include "DataStructures/SliceIterator.hpp"
17 #include "ErrorHandling/Assert.hpp"
18 #include "Utilities/Gsl.hpp"
19 #include "Utilities/TMPL.hpp"
20 #include "Utilities/TypeTraits.hpp"
21 
22 /// \cond
23 template <size_t VolumeDim>
24 class OrientationMap;
25 template <size_t>
26 class Index;
27 template <typename TagsList>
28 class Variables;
29 
30 // clang-tidy: redundant declarations
31 template <typename Tag, typename TagList>
32 constexpr typename Tag::type& get(Variables<TagList>& v) noexcept; // NOLINT
33 template <typename Tag, typename TagList>
34 constexpr const typename Tag::type& get( // NOLINT
35  const Variables<TagList>& v) noexcept;
36 /// \endcond
37 
38 // @{
39 /*!
40  * \ingroup DataStructuresGroup
41  * \brief Slices the data within `vars` to a codimension 1 slice. The
42  * slice has a constant logical coordinate in direction `sliced_dim`,
43  * slicing the volume at `fixed_index` in that dimension. For
44  * example, to get the lower boundary of `sliced_dim`, pass `0` for
45  * `fixed_index`; to get the upper boundary, pass
46  * `extents[sliced_dim] - 1`.
47  *
48  * \see add_slice_to_data
49  *
50  * returns Variables class sliced to a hypersurface.
51  */
52 template <std::size_t VolumeDim, typename TagsList>
53 void data_on_slice(const gsl::not_null<Variables<TagsList>*> interface_vars,
54  const Variables<TagsList>& vars,
55  const Index<VolumeDim>& element_extents,
56  const size_t sliced_dim, const size_t fixed_index) noexcept {
57  const size_t interface_grid_points =
58  element_extents.slice_away(sliced_dim).product();
59  const size_t volume_grid_points = vars.number_of_grid_points();
60  constexpr const size_t number_of_independent_components =
61  Variables<TagsList>::number_of_independent_components;
62 
63  if (interface_vars->number_of_grid_points() != interface_grid_points) {
64  *interface_vars = Variables<TagsList>(interface_grid_points);
65  }
66  using value_type = typename Variables<TagsList>::value_type;
67  const value_type* vars_data = vars.data();
68  value_type* interface_vars_data = interface_vars->data();
69  for (SliceIterator si(element_extents, sliced_dim, fixed_index); si; ++si) {
70  for (size_t i = 0; i < number_of_independent_components; ++i) {
71  // clang-tidy: do not use pointer arithmetic
72  interface_vars_data[si.slice_offset() + // NOLINT
73  i * interface_grid_points] = // NOLINT
74  vars_data[si.volume_offset() + i * volume_grid_points]; // NOLINT
75  }
76  }
77 }
78 
79 template <std::size_t VolumeDim, typename TagsList>
80 Variables<TagsList> data_on_slice(const Variables<TagsList>& vars,
81  const Index<VolumeDim>& element_extents,
82  const size_t sliced_dim,
83  const size_t fixed_index) noexcept {
84  Variables<TagsList> interface_vars(
85  element_extents.slice_away(sliced_dim).product());
86  data_on_slice(make_not_null(&interface_vars), vars, element_extents,
87  sliced_dim, fixed_index);
88  return interface_vars;
89 }
90 // @}
91 
92 // @{
93 /*!
94  * \ingroup DataStructuresGroup
95  * \brief Slices volume `Tensor`s into a `Variables`
96  *
97  * The slice has a constant logical coordinate in direction `sliced_dim`,
98  * slicing the volume at `fixed_index` in that dimension. For
99  * example, to get the lower boundary of `sliced_dim`, pass `0` for
100  * `fixed_index`; to get the upper boundary, pass
101  * `extents[sliced_dim] - 1`.
102  */
103 template <typename... TagsToSlice, size_t VolumeDim>
105  const gsl::not_null<Variables<tmpl::list<TagsToSlice...>>*> interface_vars,
106  const Index<VolumeDim>& element_extents, const size_t sliced_dim,
107  const size_t fixed_index,
108  const typename TagsToSlice::type&... tensors) noexcept {
109  const size_t interface_grid_points =
110  element_extents.slice_away(sliced_dim).product();
111  if (interface_vars->number_of_grid_points() != interface_grid_points) {
112  *interface_vars =
113  Variables<tmpl::list<TagsToSlice...>>(interface_grid_points);
114  }
115  for (SliceIterator si(element_extents, sliced_dim, fixed_index); si; ++si) {
116  const auto lambda = [&si](auto& interface_tensor,
117  const auto& volume_tensor) noexcept {
118  for (decltype(auto) interface_and_volume_tensor_components :
119  boost::combine(interface_tensor, volume_tensor)) {
120  boost::get<0>(
121  interface_and_volume_tensor_components)[si.slice_offset()] =
122  boost::get<1>(
123  interface_and_volume_tensor_components)[si.volume_offset()];
124  }
125  };
126  expand_pack((lambda(get<TagsToSlice>(*interface_vars), tensors),
127  cpp17::void_type{})...);
128  }
129 }
130 
131 template <typename... TagsToSlice, size_t VolumeDim>
132 Variables<tmpl::list<TagsToSlice...>> data_on_slice(
133  const Index<VolumeDim>& element_extents, const size_t sliced_dim,
134  const size_t fixed_index,
135  const typename TagsToSlice::type&... tensors) noexcept {
136  Variables<tmpl::list<TagsToSlice...>> interface_vars(
137  element_extents.slice_away(sliced_dim).product());
138  data_on_slice<TagsToSlice...>(make_not_null(&interface_vars), element_extents,
139  sliced_dim, fixed_index, tensors...);
140  return interface_vars;
141 }
142 // @}
143 
144 /*!
145  * \ingroup DataStructuresGroup
146  * \brief Adds data on a codimension 1 slice to a volume quantity. The
147  * slice has a constant logical coordinate in direction `sliced_dim`,
148  * slicing the volume at `fixed_index` in that dimension. For
149  * example, to add to the lower boundary of `sliced_dim`, pass `0` for
150  * `fixed_index`; to add to the upper boundary, pass
151  * `extents[sliced_dim] - 1`.
152  *
153  * \see data_on_slice
154  */
155 template <std::size_t VolumeDim, typename TagsList>
156 void add_slice_to_data(const gsl::not_null<Variables<TagsList>*> volume_vars,
157  const Variables<TagsList>& vars_on_slice,
158  const Index<VolumeDim>& extents, const size_t sliced_dim,
159  const size_t fixed_index) noexcept {
160  constexpr const size_t number_of_independent_components =
161  Variables<TagsList>::number_of_independent_components;
162  const size_t volume_grid_points = extents.product();
163  const size_t slice_grid_points = extents.slice_away(sliced_dim).product();
164  ASSERT(volume_vars->number_of_grid_points() == volume_grid_points,
165  "volume_vars has wrong number of grid points. Expected "
166  << volume_grid_points << ", got "
167  << volume_vars->number_of_grid_points());
168  ASSERT(vars_on_slice.number_of_grid_points() == slice_grid_points,
169  "vars_on_slice has wrong number of grid points. Expected "
170  << slice_grid_points << ", got "
171  << vars_on_slice.number_of_grid_points());
172  using value_type = typename Variables<TagsList>::value_type;
173  value_type* const volume_data = volume_vars->data();
174  const value_type* const slice_data = vars_on_slice.data();
175  for (SliceIterator si(extents, sliced_dim, fixed_index); si; ++si) {
176  for (size_t i = 0; i < number_of_independent_components; ++i) {
177  // clang-tidy: do not use pointer arithmetic
178  volume_data[si.volume_offset() + i * volume_grid_points] += // NOLINT
179  slice_data[si.slice_offset() + i * slice_grid_points]; // NOLINT
180  }
181  }
182 }
183 
184 namespace OrientVariables_detail {
185 
186 template <typename TagsList>
187 void orient_each_component(
188  const gsl::not_null<Variables<TagsList>*> oriented_variables,
189  const Variables<TagsList>& variables,
190  const std::vector<size_t>& oriented_offset) noexcept {
191  tmpl::for_each<TagsList>(
192  [&oriented_variables, &variables, &oriented_offset](auto tag) {
193  using Tag = tmpl::type_from<decltype(tag)>;
194  auto& oriented_tensor = get<Tag>(*oriented_variables);
195  const auto& tensor = get<Tag>(variables);
196  for (decltype(auto) oriented_and_tensor_components :
197  boost::combine(oriented_tensor, tensor)) {
198  DataVector& oriented_tensor_component =
199  boost::get<0>(oriented_and_tensor_components);
200  const DataVector& tensor_component =
201  boost::get<1>(oriented_and_tensor_components);
202  for (size_t s = 0; s < tensor_component.size(); ++s) {
203  oriented_tensor_component[oriented_offset[s]] = tensor_component[s];
204  }
205  }
206  });
207 }
208 
209 template <size_t VolumeDim>
210 std::vector<size_t> oriented_offset(
211  const Index<VolumeDim>& extents,
212  const OrientationMap<VolumeDim>& orientation_of_neighbor) noexcept;
213 
214 inline std::vector<size_t> oriented_offset_on_slice(
215  const Index<0>& /*slice_extents*/, const size_t /*sliced_dim*/,
216  const OrientationMap<1>& /*orientation_of_neighbor*/) noexcept {
217  // There is only one point on a slice of a 1D mesh
218  return {0};
219 }
220 
221 std::vector<size_t> oriented_offset_on_slice(
222  const Index<1>& slice_extents, size_t sliced_dim,
223  const OrientationMap<2>& orientation_of_neighbor) noexcept;
224 
225 std::vector<size_t> oriented_offset_on_slice(
226  const Index<2>& slice_extents, size_t sliced_dim,
227  const OrientationMap<3>& orientation_of_neighbor) noexcept;
228 
229 } // namespace OrientVariables_detail
230 
231 /// \ingroup DataStructuresGroup
232 /// Orient variables to the data-storage order of a neighbor element with
233 /// the given orientation.
234 /// @{
235 template <size_t VolumeDim, typename TagsList>
236 Variables<TagsList> orient_variables(
237  const Variables<TagsList>& variables, const Index<VolumeDim>& extents,
238  const OrientationMap<VolumeDim>& orientation_of_neighbor) noexcept {
239  // Skip work (aside from a copy) if neighbor is aligned
240  if (orientation_of_neighbor.is_aligned()) {
241  return variables;
242  }
243 
244  const size_t number_of_grid_points = extents.product();
245  ASSERT(variables.number_of_grid_points() == number_of_grid_points,
246  "Inconsistent `variables` and `extents`:\n"
247  " variables.number_of_grid_points() = "
248  << variables.number_of_grid_points()
249  << "\n"
250  " extents.product() = "
251  << extents.product());
252  Variables<TagsList> oriented_variables(number_of_grid_points);
253  const auto oriented_offset =
254  OrientVariables_detail::oriented_offset(extents, orientation_of_neighbor);
255  OrientVariables_detail::orient_each_component(
256  make_not_null(&oriented_variables), variables, oriented_offset);
257 
258  return oriented_variables;
259 }
260 
261 template <size_t VolumeDim, typename TagsList>
262 Variables<TagsList> orient_variables_on_slice(
263  const Variables<TagsList>& variables_on_slice,
264  const Index<VolumeDim - 1>& slice_extents, const size_t sliced_dim,
265  const OrientationMap<VolumeDim>& orientation_of_neighbor) noexcept {
266  // Skip work (aside from a copy) if neighbor slice is aligned
267  if (orientation_of_neighbor.is_aligned()) {
268  return variables_on_slice;
269  }
270 
271  const size_t number_of_grid_points = slice_extents.product();
272  ASSERT(variables_on_slice.number_of_grid_points() == number_of_grid_points,
273  "Inconsistent `variables_on_slice` and `slice_extents`:\n"
274  " variables_on_slice.number_of_grid_points() = "
275  << variables_on_slice.number_of_grid_points()
276  << "\n"
277  " slice_extents.product() = "
278  << slice_extents.product());
279  Variables<TagsList> oriented_variables(number_of_grid_points);
280  const auto oriented_offset = OrientVariables_detail::oriented_offset_on_slice(
281  slice_extents, sliced_dim, orientation_of_neighbor);
282  OrientVariables_detail::orient_each_component(
283  make_not_null(&oriented_variables), variables_on_slice, oriented_offset);
284 
285  return oriented_variables;
286 }
287 /// }@
Mark a return type as being "void". In C++17 void is a regular type under certain circumstances...
Definition: TypeTraits.hpp:47
Variables< TagsList > orient_variables(const Variables< TagsList > &variables, const Index< VolumeDim > &extents, const OrientationMap< VolumeDim > &orientation_of_neighbor) noexcept
Orient variables to the data-storage order of a neighbor element with the given orientation.
Definition: VariablesHelpers.hpp:236
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
void data_on_slice(const gsl::not_null< Variables< TagsList > *> interface_vars, const Variables< TagsList > &vars, const Index< VolumeDim > &element_extents, const size_t sliced_dim, const size_t fixed_index) noexcept
Slices the data within vars to a codimension 1 slice. The slice has a constant logical coordinate in ...
Definition: VariablesHelpers.hpp:53
Definition: VariablesHelpers.cpp:202
Defines macro ASSERT.
Variables< TagsList > orient_variables_on_slice(const Variables< TagsList > &variables_on_slice, const Index< VolumeDim - 1 > &slice_extents, const size_t sliced_dim, const OrientationMap< VolumeDim > &orientation_of_neighbor) noexcept
Orient variables to the data-storage order of a neighbor element with the given orientation.
Definition: VariablesHelpers.hpp:262
Stores a collection of function values.
Definition: DataVector.hpp:46
Wraps the template metaprogramming library used (brigand)
An integer multi-index.
Definition: Index.hpp:28
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
Defines type traits, some of which are future STL type_traits header.
void add_slice_to_data(const gsl::not_null< Variables< TagsList > *> volume_vars, const Variables< TagsList > &vars_on_slice, const Index< VolumeDim > &extents, const size_t sliced_dim, const size_t fixed_index) noexcept
Adds data on a codimension 1 slice to a volume quantity. The slice has a constant logical coordinate ...
Definition: VariablesHelpers.hpp:156
Require a pointer to not be a nullptr
Definition: ConservativeFromPrimitive.hpp:12
Iterate over a (dim-1)-dimensional slice.
Definition: SliceIterator.hpp:24
constexpr void expand_pack(Ts &&...) noexcept
Allows zero-cost unordered expansion of a parameter.
Definition: TMPL.hpp:545