OrientationMapHelpers.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <boost/range/combine.hpp>
7 #include <boost/tuple/tuple.hpp>
8 #include <cstddef>
9 #include <vector>
10 
11 #include "ErrorHandling/Assert.hpp"
12 #include "Utilities/Gsl.hpp"
13 #include "Utilities/TMPL.hpp"
14 
15 /// \cond
16 template <size_t VolumeDim>
17 class OrientationMap;
18 template <size_t>
19 class Index;
20 template <typename TagsList>
21 class Variables;
22 
23 // clang-tidy: redundant declarations
24 template <typename Tag, typename TagList>
25 constexpr typename Tag::type& get(Variables<TagList>& v) noexcept; // NOLINT
26 template <typename Tag, typename TagList>
27 constexpr const typename Tag::type& get( // NOLINT
28  const Variables<TagList>& v) noexcept;
29 /// \endcond
30 
32 
33 template <typename TagsList>
34 void orient_each_component(
35  const gsl::not_null<Variables<TagsList>*> oriented_variables,
36  const Variables<TagsList>& variables,
37  const std::vector<size_t>& oriented_offset) noexcept {
38  using VectorType = typename Variables<TagsList>::vector_type;
39  tmpl::for_each<TagsList>(
40  [&oriented_variables, &variables, &oriented_offset](auto tag) {
41  using Tag = tmpl::type_from<decltype(tag)>;
42  auto& oriented_tensor = get<Tag>(*oriented_variables);
43  const auto& tensor = get<Tag>(variables);
44  for (decltype(auto) oriented_and_tensor_components :
45  boost::combine(oriented_tensor, tensor)) {
46  VectorType& oriented_tensor_component =
47  boost::get<0>(oriented_and_tensor_components);
48  const VectorType& tensor_component =
49  boost::get<1>(oriented_and_tensor_components);
50  for (size_t s = 0; s < tensor_component.size(); ++s) {
51  oriented_tensor_component[oriented_offset[s]] = tensor_component[s];
52  }
53  }
54  });
55 }
56 
57 template <size_t VolumeDim>
58 std::vector<size_t> oriented_offset(
59  const Index<VolumeDim>& extents,
60  const OrientationMap<VolumeDim>& orientation_of_neighbor) noexcept;
61 
62 inline std::vector<size_t> oriented_offset_on_slice(
63  const Index<0>& /*slice_extents*/, const size_t /*sliced_dim*/,
64  const OrientationMap<1>& /*orientation_of_neighbor*/) noexcept {
65  // There is only one point on a slice of a 1D mesh
66  return {0};
67 }
68 
69 std::vector<size_t> oriented_offset_on_slice(
70  const Index<1>& slice_extents, size_t sliced_dim,
71  const OrientationMap<2>& orientation_of_neighbor) noexcept;
72 
73 std::vector<size_t> oriented_offset_on_slice(
74  const Index<2>& slice_extents, size_t sliced_dim,
75  const OrientationMap<3>& orientation_of_neighbor) noexcept;
76 
77 } // namespace OrientationMapHelpers_detail
78 
79 /// \ingroup ComputationalDomainGroup
80 /// Orient variables to the data-storage order of a neighbor element with
81 /// the given orientation.
82 /// @{
83 template <size_t VolumeDim, typename TagsList>
84 Variables<TagsList> orient_variables(
85  const Variables<TagsList>& variables, const Index<VolumeDim>& extents,
86  const OrientationMap<VolumeDim>& orientation_of_neighbor) noexcept {
87  // Skip work (aside from a copy) if neighbor is aligned
88  if (orientation_of_neighbor.is_aligned()) {
89  return variables;
90  }
91 
92  const size_t number_of_grid_points = extents.product();
93  ASSERT(variables.number_of_grid_points() == number_of_grid_points,
94  "Inconsistent `variables` and `extents`:\n"
95  " variables.number_of_grid_points() = "
96  << variables.number_of_grid_points()
97  << "\n"
98  " extents.product() = "
99  << extents.product());
100  Variables<TagsList> oriented_variables(number_of_grid_points);
101  const auto oriented_offset = OrientationMapHelpers_detail::oriented_offset(
102  extents, orientation_of_neighbor);
103  OrientationMapHelpers_detail::orient_each_component(
104  make_not_null(&oriented_variables), variables, oriented_offset);
105 
106  return oriented_variables;
107 }
108 
109 template <size_t VolumeDim, typename TagsList>
110 Variables<TagsList> orient_variables_on_slice(
111  const Variables<TagsList>& variables_on_slice,
112  const Index<VolumeDim - 1>& slice_extents, const size_t sliced_dim,
113  const OrientationMap<VolumeDim>& orientation_of_neighbor) noexcept {
114  // Skip work (aside from a copy) if neighbor slice is aligned
115  if (orientation_of_neighbor.is_aligned()) {
116  return variables_on_slice;
117  }
118 
119  const size_t number_of_grid_points = slice_extents.product();
120  ASSERT(variables_on_slice.number_of_grid_points() == number_of_grid_points,
121  "Inconsistent `variables_on_slice` and `slice_extents`:\n"
122  " variables_on_slice.number_of_grid_points() = "
123  << variables_on_slice.number_of_grid_points()
124  << "\n"
125  " slice_extents.product() = "
126  << slice_extents.product());
127  Variables<TagsList> oriented_variables(number_of_grid_points);
128  const auto oriented_offset =
129  OrientationMapHelpers_detail::oriented_offset_on_slice(
130  slice_extents, sliced_dim, orientation_of_neighbor);
131  OrientationMapHelpers_detail::orient_each_component(
132  make_not_null(&oriented_variables), variables_on_slice, oriented_offset);
133 
134  return oriented_variables;
135 }
136 /// }@
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: OrientationMapHelpers.hpp:84
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
Definition: OrientationMapHelpers.cpp:203
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: OrientationMapHelpers.hpp:110
Defines macro ASSERT.
Wraps the template metaprogramming library used (brigand)
An integer multi-index.
Definition: Index.hpp:30
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
Require a pointer to not be a nullptr
Definition: ConservativeFromPrimitive.hpp:12