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