Line data Source code
1 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"
17 : #include "DataStructures/Tensor/Tensor.hpp"
18 : #include "Domain/CoordinateMaps/Distribution.hpp"
19 : #include "Domain/Structure/Direction.hpp"
20 : #include "Domain/Structure/Side.hpp"
21 : #include "Utilities/ConstantExpressions.hpp"
22 : #include "Utilities/Gsl.hpp"
23 : #include "Utilities/MakeArray.hpp"
24 :
25 : /// \cond
26 : template <size_t VolumeDim>
27 : class BlockNeighbor;
28 : namespace domain {
29 : template <typename SourceFrame, typename TargetFrame, size_t Dim>
30 : class CoordinateMapBase;
31 : } // namespace domain
32 : template <size_t VolumeDim, typename T>
33 : class DirectionMap;
34 : template <size_t VolumeDim>
35 : class Domain;
36 : template <size_t VolumeDim>
37 : class OrientationMap;
38 : namespace Options {
39 : class Option;
40 : template <typename T>
41 : struct create_from_yaml;
42 : } // namespace Options
43 : namespace domain::CoordinateMaps {
44 : template <typename Map1, typename Map2>
45 : class ProductOf2Maps;
46 : template <typename Map1, typename Map2, typename Map3>
47 : class ProductOf3Maps;
48 : class Interval;
49 : template <size_t Dim>
50 : class Wedge;
51 : class Frustum;
52 : } // namespace domain::CoordinateMaps
53 : /// \endcond
54 :
55 : /// \ingroup ComputationalDomainGroup
56 : /// Each member in `PairOfFaces` holds the global corner ids of a block face.
57 : /// `PairOfFaces` is used in setting up periodic boundary conditions by
58 : /// identifying the two faces with each other.
59 : /// \requires The pair of faces must belong to a single block.
60 1 : struct PairOfFaces {
61 0 : std::vector<size_t> first;
62 0 : std::vector<size_t> second;
63 : };
64 :
65 : /// \ingroup ComputationalDomainGroup
66 : /// Sets up the BlockNeighbors using the corner numbering scheme
67 : /// provided by the user to deduce the correct neighbors and
68 : /// orientations. Does not set up periodic boundary conditions.
69 : template <size_t VolumeDim>
70 1 : void set_internal_boundaries(
71 : gsl::not_null<
72 : std::vector<DirectionMap<VolumeDim, BlockNeighbor<VolumeDim>>>*>
73 : neighbors_of_all_blocks,
74 : const std::vector<std::array<size_t, two_to_the(VolumeDim)>>&
75 : corners_of_all_blocks);
76 :
77 : /// \ingroup ComputationalDomainGroup
78 : /// Sets up the BlockNeighbors using the corner numbering scheme
79 : /// implied by the maps provided by the user to deduce the correct
80 : /// neighbors and orientations.
81 : /// \warning Does not set up periodic boundary conditions.
82 : template <size_t VolumeDim>
83 1 : void set_internal_boundaries(
84 : gsl::not_null<
85 : std::vector<DirectionMap<VolumeDim, BlockNeighbor<VolumeDim>>>*>
86 : neighbors_of_all_blocks,
87 : const std::vector<std::unique_ptr<domain::CoordinateMapBase<
88 : Frame::BlockLogical, Frame::Inertial, VolumeDim>>>& maps);
89 :
90 : /// \ingroup ComputationalDomainGroup
91 : /// Sets up additional BlockNeighbors corresponding to any
92 : /// identifications of faces provided by the user. Can be used
93 : /// for manually setting up periodic boundary conditions.
94 : template <size_t VolumeDim>
95 1 : void set_identified_boundaries(
96 : const std::vector<PairOfFaces>& identifications,
97 : const std::vector<std::array<size_t, two_to_the(VolumeDim)>>&
98 : corners_of_all_blocks,
99 : gsl::not_null<
100 : std::vector<DirectionMap<VolumeDim, BlockNeighbor<VolumeDim>>>*>
101 : neighbors_of_all_blocks);
102 :
103 : /// \ingroup ComputationalDomainGroup
104 : /// \brief The multi-indices that identify the individual Blocks in the lattice
105 : template <size_t VolumeDim>
106 1 : auto indices_for_rectilinear_domains(
107 : const Index<VolumeDim>& domain_extents,
108 : const std::vector<Index<VolumeDim>>& block_indices_to_exclude = {})
109 : -> std::vector<Index<VolumeDim>>;
110 :
111 : /// \ingroup ComputationalDomainGroup
112 : /// \brief The corners for a rectilinear domain made of n-cubes.
113 : ///
114 : /// The `domain_extents` argument holds the number of blocks to have
115 : /// in each dimension. The blocks all have aligned orientations by
116 : /// construction. The `block_indices_to_exclude` argument allows the user
117 : /// to selectively exclude blocks from the resulting domain. This allows
118 : /// for the creation of non-trivial shapes such as the net for a tesseract.
119 : template <size_t VolumeDim>
120 1 : auto corners_for_rectilinear_domains(
121 : const Index<VolumeDim>& domain_extents,
122 : const std::vector<Index<VolumeDim>>& block_indices_to_exclude = {})
123 : -> std::vector<std::array<size_t, two_to_the(VolumeDim)>>;
124 :
125 : /// \ingroup ComputationalDomainGroup
126 : /// The number of wedges to include in the Sphere domain.
127 1 : enum class ShellWedges {
128 : /// Use the entire shell
129 : All,
130 : /// Use only the four equatorial wedges
131 : FourOnEquator,
132 : /// Use only the single wedge along -x
133 : OneAlongMinusX
134 : };
135 :
136 : /// \ingroup ComputationalDomainGroup
137 : /// The first index in the list "UpperZ", "LowerZ", "UpperY", "LowerY", "UpperX"
138 : /// "LowerX" that is included in `which_wedges`. It is 0 for `ShellWedges::All`,
139 : /// 2 for `ShellWedges::FourOnEquator`, and 5 for `ShellWedges::OneAlongMinusX`.
140 1 : size_t which_wedge_index(const ShellWedges& which_wedges);
141 :
142 : /*!
143 : * \ingroup ComputationalDomainGroup
144 : * These are the CoordinateMaps of the Wedge<3>s used in the Sphere and
145 : * binary compact object DomainCreators. This function can also be used to
146 : * wrap the Sphere in a cube made of six Wedge<3>s.
147 : *
148 : * \param inner_radius Radius of the inner boundary of the shell, or the
149 : * radius circumscribing the inner cube of a sphere.
150 : * \param outer_radius Outer radius of the shell or sphere.
151 : * \param inner_sphericity Specifies if the wedges form a spherical inner
152 : * boundary (1.0) or a cubical inner boundary (0.0).
153 : * \param outer_sphericity Specifies if the wedges form a spherical outer
154 : * boundary (1.0) or a cubical outer boundary (0.0).
155 : * \param use_equiangular_map Toggles the equiangular map of the Wedge map.
156 : * \param use_half_wedges When `true`, the wedges in the +z,-z,+y,-y directions
157 : * are cut in half along their xi-axes. The resulting ten CoordinateMaps are
158 : * used for the outermost Blocks of the BBH Domain.
159 : * \param radial_partitioning Specifies the radial boundaries of sub-shells
160 : * between `inner_radius` and `outer_radius`. If the inner and outer
161 : * sphericities are different, the innermost shell does the transition.
162 : * \param radial_distribution Select the radial distribution of grid points in
163 : * the spherical shells.
164 : * \param which_wedges Select a subset of wedges.
165 : * \param opening_angle sets the combined opening angle of the two half wedges
166 : * that open up along the y-z plane. The endcap wedges are then given an angle
167 : * of pi minus this opening angle. This parameter only has an effect if
168 : * `use_half_wedges` is set to `true`.
169 : */
170 1 : std::vector<domain::CoordinateMaps::Wedge<3>> sph_wedge_coordinate_maps(
171 : double inner_radius, double outer_radius, double inner_sphericity,
172 : double outer_sphericity, bool use_equiangular_map,
173 : bool use_half_wedges = false,
174 : const std::vector<double>& radial_partitioning = {},
175 : const std::vector<domain::CoordinateMaps::Distribution>&
176 : radial_distribution = {domain::CoordinateMaps::Distribution::Linear},
177 : ShellWedges which_wedges = ShellWedges::All, double opening_angle = M_PI_2);
178 :
179 : /// \ingroup ComputationalDomainGroup
180 : /// These are the ten Frustums used in the DomainCreators for binary compact
181 : /// objects. The Frustums partition the volume defined by two bounding
182 : /// surfaces: The inner surface is the surface of the two joined inner cubes
183 : /// enveloping the two compact objects, while the outer is the surface of the
184 : /// outer cube. The cubes enveloping the two Shells each have a side length of
185 : /// `length_inner_cube`. The outer cube has a side length of
186 : /// `length_outer_cube`. `origin_preimage` is a parameter
187 : /// that moves the center of the two joined inner cubes away from the origin
188 : /// and to `-origin_preimage`. `projective_scale_factor` acts to change the
189 : /// gridpoint distribution in the radial direction. \see Frustum for details.
190 : /// The value for `sphericity` determines whether the outer surface is a cube
191 : /// (value of 0), a sphere (value of 1) or somewhere in between.
192 : /// The value for `opening_angle` determines the gridpoint distribution used
193 : /// in the Frustums such that they conform to the outer sphere of Wedges with
194 : /// the same value for `opening_angle`.
195 1 : std::vector<domain::CoordinateMaps::Frustum> frustum_coordinate_maps(
196 : double length_inner_cube, double length_outer_cube,
197 : bool use_equiangular_map,
198 : const std::array<double, 3>& origin_preimage = {{0.0, 0.0, 0.0}},
199 : double projective_scale_factor = 1.0, double sphericity = 0.0,
200 : double opening_angle = M_PI_2);
201 :
202 : /// \ingroup ComputationalDomainGroup
203 : /// \brief The corners for a domain with radial layers.
204 : ///
205 : /// Generates the corners for a Domain which is made of one or more layers
206 : /// of Blocks fully enveloping an interior volume, e.g. Sphere.
207 : ///
208 : /// \param number_of_layers specifies how many layers of Blocks to have
209 : /// in the final domain.
210 : /// \param include_central_block set to `true` where the interior
211 : /// volume is filled with a central Block, and `false` where the
212 : /// interior volume is left empty.
213 : /// \param central_block_corners are used as seed values to generate the corners
214 : /// for the surrounding Blocks.
215 : /// \param which_wedges can be used to exclude a subset of the wedges.
216 1 : std::vector<std::array<size_t, 8>> corners_for_radially_layered_domains(
217 : size_t number_of_layers, bool include_central_block,
218 : const std::array<size_t, 8>& central_block_corners = {{1, 2, 3, 4, 5, 6, 7,
219 : 8}},
220 : ShellWedges which_wedges = ShellWedges::All);
221 :
222 : /// \ingroup ComputationalDomainGroup
223 : /// \brief The corners for a domain with biradial layers.
224 : ///
225 : /// Generates the corners for a BBH-like Domain which is made of one or more
226 : /// layers of Blocks fully enveloping two interior volumes. The
227 : /// `number_of_radial_layers` gives the number of layers that fully envelop
228 : /// each interior volume with six Blocks each. The `number_of_biradial_layers`
229 : /// gives the number of layers that fully envelop both volumes at once, using
230 : /// ten Blocks per layer as opposed to six. The `central_block_corners_lhs`
231 : /// are used as seed values to generate the corners for the surrounding
232 : /// Blocks.
233 1 : std::vector<std::array<size_t, 8>> corners_for_biradially_layered_domains(
234 : size_t number_of_radial_layers, size_t number_of_biradial_layers,
235 : bool include_central_block_lhs, bool include_central_block_rhs,
236 : const std::array<size_t, 8>& central_block_corners_lhs = {
237 : {1, 2, 3, 4, 5, 6, 7, 8}});
238 :
239 : /// \ingroup ComputationalDomainGroup
240 : /// These are the CoordinateMaps used in the Cylinder DomainCreator.
241 : ///
242 : /// The `radial_partitioning` specifies the radial boundaries of sub-shells
243 : /// between `inner_radius` and `outer_radius`, while `partitioning_in_z`
244 : /// specifies the z-boundaries, splitting the cylinder into stacked
245 : /// 3-dimensional disks. The circularity of the shell wedges changes from 0 to 1
246 : /// within the innermost sub-shell.
247 : ///
248 : /// Set the `radial_distribution` to select the radial distribution of grid
249 : /// points in the cylindrical shells. The innermost shell must have
250 : /// `domain::CoordinateMaps::Distribution::Linear` because it changes the
251 : /// circularity. The distribution along the z-axis for each circular
252 : /// disc is specified through `distribution_in_z`.
253 : template <typename TargetFrame>
254 1 : auto cyl_wedge_coordinate_maps(
255 : double inner_radius, double outer_radius, double lower_z_bound,
256 : double upper_z_bound, bool use_equiangular_map,
257 : const std::vector<double>& radial_partitioning = {},
258 : const std::vector<double>& partitioning_in_z = {},
259 : const std::vector<domain::CoordinateMaps::Distribution>&
260 : radial_distribution = {domain::CoordinateMaps::Distribution::Linear},
261 : const std::vector<domain::CoordinateMaps::Distribution>& distribution_in_z =
262 : {domain::CoordinateMaps::Distribution::Linear})
263 : -> std::vector<std::unique_ptr<
264 : domain::CoordinateMapBase<Frame::BlockLogical, TargetFrame, 3>>>;
265 :
266 0 : enum class CylindricalDomainParityFlip { none, z_direction };
267 :
268 : /// \ingroup ComputationalDomainGroup
269 : /// Same as `cyl_wedge_coordinate_maps`, but only the center square blocks,
270 : ///
271 : /// If `CylindricalDomainParityFlip::z_direction` is specified, then
272 : /// the returned maps describe a cylinder with `lower_z_bound`
273 : /// corresponding to logical coordinate `upper_zeta` and `upper_z_bound`
274 : /// corresponding to logical coordinate `lower_zeta`, and thus the
275 : /// resulting maps are left-handed.
276 : /// `CylindricalDomainParityFlip::z_direction` is therefore useful
277 : /// only when composing with another map that is also left-handed, so
278 : /// that the composed coordinate system is right-handed.
279 : ///
280 : /// Returned as a vector of the coordinate maps so that they can
281 : /// be composed with other maps later.
282 1 : auto cyl_wedge_coord_map_center_blocks(
283 : double inner_radius, double lower_z_bound, double upper_z_bound,
284 : bool use_equiangular_map, const std::vector<double>& partitioning_in_z = {},
285 : const std::vector<domain::CoordinateMaps::Distribution>& distribution_in_z =
286 : {domain::CoordinateMaps::Distribution::Linear},
287 : CylindricalDomainParityFlip parity_flip = CylindricalDomainParityFlip::none)
288 : -> std::vector<domain::CoordinateMaps::ProductOf3Maps<
289 : domain::CoordinateMaps::Interval, domain::CoordinateMaps::Interval,
290 : domain::CoordinateMaps::Interval>>;
291 :
292 : /// \ingroup ComputationalDomainGroup
293 : /// Same as cyl_wedge_coordinate_maps, but only the surrounding wedge blocks.
294 : ///
295 : /// If `CylindricalDomainParityFlip::z_direction` is specified, then
296 : /// the returned maps describe a cylinder with `lower_z_bound`
297 : /// corresponding to logical coordinate `upper_zeta` and `upper_z_bound`
298 : /// corresponding to logical coordinate `lower_zeta`, and thus the
299 : /// resulting maps are left-handed.
300 : /// `CylindricalDomainParityFlip::z_direction` is therefore useful
301 : /// only when composing with another map that is also left-handed, so
302 : /// that the composed coordinate system is right-handed.
303 : ///
304 : /// Returned as a vector of the coordinate maps so that they can
305 : /// be composed with other maps later.
306 1 : auto cyl_wedge_coord_map_surrounding_blocks(
307 : double inner_radius, double outer_radius, double lower_z_bound,
308 : double upper_z_bound, bool use_equiangular_map, double inner_circularity,
309 : const std::vector<double>& radial_partitioning = {},
310 : const std::vector<double>& partitioning_in_z = {},
311 : const std::vector<domain::CoordinateMaps::Distribution>&
312 : radial_distribution = {domain::CoordinateMaps::Distribution::Linear},
313 : const std::vector<domain::CoordinateMaps::Distribution>& distribution_in_z =
314 : {domain::CoordinateMaps::Distribution::Linear},
315 : CylindricalDomainParityFlip parity_flip = CylindricalDomainParityFlip::none)
316 : -> std::vector<domain::CoordinateMaps::ProductOf2Maps<
317 : domain::CoordinateMaps::Wedge<2>, domain::CoordinateMaps::Interval>>;
318 :
319 : /// \ingroup ComputationalDomainGroup
320 : /// \brief The corners for a cylindrical domain split into discs with radial
321 : /// shells.
322 : ///
323 : /// Generates the corners for a Domain which is made of one or more stacked
324 : /// discs consisting of layers of Blocks enveloping an interior square prism.
325 : /// The `number_of_shells` specifies how many of these layers of Blocks to have
326 : /// in each disc.
327 : ///
328 : /// The `number_of_discs` specifies how many discs make up the domain.
329 : /// The very basic cylinder with one shell and one layer serves as a base
330 : /// to generate the corners for subsequent shells first and discs second.
331 1 : std::vector<std::array<size_t, 8>> corners_for_cylindrical_layered_domains(
332 : size_t number_of_shells, size_t number_of_discs);
333 :
334 : /// \ingroup ComputationalDomainGroup
335 : /// \brief Permutes the corner numbers of an n-cube.
336 : ///
337 : /// Returns the correct ordering of global corner numbers for a rotated block
338 : /// in an otherwise aligned edifice of blocks, given the OrientationMap a
339 : /// block aligned with the edifice has relative to this one, and given the
340 : /// corner numbering the rotated block would have if it were aligned.
341 : /// This is useful in creating domains for testing purposes, e.g.
342 : /// RotatedIntervals, RotatedRectangles, and RotatedBricks.
343 : template <size_t VolumeDim>
344 1 : std::array<size_t, two_to_the(VolumeDim)> discrete_rotation(
345 : const OrientationMap<VolumeDim>& orientation,
346 : const std::array<size_t, two_to_the(VolumeDim)>& corners_of_aligned);
347 :
348 : /// \ingroup ComputationalDomainGroup
349 : /// \brief The CoordinateMaps for a rectilinear domain of n-cubes.
350 : ///
351 : /// Allows for both Affine and Equiangular maps.
352 : template <typename TargetFrame, size_t VolumeDim>
353 1 : auto maps_for_rectilinear_domains(
354 : const Index<VolumeDim>& domain_extents,
355 : const std::array<std::vector<double>, VolumeDim>& block_demarcations,
356 : const std::vector<Index<VolumeDim>>& block_indices_to_exclude = {},
357 : const std::vector<OrientationMap<VolumeDim>>& orientations_of_all_blocks =
358 : {},
359 : bool use_equiangular_map = false)
360 : -> std::vector<std::unique_ptr<domain::CoordinateMapBase<
361 : Frame::BlockLogical, TargetFrame, VolumeDim>>>;
362 :
363 : /// \ingroup ComputationalDomainGroup
364 : /// \brief Create a rectilinear Domain of multicubes.
365 : ///
366 : /// \details Useful for constructing domains for testing non-trivially
367 : /// connected rectilinear domains made up of cubes. We refer to a domain of
368 : /// this type as an edifice. The `domain_extents` provides the size (in the
369 : /// number of blocks) of the initial aligned edifice to construct. The
370 : /// `block_indices_to_exclude` parameter is used in refining the shape of
371 : /// the edifice from a cube to sometime more non-trivial, such as an L-shape
372 : /// or the net of a tesseract. The `block_demarcations` and
373 : /// `use_equiangular_map` parameters determine the CoordinateMaps to be used.
374 : /// `orientations_of_all_blocks` contains the OrientationMap of the edifice
375 : /// relative to each block.
376 : ///
377 : /// The `identifications` parameter is used when identifying the faces of
378 : /// blocks in an edifice. This is used to identify the 1D boundaries in the 2D
379 : /// net for a 3D cube to construct a domain with topology S2. Note: If the user
380 : /// wishes to rotate the blocks as well as manually identify their faces, the
381 : /// user must provide the PairOfFaces corresponding to the rotated corners.
382 : template <size_t VolumeDim>
383 1 : Domain<VolumeDim> rectilinear_domain(
384 : const Index<VolumeDim>& domain_extents,
385 : const std::array<std::vector<double>, VolumeDim>& block_demarcations,
386 : const std::vector<Index<VolumeDim>>& block_indices_to_exclude = {},
387 : const std::vector<OrientationMap<VolumeDim>>& orientations_of_all_blocks =
388 : {},
389 : const std::array<bool, VolumeDim>& dimension_is_periodic =
390 : make_array<VolumeDim>(false),
391 : const std::vector<PairOfFaces>& identifications = {},
392 : bool use_equiangular_map = false);
393 :
394 : /// \ingroup ComputationalDomainGroup
395 : /// Iterates over the corners of a VolumeDim-dimensional cube.
396 : template <size_t VolumeDim>
397 1 : class VolumeCornerIterator {
398 : public:
399 0 : VolumeCornerIterator() { setup_from_local_corner_number(); }
400 :
401 0 : explicit VolumeCornerIterator(size_t initial_local_corner_number)
402 : : local_corner_number_(initial_local_corner_number) {
403 : setup_from_local_corner_number();
404 : }
405 0 : VolumeCornerIterator(
406 : // The block index is also global corner
407 : // index of the lowest corner of the block.
408 : Index<VolumeDim> block_index, Index<VolumeDim> global_corner_extents)
409 : : global_corner_number_(
410 : collapsed_index(block_index, global_corner_extents)),
411 : global_corner_index_(block_index),
412 : global_corner_extents_(global_corner_extents) {}
413 :
414 0 : void operator++() {
415 : ++local_corner_number_;
416 : setup_from_local_corner_number();
417 : }
418 :
419 0 : explicit operator bool() const {
420 : return local_corner_number_ < two_to_the(VolumeDim);
421 : }
422 :
423 0 : size_t local_corner_number() const { return local_corner_number_; }
424 :
425 0 : size_t global_corner_number() const {
426 : std::array<size_t, VolumeDim> new_indices{};
427 : for (size_t i = 0; i < VolumeDim; i++) {
428 : gsl::at(new_indices, i) =
429 : global_corner_index_[i] +
430 : (gsl::at(array_sides_, i) == Side::Upper ? 1 : 0);
431 : }
432 : const Index<VolumeDim> interior_multi_index(new_indices);
433 : return collapsed_index(interior_multi_index, global_corner_extents_);
434 : }
435 :
436 0 : const std::array<Side, VolumeDim>& operator()() const { return array_sides_; }
437 :
438 0 : const std::array<Side, VolumeDim>& operator*() const { return array_sides_; }
439 :
440 0 : const std::array<double, VolumeDim>& coords_of_corner() const {
441 : return coords_of_corner_;
442 : }
443 :
444 0 : const std::array<Direction<VolumeDim>, VolumeDim>& directions_of_corner()
445 : const {
446 : return array_directions_;
447 : }
448 :
449 0 : void setup_from_local_corner_number() {
450 : for (size_t i = 0; i < VolumeDim; i++) {
451 : gsl::at(coords_of_corner_, i) =
452 : 2.0 * get_nth_bit(local_corner_number_, i) - 1.0;
453 : gsl::at(array_sides_, i) =
454 : 2 * get_nth_bit(local_corner_number_, i) - 1 == 1 ? Side::Upper
455 : : Side::Lower;
456 : gsl::at(array_directions_, i) =
457 : Direction<VolumeDim>(i, gsl::at(array_sides_, i));
458 : }
459 : }
460 :
461 : private:
462 0 : size_t local_corner_number_ = 0;
463 0 : size_t global_corner_number_{std::numeric_limits<size_t>::max()};
464 0 : Index<VolumeDim> global_corner_index_{};
465 0 : Index<VolumeDim> global_corner_extents_{};
466 0 : std::array<Side, VolumeDim> array_sides_ = make_array<VolumeDim>(Side::Lower);
467 0 : std::array<Direction<VolumeDim>, VolumeDim> array_directions_{};
468 0 : std::array<double, VolumeDim> coords_of_corner_ = make_array<VolumeDim>(-1.0);
469 : };
470 :
471 : /// \ingroup ComputationalDomainGroup
472 : /// Iterates over the 2^(VolumeDim-1) logical corners of the face of a
473 : /// VolumeDim-dimensional cube in the given direction.
474 : template <size_t VolumeDim>
475 1 : class FaceCornerIterator {
476 : public:
477 0 : explicit FaceCornerIterator(Direction<VolumeDim> direction);
478 :
479 0 : void operator++() {
480 : face_index_++;
481 : do {
482 : index_++;
483 : } while (get_nth_bit(index_, direction_.dimension()) ==
484 : (direction_.side() == Side::Upper ? 0 : 1));
485 : for (size_t i = 0; i < VolumeDim; ++i) {
486 : corner_[i] = 2 * static_cast<int>(get_nth_bit(index_, i)) - 1;
487 : }
488 : }
489 :
490 0 : explicit operator bool() const {
491 : return face_index_ < two_to_the(VolumeDim - 1);
492 : }
493 :
494 0 : tnsr::I<double, VolumeDim, Frame::BlockLogical> operator()() const {
495 : return corner_;
496 : }
497 :
498 0 : tnsr::I<double, VolumeDim, Frame::BlockLogical> operator*() const {
499 : return corner_;
500 : }
501 :
502 : // Returns the value used to construct the logical corner.
503 0 : size_t volume_index() const { return index_; }
504 :
505 : // Returns the number of times operator++ has been called.
506 0 : size_t face_index() const { return face_index_; }
507 :
508 : private:
509 0 : const Direction<VolumeDim> direction_;
510 0 : size_t index_;
511 0 : size_t face_index_ = 0;
512 0 : tnsr::I<double, VolumeDim, Frame::BlockLogical> corner_;
513 : };
514 :
515 : template <size_t VolumeDim>
516 : FaceCornerIterator<VolumeDim>::FaceCornerIterator(
517 : Direction<VolumeDim> direction)
518 : : direction_(std::move(direction)),
519 : index_(direction_.side() == Side::Upper
520 : ? two_to_the(direction_.dimension())
521 : : 0) {
522 : for (size_t i = 0; i < VolumeDim; ++i) {
523 : corner_[i] = 2 * static_cast<int>(get_nth_bit(index_, i)) - 1;
524 : }
525 : }
526 :
527 0 : std::ostream& operator<<(std::ostream& os, const ShellWedges& which_wedges);
528 :
529 : template <>
530 0 : struct Options::create_from_yaml<ShellWedges> {
531 : template <typename Metavariables>
532 0 : static ShellWedges create(const Options::Option& options) {
533 : return create<void>(options);
534 : }
535 : };
536 : template <>
537 0 : ShellWedges Options::create_from_yaml<ShellWedges>::create<void>(
538 : const Options::Option& options);
|
