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