DomainHelpers.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 DomainHelper functions
6 
7 #pragma once
8 
9 #include <array>
10 #include <cstddef>
11 #include <iosfwd>
12 #include <limits>
13 #include <memory>
14 #include <vector>
15 
16 #include "DataStructures/Index.hpp"
18 #include "Domain/Direction.hpp"
19 #include "Domain/Side.hpp"
21 #include "Utilities/Gsl.hpp"
22 #include "Utilities/MakeArray.hpp"
23 
24 /// \cond
25 template <size_t VolumeDim>
26 class BlockNeighbor;
27 template <typename SourceFrame, typename TargetFrame, size_t Dim>
28 class CoordinateMapBase;
29 template <size_t VolumeDim, typename T>
30 class DirectionMap;
31 template <size_t VolumeDim, typename TargetFrame>
32 class Domain;
33 template <size_t VolumeDim>
34 class OrientationMap;
35 class Option;
36 template <typename T>
37 struct create_from_yaml;
38 /// \endcond
39 
40 /// \ingroup ComputationalDomainGroup
41 /// Each member in `PairOfFaces` holds the global corner ids of a block face.
42 /// `PairOfFaces` is used in setting up periodic boundary conditions by
43 /// identifying the two faces with each other.
44 /// \requires The pair of faces must belong to a single block.
45 struct PairOfFaces {
46  std::vector<size_t> first;
47  std::vector<size_t> second;
48 };
49 
50 /// \ingroup ComputationalDomainGroup
51 /// Sets up the BlockNeighbors using the corner numbering scheme
52 /// to deduce the correct neighbors and orientations. Does not set
53 /// up periodic boundary conditions.
54 template <size_t VolumeDim>
56  const std::vector<std::array<size_t, two_to_the(VolumeDim)>>&
57  corners_of_all_blocks,
60  neighbors_of_all_blocks) noexcept;
61 
62 /// \ingroup ComputationalDomainGroup
63 /// Sets up additional BlockNeighbors corresponding to any
64 /// identifications of faces provided by the user. Can be used
65 /// for manually setting up periodic boundary conditions.
66 template <size_t VolumeDim>
68  const std::vector<PairOfFaces>& identifications,
69  const std::vector<std::array<size_t, two_to_the(VolumeDim)>>&
70  corners_of_all_blocks,
73  neighbors_of_all_blocks) noexcept;
74 
75 /// \ingroup ComputationalDomainGroup
76 /// \brief The corners for a rectilinear domain made of n-cubes.
77 ///
78 /// The `domain_extents` argument holds the number of blocks to have
79 /// in each dimension. The blocks all have aligned orientations by
80 /// construction. The `block_indices_to_exclude` argument allows the user
81 /// to selectively exclude blocks from the resulting domain. This allows
82 /// for the creation of non-trivial shapes such as the net for a tesseract.
83 template <size_t VolumeDim>
85  const Index<VolumeDim>& domain_extents,
86  const std::vector<Index<VolumeDim>>& block_indices_to_exclude = {}) noexcept
87  -> std::vector<std::array<size_t, two_to_the(VolumeDim)>>;
88 
89 /// \ingroup ComputationalDomainGroup
90 /// The number of wedges to include in the Shell domain.
91 enum class ShellWedges {
92  /// Use the entire shell
93  All,
94  /// Use only the four equatorial wedges
96  /// Use only the single wedge along -x
98 };
99 
100 /// \ingroup ComputationalDomainGroup
101 /// These are the CoordinateMaps of the Wedge3Ds used in the Sphere, Shell, and
102 /// binary compact object DomainCreators. This function can also be used to
103 /// wrap the Sphere or Shell in a cube made of six Wedge3Ds.
104 /// The argument `x_coord_of_shell_center` specifies a translation of the Shell
105 /// in the x-direction in the TargetFrame. For example, the BBH DomainCreator
106 /// uses this to set the position of each BH.
107 /// When the argument `use_half_wedges` is set to `true`, the wedges in the
108 /// +z,-z,+y,-y directions are cut in half along their xi-axes. The resulting
109 /// ten CoordinateMaps are used for the outermost Blocks of the BBH Domain.
110 /// The argument `aspect_ratio` sets the equatorial compression factor,
111 /// used by the EquatorialCompression maps which get composed with the Wedges.
112 /// This is done if `aspect_ratio` is set to something other than the default
113 /// value of one. When the argument `use_logarithmic_map` is set to `true`,
114 /// the radial gridpoints of the wedge map are set to be spaced logarithmically.
115 /// The `number_of_layers` is used when the user wants to have multiple layers
116 /// of Blocks in the radial direction.
117 template <typename TargetFrame>
118 auto wedge_coordinate_maps(double inner_radius, double outer_radius,
119  double inner_sphericity, double outer_sphericity,
120  bool use_equiangular_map,
121  double x_coord_of_shell_center = 0.0,
122  bool use_half_wedges = false,
123  double aspect_ratio = 1.0,
124  bool use_logarithmic_map = false,
125  ShellWedges which_wedges = ShellWedges::All,
126  size_t number_of_layers = 1) noexcept
127  -> std::vector<
129 
130 /// \ingroup ComputationalDomainGroup
131 /// These are the ten Frustums used in the DomainCreators for binary compact
132 /// objects. The cubes enveloping the two Shells each have a side length of
133 /// `length_inner_cube`. The ten frustums also make up a cube of their own,
134 /// of side length `length_outer_cube`.
135 template <typename TargetFrame>
136 auto frustum_coordinate_maps(double length_inner_cube, double length_outer_cube,
137  bool use_equiangular_map) noexcept
138  -> std::vector<
139  std::unique_ptr<CoordinateMapBase<Frame::Logical, TargetFrame, 3>>>;
140 
141 /// \ingroup ComputationalDomainGroup
142 /// \brief The corners for a domain with radial layers.
143 ///
144 /// Generates the corners for a Domain which is made of one or more layers
145 /// of Blocks fully enveloping an interior volume, e.g. Shell or Sphere. The
146 /// `number_of_layers` specifies how many of these layers of Blocks to have
147 /// in the final domain.
148 /// `include_central_block` is set to `true` in Sphere, where the interior
149 /// volume is filled with a central Block, and `false` in Shell, where the
150 /// interior volume is left empty.
151 /// The `central_block_corners` are used as seed values to generate the corners
152 /// for the surrounding Blocks.
154  size_t number_of_layers, bool include_central_block,
155  const std::array<size_t, 8>& central_block_corners = {{1, 2, 3, 4, 5, 6, 7,
156  8}},
157  ShellWedges which_wedges = ShellWedges::All) noexcept;
158 
159 /// \ingroup ComputationalDomainGroup
160 /// \brief The corners for a domain with biradial layers.
161 ///
162 /// Generates the corners for a BBH-like Domain which is made of one or more
163 /// layers of Blocks fully enveloping two interior volumes. The
164 /// `number_of_radial_layers` gives the number of layers that fully envelop
165 /// each interior volume with six Blocks each. The `number_of_biradial_layers`
166 /// gives the number of layers that fully envelop both volumes at once, using
167 /// ten Blocks per layer as opposed to six. The `central_block_corners_lhs`
168 /// are used as seed values to generate the corners for the surrounding
169 /// Blocks.
171  size_t number_of_radial_layers, size_t number_of_biradial_layers,
172  bool include_central_block_lhs, bool include_central_block_rhs,
173  const std::array<size_t, 8>& central_block_corners_lhs = {
174  {1, 2, 3, 4, 5, 6, 7, 8}}) noexcept;
175 
176 /// \ingroup ComputationalDomainGroup
177 /// \brief Permutes the corner numbers of an n-cube.
178 ///
179 /// Returns the correct ordering of global corner numbers for a rotated block
180 /// in an otherwise aligned edifice of blocks, given the OrientationMap a
181 /// block aligned with the edifice has relative to this one, and given the
182 /// corner numbering the rotated block would have if it were aligned.
183 /// This is useful in creating domains for testing purposes, e.g.
184 /// RotatedIntervals, RotatedRectangles, and RotatedBricks.
185 template <size_t VolumeDim>
187  const OrientationMap<VolumeDim>& orientation,
188  const std::array<size_t, two_to_the(VolumeDim)>&
189  corners_of_aligned) noexcept;
190 
191 /// \ingroup ComputationalDomainGroup
192 /// \brief The CoordinateMaps for a rectilinear domain of n-cubes.
193 ///
194 /// Allows for both Affine and Equiangular maps.
195 template <typename TargetFrame, size_t VolumeDim>
197  const Index<VolumeDim>& domain_extents,
198  const std::array<std::vector<double>, VolumeDim>& block_demarcations,
199  const std::vector<Index<VolumeDim>>& block_indices_to_exclude = {},
200  const std::vector<OrientationMap<VolumeDim>>& orientations_of_all_blocks =
201  {},
202  bool use_equiangular_map = false) noexcept
205 
206 /// \ingroup ComputationalDomainGroup
207 /// \brief Create a rectilinear Domain of multicubes.
208 ///
209 /// \details Useful for constructing domains for testing non-trivially
210 /// connected rectilinear domains made up of cubes. We refer to a domain of
211 /// this type as an edifice. The `domain_extents` provides the size (in the
212 /// number of blocks) of the initial aligned edifice to construct. The
213 /// `block_indices_to_exclude` parameter is used in refining the shape of
214 /// the edifice from a cube to sometime more non-trivial, such as an L-shape
215 /// or the net of a tesseract. The `block_demarcations` and
216 /// `use_equiangular_map` parameters determine the CoordinateMaps to be used.
217 /// `orientations_of_all_blocks` contains the OrientationMap of the edifice
218 /// relative to each block.
219 ///
220 /// The `identifications` parameter is used when identifying the faces of
221 /// blocks in an edifice. This is used to identify the 1D boundaries in the 2D
222 /// net for a 3D cube to construct a domain with topology S2. Note: If the user
223 /// wishes to rotate the blocks as well as manually identify their faces, the
224 /// user must provide the PairOfFaces corresponding to the rotated corners.
225 template <size_t VolumeDim, typename TargetFrame>
227  const Index<VolumeDim>& domain_extents,
228  const std::array<std::vector<double>, VolumeDim>& block_demarcations,
229  const std::vector<Index<VolumeDim>>& block_indices_to_exclude = {},
230  const std::vector<OrientationMap<VolumeDim>>& orientations_of_all_blocks =
231  {},
232  const std::array<bool, VolumeDim>& dimension_is_periodic =
233  make_array<VolumeDim>(false),
234  const std::vector<PairOfFaces>& identifications = {},
235  bool use_equiangular_map = false) noexcept;
236 
237 /// \ingroup ComputationalDomainGroup
238 /// Iterates over the corners of a VolumeDim-dimensional cube.
239 template <size_t VolumeDim>
241  public:
242  VolumeCornerIterator() noexcept { setup_from_local_corner_number(); }
243 
244  explicit VolumeCornerIterator(size_t initial_local_corner_number) noexcept
245  : local_corner_number_(initial_local_corner_number) {
246  setup_from_local_corner_number();
247  }
249  // The block index is also global corner
250  // index of the lowest corner of the block.
251  Index<VolumeDim> block_index,
252  Index<VolumeDim> global_corner_extents) noexcept
253  : global_corner_number_(
254  collapsed_index(block_index, global_corner_extents)),
255  global_corner_index_(block_index),
256  global_corner_extents_(global_corner_extents) {}
257 
258  void operator++() noexcept {
259  ++local_corner_number_;
260  setup_from_local_corner_number();
261  }
262 
263  explicit operator bool() const noexcept {
264  return local_corner_number_ < two_to_the(VolumeDim);
265  }
266 
267  const size_t& local_corner_number() const noexcept {
268  return local_corner_number_;
269  }
270 
271  size_t global_corner_number() const noexcept {
272  std::array<size_t, VolumeDim> new_indices{};
273  for (size_t i = 0; i < VolumeDim; i++) {
274  gsl::at(new_indices, i) =
275  global_corner_index_[i] +
276  (gsl::at(array_sides_, i) == Side::Upper ? 1 : 0);
277  }
278  const Index<VolumeDim> interior_multi_index(new_indices);
279  return collapsed_index(interior_multi_index, global_corner_extents_);
280  }
281 
282  const std::array<Side, VolumeDim>& operator()() const noexcept {
283  return array_sides_;
284  }
285 
286  const std::array<Side, VolumeDim>& operator*() const noexcept {
287  return array_sides_;
288  }
289 
290  const std::array<double, VolumeDim>& coords_of_corner() const noexcept {
291  return coords_of_corner_;
292  }
293 
294  const std::array<Direction<VolumeDim>, VolumeDim>& directions_of_corner()
295  const noexcept {
296  return array_directions_;
297  }
298 
299  void setup_from_local_corner_number() noexcept {
300  for (size_t i = 0; i < VolumeDim; i++) {
301  gsl::at(coords_of_corner_, i) =
302  2.0 * get_nth_bit(local_corner_number_, i) - 1.0;
303  gsl::at(array_sides_, i) =
304  2 * get_nth_bit(local_corner_number_, i) - 1 == 1 ? Side::Upper
305  : Side::Lower;
306  gsl::at(array_directions_, i) =
307  Direction<VolumeDim>(i, gsl::at(array_sides_, i));
308  }
309  }
310 
311  private:
312  size_t local_corner_number_ = 0;
313  size_t global_corner_number_{std::numeric_limits<size_t>::max()};
314  Index<VolumeDim> global_corner_index_{};
315  Index<VolumeDim> global_corner_extents_{};
316  std::array<Side, VolumeDim> array_sides_ = make_array<VolumeDim>(Side::Lower);
317  std::array<Direction<VolumeDim>, VolumeDim> array_directions_{};
318  std::array<double, VolumeDim> coords_of_corner_ = make_array<VolumeDim>(-1.0);
319 };
320 
321 /// \ingroup ComputationalDomainGroup
322 /// Iterates over the 2^(VolumeDim-1) logical corners of the face of a
323 /// VolumeDim-dimensional cube in the given direction.
324 template <size_t VolumeDim>
326  public:
327  explicit FaceCornerIterator(Direction<VolumeDim> direction) noexcept;
328 
329  void operator++() noexcept {
330  face_index_++;
331  do {
332  index_++;
333  } while (get_nth_bit(index_, direction_.dimension()) ==
334  (direction_.side() == Side::Upper ? 0 : 1));
335  for (size_t i = 0; i < VolumeDim; ++i) {
336  corner_[i] = 2 * static_cast<int>(get_nth_bit(index_, i)) - 1;
337  }
338  }
339 
340  explicit operator bool() const noexcept {
341  return face_index_ < two_to_the(VolumeDim - 1);
342  }
343 
344  tnsr::I<double, VolumeDim, Frame::Logical> operator()() const noexcept {
345  return corner_;
346  }
347 
348  tnsr::I<double, VolumeDim, Frame::Logical> operator*() const noexcept {
349  return corner_;
350  }
351 
352  // Returns the value used to construct the logical corner.
353  size_t volume_index() const noexcept { return index_; }
354 
355  // Returns the number of times operator++ has been called.
356  size_t face_index() const noexcept { return face_index_; }
357 
358  private:
359  const Direction<VolumeDim> direction_;
360  size_t index_;
361  size_t face_index_ = 0;
362  tnsr::I<double, VolumeDim, Frame::Logical> corner_;
363 };
364 
365 template <size_t VolumeDim>
367  Direction<VolumeDim> direction) noexcept
368  : direction_(std::move(direction)),
369  index_(direction.side() == Side::Upper
370  ? two_to_the(direction_.dimension())
371  : 0) {
372  for (size_t i = 0; i < VolumeDim; ++i) {
373  corner_[i] = 2 * static_cast<int>(get_nth_bit(index_, i)) - 1;
374  }
375 }
376 
377 std::ostream& operator<<(std::ostream& os,
378  const ShellWedges& which_wedges) noexcept;
379 
380 template <>
382  static ShellWedges create(const Option& options);
383 };
Defines class template Direction.
auto corners_for_rectilinear_domains(const Index< VolumeDim > &domain_extents, const std::vector< Index< VolumeDim >> &block_indices_to_exclude={}) noexcept -> std::vector< std::array< size_t, two_to_the(VolumeDim)>>
The corners for a rectilinear domain made of n-cubes.
Definition: DomainHelpers.cpp:794
Use the entire shell.
auto wedge_coordinate_maps(double inner_radius, double outer_radius, double inner_sphericity, double outer_sphericity, bool use_equiangular_map, double x_coord_of_shell_center=0.0, bool use_half_wedges=false, double aspect_ratio=1.0, bool use_logarithmic_map=false, ShellWedges which_wedges=ShellWedges::All, size_t number_of_layers=1) noexcept -> std::vector< std::unique_ptr< CoordinateMapBase< Frame::Logical, TargetFrame, 3 >>>
These are the CoordinateMaps of the Wedge3Ds used in the Sphere, Shell, and binary compact object Dom...
Definition: DomainHelpers.cpp:499
constexpr size_t get_nth_bit(const size_t i, const size_t N)
Get the nth bit from the right, counting from zero.
Definition: ConstantExpressions.hpp:44
The type that options are passed around as. Contains YAML node data and an OptionContext.
Definition: Options.hpp:103
Use only the four equatorial wedges.
Each member in PairOfFaces holds the global corner ids of a block face. PairOfFaces is used in settin...
Definition: DomainHelpers.hpp:45
size_t collapsed_index(const Index< N > &index, const Index< N > &extents) noexcept
Get the collapsed index into a 1D array of the data corresponding to this Index. Note that the first ...
Use only the single wedge along -x.
Iterates over the corners of a VolumeDim-dimensional cube.
Definition: DomainHelpers.hpp:240
void set_identified_boundaries(const std::vector< PairOfFaces > &identifications, const std::vector< std::array< size_t, two_to_the(VolumeDim)>> &corners_of_all_blocks, gsl::not_null< std::vector< DirectionMap< VolumeDim, BlockNeighbor< VolumeDim >>> *> neighbors_of_all_blocks) noexcept
Sets up additional BlockNeighbors corresponding to any identifications of faces provided by the user...
Definition: DomainHelpers.cpp:400
Used by the parser to create an object. The default action is to parse options using T::options...
Definition: Options.hpp:143
std::vector< std::array< size_t, 8 > > corners_for_biradially_layered_domains(size_t number_of_radial_layers, size_t number_of_biradial_layers, bool include_central_block_lhs, bool include_central_block_rhs, const std::array< size_t, 8 > &central_block_corners_lhs={ {1, 2, 3, 4, 5, 6, 7, 8}}) noexcept
The corners for a domain with biradial layers.
Definition: DomainHelpers.cpp:719
Defines function make_array.
ShellWedges
The number of wedges to include in the Shell domain.
Definition: DomainHelpers.hpp:91
An optimized map with Direction keys.
Definition: DirectionMap.hpp:15
Side side() const noexcept
The side of the Direction.
Definition: Direction.hpp:46
Abstract base class for CoordinateMap.
Definition: CoordinateMap.hpp:49
auto operator*(const TensorExpression< T1, X, Symm1, IndexList1, Args1 > &t1, const TensorExpression< T2, X, Symm2, IndexList2, Args2 > &t2)
Definition: Product.hpp:89
constexpr auto create(Args &&... args)
Create a new DataBox.
Definition: DataBox.hpp:1259
Defines class template Index.
A particular Side along a particular coordinate Axis.
Definition: Direction.hpp:23
std::vector< std::array< size_t, 8 > > corners_for_radially_layered_domains(size_t number_of_layers, bool include_central_block, const std::array< size_t, 8 > &central_block_corners={{1, 2, 3, 4, 5, 6, 7, 8}}, ShellWedges which_wedges=ShellWedges::All) noexcept
The corners for a domain with radial layers.
Definition: DomainHelpers.cpp:686
std::array< size_t, two_to_the(VolumeDim)> discrete_rotation(const OrientationMap< VolumeDim > &orientation, const std::array< size_t, two_to_the(VolumeDim)> &corners_of_aligned) noexcept
Permutes the corner numbers of an n-cube.
Definition: DomainHelpers.cpp:935
T max(T... args)
constexpr T two_to_the(T n)
Compute 2 to the n for integral types.
Definition: ConstantExpressions.hpp:34
Define simple functions for constant expressions.
Domain< VolumeDim, TargetFrame > rectilinear_domain(const Index< VolumeDim > &domain_extents, const std::array< std::vector< double >, VolumeDim > &block_demarcations, const std::vector< Index< VolumeDim >> &block_indices_to_exclude={}, const std::vector< OrientationMap< VolumeDim >> &orientations_of_all_blocks={}, const std::array< bool, VolumeDim > &dimension_is_periodic=make_array< VolumeDim >(false), const std::vector< PairOfFaces > &identifications={}, bool use_equiangular_map=false) noexcept
Create a rectilinear Domain of multicubes.
Definition: DomainHelpers.cpp:975
Defines classes for Tensor.
Iterates over the 2^(VolumeDim-1) logical corners of the face of a VolumeDim-dimensional cube in the ...
Definition: DomainHelpers.hpp:325
auto frustum_coordinate_maps(double length_inner_cube, double length_outer_cube, bool use_equiangular_map) noexcept -> std::vector< std::unique_ptr< CoordinateMapBase< Frame::Logical, TargetFrame, 3 >>>
These are the ten Frustums used in the DomainCreators for binary compact objects. The cubes envelopin...
Definition: DomainHelpers.cpp:599
A wrapper around a vector of Blocks that represent the computational domain.
Definition: Domain.hpp:36
Information about the neighbor of a Block in a particular direction.
Definition: BlockNeighbor.hpp:25
Defines functions and classes from the GSL.
void set_internal_boundaries(const std::vector< std::array< size_t, two_to_the(VolumeDim)>> &corners_of_all_blocks, gsl::not_null< std::vector< DirectionMap< VolumeDim, BlockNeighbor< VolumeDim >>> *> neighbors_of_all_blocks) noexcept
Sets up the BlockNeighbors using the corner numbering scheme to deduce the correct neighbors and orie...
Definition: DomainHelpers.cpp:369
Require a pointer to not be a nullptr
Definition: ConservativeFromPrimitive.hpp:12
auto maps_for_rectilinear_domains(const Index< VolumeDim > &domain_extents, const std::array< std::vector< double >, VolumeDim > &block_demarcations, const std::vector< Index< VolumeDim >> &block_indices_to_exclude={}, const std::vector< OrientationMap< VolumeDim >> &orientations_of_all_blocks={}, bool use_equiangular_map=false) noexcept -> std::vector< std::unique_ptr< CoordinateMapBase< Frame::Logical, TargetFrame, VolumeDim >>>
The CoordinateMaps for a rectilinear domain of n-cubes.
Definition: DomainHelpers.cpp:864
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
Defines enum class Side.