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