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