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 : /*!
207 : * \ingroup ComputationalDomainGroup
208 : * \brief Redistribute points in the radial ($\zeta$) direction of a domain
209 : * layer (e.g. a shell or cube around an object) between the radii
210 : * \f$r_\mathrm{in}\f$ and \f$r_\mathrm{out}\f$.
211 : *
212 : * The map takes a logical coordinate in \f$[-1, 1]\f$ to a
213 : * logarithmically-redistributed logical coordinate in \f$[-1, 1]\f$. It is
214 : * meant to be applied before a map that is affine along each radial ray, such
215 : * as a `domain::CoordinateMaps::Wedge` map. A `log_map_strength` of 1.0 places
216 : * the log-map physical origin ("singularity position") at radius 0; larger
217 : * values move it toward the inner radius, concentrating grid points near the
218 : * object.
219 : *
220 : * We define three radial coordinates:
221 : * - \f$\zeta\f$: input logical coordinate in \f$[-1, 1]\f$;
222 : * - \f$\zeta_{\log}\f$: logarithmically-redistributed logical coordinate in
223 : * \f$[-1, 1]\f$;
224 : * - \f$r\f$: physical radius in \f$[r_\mathrm{in}, r_\mathrm{out}]\f$.
225 : *
226 : * Hence, we define two maps:
227 : * - Log map \f$\zeta_{\log}(\zeta)\f$:
228 : * \f{align}
229 : * \ln(\zeta_{\log} - \zeta_0)
230 : * &= \frac{1 + \zeta}{2} \ln(1 - \zeta_0) \nonumber \\
231 : * &+ \frac{1 - \zeta}{2} \ln(-1 - \zeta_0),
232 : * \f}
233 : * where \f$\zeta_0 < -1\f$ is the "logical position" of the log map's
234 : * origin.
235 : * - Radius map \f$r(\zeta_{\log})\f$, affine along each radial ray:
236 : * \f{align}
237 : * r(\zeta_{\log}) = \frac{1 - \zeta_{\log}}{2} r_\mathrm{in}
238 : * + \frac{1 + \zeta_{\log}}{2} r_\mathrm{out}.
239 : * \f}
240 : * For example, for a `domain::CoordinateMaps::Wedge` map between two
241 : * spherical surfaces, the redistribution is isotropic. With a flat outer
242 : * surface, \f$r_\mathrm{out}\f$ is the distance to the corners, so the
243 : * redistribution varies with direction.
244 : *
245 : * Composing these two maps, we get \f$r(\zeta) = r(\zeta_{\log}(\zeta))\f$:
246 : *
247 : * \f{align}
248 : * \ln(r - r_0) &= \frac{1 + \zeta}{2} \ln(r_\mathrm{out} - r_0) \nonumber \\
249 : * &+ \frac{1 - \zeta}{2} \ln(r_\mathrm{in} - r_0),
250 : * \f}
251 : *
252 : * where \f$r_0 < r_\mathrm{in}\f$ is the "physical position" of the log map's
253 : * origin. For the log map \f$\zeta_{\log}(\zeta)\f$, we need to specify
254 : * \f$\zeta_0\f$, but it is easier to specify \f$r_0\f$ because the
255 : * ratio of the radial Jacobian (\f$J := dr/d\zeta\f$) at the outer
256 : * (\f$\zeta=+1\f$) vs inner (\f$\zeta=-1\f$) boundaries is
257 : *
258 : * \f{align}
259 : * \frac{J_\mathrm{out}}{J_\mathrm{in}}
260 : * = \frac{r_\mathrm{out} - r_0}{r_\mathrm{in} - r_0}
261 : * := \alpha \frac{r_\mathrm{out}}{r_\mathrm{in}},
262 : * \f}
263 : *
264 : * where \f$\alpha \geq 1\f$ is `log_map_strength`, which we use to
265 : * parametrize \f$r_0\f$. Once we specify \f$\alpha\f$, we can invert the
266 : * Jacobian ratio above to get
267 : *
268 : * \f{align}
269 : * r_0 = \frac{r_\mathrm{in} r_\mathrm{out} (1 - \alpha)}
270 : * {r_\mathrm{in} - \alpha r_\mathrm{out}}.
271 : * \f}
272 : *
273 : * Then, we can invert the radius map \f$r(\zeta_{\log})\f$ to get
274 : *
275 : * \f{align}
276 : * \zeta_0 = \frac{2 r_0 - (r_\mathrm{out} + r_\mathrm{in})}
277 : * {r_\mathrm{out} - r_\mathrm{in}}.
278 : * \f}
279 : *
280 : * \param inner_radius The inner radius \f$r_\mathrm{in}\f$ of the layer.
281 : * \param outer_radius The outer radius \f$r_\mathrm{out}\f$ of the layer.
282 : * \param log_map_strength The log-map strength \f$\alpha \geq 1\f$.
283 : */
284 1 : domain::CoordinateMaps::Interval radial_distribution_log_map(
285 : double inner_radius, double outer_radius, double log_map_strength);
286 :
287 : /// \ingroup ComputationalDomainGroup
288 : /// These are the ten Frustums used in the DomainCreators for binary compact
289 : /// objects. The Frustums partition the volume defined by two bounding
290 : /// surfaces: The inner surface is the surface of the two joined inner cubes
291 : /// enveloping the two compact objects, while the outer is the surface of the
292 : /// outer cube.
293 : ///
294 : /// When the sphericity is 0, the \p length_inner_cube must be less than $1/2$
295 : /// \p length_outer_cube while when the sphericity is 1 it must be less than
296 : /// $\sqrt{3}/2$ \p length_outer_cube.
297 : ///
298 : /// \param length_inner_cube The side length of the cubes enveloping the two
299 : /// shells.
300 : /// \param length_outer_cube The side length of the outer cube.
301 : /// \param equiangular_map_at_outer Whether to apply a tangent map in the
302 : /// angular directions at the outer boundary.
303 : /// \param equiangular_map_at_inner Whether to apply a tangent map in the
304 : /// angular directions at the inner boundary.
305 : /// \param origin_preimage The center of the two joined inner cubes is moved
306 : /// away from the origin and to this point, origin_preimage.
307 : /// \param radial_distribution The gridpoint distribution in the radial
308 : /// direction, possibly dependent on the value passed to `distribution_value`.
309 : /// \param distribution_value Used by `radial_distribution`. \see Frustum for
310 : /// details.
311 : /// \param sphericity Determines whether the outer surface is a cube
312 : /// (value of 0), a sphere (value of 1) or somewhere in between.
313 : /// \param opening_angle determines the gridpoint distribution used
314 : /// in the Frustums such that they conform to the outer sphere of Wedges with
315 : /// the same value for `opening_angle`.
316 1 : std::vector<domain::CoordinateMaps::Frustum> frustum_coordinate_maps(
317 : double length_inner_cube, double length_outer_cube,
318 : bool equiangular_map_at_outer, bool equiangular_map_at_inner,
319 : const std::array<double, 3>& origin_preimage = {{0.0, 0.0, 0.0}},
320 : domain::CoordinateMaps::Distribution radial_distribution =
321 : domain::CoordinateMaps::Distribution::Linear,
322 : std::optional<double> distribution_value = std::nullopt,
323 : double sphericity = 0.0, double opening_angle = M_PI_2);
324 :
325 : /// \ingroup ComputationalDomainGroup
326 : /// \brief The corners for a domain with radial layers.
327 : ///
328 : /// Generates the corners for a Domain which is made of one or more layers
329 : /// of Blocks fully enveloping an interior volume, e.g. Sphere.
330 : ///
331 : /// \param number_of_layers specifies how many layers of Blocks to have
332 : /// in the final domain.
333 : /// \param include_central_block set to `true` where the interior
334 : /// volume is filled with a central Block, and `false` where the
335 : /// interior volume is left empty.
336 : /// \param central_block_corners are used as seed values to generate the corners
337 : /// for the surrounding Blocks.
338 : /// \param which_wedges can be used to exclude a subset of the wedges.
339 1 : std::vector<std::array<size_t, 8>> corners_for_radially_layered_domains(
340 : size_t number_of_layers, bool include_central_block,
341 : const std::array<size_t, 8>& central_block_corners = {{1, 2, 3, 4, 5, 6, 7,
342 : 8}},
343 : ShellWedges which_wedges = ShellWedges::All);
344 :
345 : /// \ingroup ComputationalDomainGroup
346 : /// \brief The corners for a domain with biradial layers.
347 : ///
348 : /// Generates the corners for a BBH-like Domain which is made of one or more
349 : /// layers of Blocks fully enveloping two interior volumes. The
350 : /// `number_of_radial_layers` gives the number of layers that fully envelop
351 : /// each interior volume with six Blocks each. The `number_of_biradial_layers`
352 : /// gives the number of layers that fully envelop both volumes at once, using
353 : /// ten Blocks per layer as opposed to six. The `central_block_corners_lhs`
354 : /// are used as seed values to generate the corners for the surrounding
355 : /// Blocks.
356 1 : std::vector<std::array<size_t, 8>> corners_for_biradially_layered_domains(
357 : size_t number_of_radial_layers, size_t number_of_biradial_layers,
358 : bool include_central_block_lhs, bool include_central_block_rhs,
359 : const std::array<size_t, 8>& central_block_corners_lhs = {
360 : {1, 2, 3, 4, 5, 6, 7, 8}});
361 :
362 : /// \ingroup ComputationalDomainGroup
363 : /// These are the CoordinateMaps used in the Cylinder DomainCreator for when
364 : /// cylinders are built using cubes instead of a single cylindrical block.
365 : ///
366 : /// The `radial_partitioning` specifies the radial boundaries of sub-shells
367 : /// between `inner_radius` and `outer_radius`, while `partitioning_in_z`
368 : /// specifies the z-boundaries, splitting the cylinder into stacked
369 : /// 3-dimensional disks. The circularity of the shell wedges changes from 0 to 1
370 : /// within the innermost sub-shell.
371 : ///
372 : /// Set the `radial_distribution` to select the radial distribution of grid
373 : /// points in the cylindrical shells. The innermost shell must have
374 : /// `domain::CoordinateMaps::Distribution::Linear` because it changes the
375 : /// circularity. The distribution along the z-axis for each circular
376 : /// disc is specified through `distribution_in_z`.
377 : template <typename TargetFrame>
378 1 : auto cyl_wedge_coordinate_maps(
379 : double inner_radius, double outer_radius, double lower_z_bound,
380 : double upper_z_bound, bool use_equiangular_map,
381 : const std::vector<double>& radial_partitioning = {},
382 : const std::vector<double>& partitioning_in_z = {},
383 : const std::vector<domain::CoordinateMaps::Distribution>&
384 : radial_distribution = {domain::CoordinateMaps::Distribution::Linear},
385 : const std::vector<domain::CoordinateMaps::Distribution>& distribution_in_z =
386 : {domain::CoordinateMaps::Distribution::Linear})
387 : -> std::vector<std::unique_ptr<
388 : domain::CoordinateMapBase<Frame::BlockLogical, TargetFrame, 3>>>;
389 :
390 0 : enum class CylindricalDomainParityFlip { none, z_direction };
391 :
392 : /// \ingroup ComputationalDomainGroup
393 : /// Same as `cyl_wedge_coordinate_maps`, but only the center square blocks,
394 : ///
395 : /// If `CylindricalDomainParityFlip::z_direction` is specified, then
396 : /// the returned maps describe a cylinder with `lower_z_bound`
397 : /// corresponding to logical coordinate `upper_zeta` and `upper_z_bound`
398 : /// corresponding to logical coordinate `lower_zeta`, and thus the
399 : /// resulting maps are left-handed.
400 : /// `CylindricalDomainParityFlip::z_direction` is therefore useful
401 : /// only when composing with another map that is also left-handed, so
402 : /// that the composed coordinate system is right-handed.
403 : ///
404 : /// Returned as a vector of the coordinate maps so that they can
405 : /// be composed with other maps later.
406 1 : auto cyl_wedge_coord_map_center_blocks(
407 : double inner_radius, double lower_z_bound, double upper_z_bound,
408 : bool use_equiangular_map, const std::vector<double>& partitioning_in_z = {},
409 : const std::vector<domain::CoordinateMaps::Distribution>& distribution_in_z =
410 : {domain::CoordinateMaps::Distribution::Linear},
411 : CylindricalDomainParityFlip parity_flip = CylindricalDomainParityFlip::none)
412 : -> std::vector<domain::CoordinateMaps::ProductOf3Maps<
413 : domain::CoordinateMaps::Interval, domain::CoordinateMaps::Interval,
414 : domain::CoordinateMaps::Interval>>;
415 :
416 : /// \ingroup ComputationalDomainGroup
417 : /// Same as cyl_wedge_coordinate_maps, but only the surrounding wedge blocks.
418 : ///
419 : /// If `CylindricalDomainParityFlip::z_direction` is specified, then
420 : /// the returned maps describe a cylinder with `lower_z_bound`
421 : /// corresponding to logical coordinate `upper_zeta` and `upper_z_bound`
422 : /// corresponding to logical coordinate `lower_zeta`, and thus the
423 : /// resulting maps are left-handed.
424 : /// `CylindricalDomainParityFlip::z_direction` is therefore useful
425 : /// only when composing with another map that is also left-handed, so
426 : /// that the composed coordinate system is right-handed.
427 : ///
428 : /// Returned as a vector of the coordinate maps so that they can
429 : /// be composed with other maps later.
430 1 : auto cyl_wedge_coord_map_surrounding_blocks(
431 : double inner_radius, double outer_radius, double lower_z_bound,
432 : double upper_z_bound, bool use_equiangular_map, double inner_circularity,
433 : const std::vector<double>& radial_partitioning = {},
434 : const std::vector<double>& partitioning_in_z = {},
435 : const std::vector<domain::CoordinateMaps::Distribution>&
436 : radial_distribution = {domain::CoordinateMaps::Distribution::Linear},
437 : const std::vector<domain::CoordinateMaps::Distribution>& distribution_in_z =
438 : {domain::CoordinateMaps::Distribution::Linear},
439 : CylindricalDomainParityFlip parity_flip = CylindricalDomainParityFlip::none)
440 : -> std::vector<domain::CoordinateMaps::ProductOf2Maps<
441 : domain::CoordinateMaps::Wedge<2>, domain::CoordinateMaps::Interval>>;
442 :
443 : /// \ingroup ComputationalDomainGroup
444 : /// This is the CoordinateMap used in the Cylinder DomainCreator for when
445 : /// cylinders are built using a single cylindrical block instead of cubes.
446 : ///
447 : /// Returns a unit cylinder with the given `inner_radius`, `outer_radius`,
448 : /// `lower_z_bound`, and `upper_z_bound`. The returned unit cylinder's
449 : /// intended use is to compose it with `UniformCylindricalEndCap`,
450 : /// `UniformCylindricalFlatEndCap`, or `UniformCylindricalSide` to create the
451 : /// different cylinders needed by `domain::CylindricalBinaryCompactObject`.
452 : ::domain::CoordinateMap<
453 : Frame::BlockLogical, Frame::Inertial,
454 : ::domain::CoordinateMaps::ProductOf3Maps<
455 : ::domain::CoordinateMaps::Affine, ::domain::CoordinateMaps::Identity<1>,
456 : ::domain::CoordinateMaps::Interval>,
457 : ::domain::CoordinateMaps::ProductOf2Maps<
458 : ::domain::CoordinateMaps::PolarToCartesian,
459 : ::domain::CoordinateMaps::Identity<1>>>
460 1 : cyl_coordinate_map(double inner_radius, double outer_radius,
461 : double lower_z_bound, double upper_z_bound);
462 :
463 : /// \ingroup ComputationalDomainGroup
464 : /// \brief The corners for a cylindrical domain split into discs with radial
465 : /// shells.
466 : ///
467 : /// Generates the corners for a Domain which is made of one or more stacked
468 : /// discs consisting of layers of Blocks enveloping an interior square prism.
469 : /// The `number_of_shells` specifies how many of these layers of Blocks to have
470 : /// in each disc.
471 : ///
472 : /// The `number_of_discs` specifies how many discs make up the domain.
473 : /// The very basic cylinder with one shell and one layer serves as a base
474 : /// to generate the corners for subsequent shells first and discs second.
475 1 : std::vector<std::array<size_t, 8>> corners_for_cylindrical_layered_domains(
476 : size_t number_of_shells, size_t number_of_discs);
477 :
478 : /// \ingroup ComputationalDomainGroup
479 : /// \brief Permutes the corner numbers of an n-cube.
480 : ///
481 : /// Returns the correct ordering of global corner numbers for a rotated block
482 : /// in an otherwise aligned edifice of blocks, given the OrientationMap a
483 : /// block aligned with the edifice has relative to this one, and given the
484 : /// corner numbering the rotated block would have if it were aligned.
485 : /// This is useful in creating domains for testing purposes, e.g.
486 : /// RotatedIntervals, RotatedRectangles, and RotatedBricks.
487 : template <size_t VolumeDim>
488 1 : std::array<size_t, two_to_the(VolumeDim)> discrete_rotation(
489 : const OrientationMap<VolumeDim>& orientation,
490 : const std::array<size_t, two_to_the(VolumeDim)>& corners_of_aligned);
491 :
492 : /// \ingroup ComputationalDomainGroup
493 : /// \brief The CoordinateMaps for a rectilinear domain of n-cubes.
494 : ///
495 : /// Allows for both Affine and Equiangular maps.
496 : template <typename TargetFrame, size_t VolumeDim>
497 1 : auto maps_for_rectilinear_domains(
498 : const Index<VolumeDim>& domain_extents,
499 : const std::array<std::vector<double>, VolumeDim>& block_demarcations,
500 : const std::vector<Index<VolumeDim>>& block_indices_to_exclude = {},
501 : const std::vector<OrientationMap<VolumeDim>>& orientations_of_all_blocks =
502 : {},
503 : const std::array<std::vector<domain::CoordinateMaps::Distribution>,
504 : VolumeDim>& distributions = {},
505 : const std::array<std::vector<std::optional<double>>, VolumeDim>&
506 : singularity_positions = {},
507 : bool use_equiangular_map = false)
508 : -> std::vector<std::unique_ptr<domain::CoordinateMapBase<
509 : Frame::BlockLogical, TargetFrame, VolumeDim>>>;
510 :
511 : /// \ingroup ComputationalDomainGroup
512 : /// \brief Create a rectilinear Domain of multicubes.
513 : ///
514 : /// \details Useful for constructing domains for testing non-trivially
515 : /// connected rectilinear domains made up of cubes. We refer to a domain of
516 : /// this type as an edifice. The `domain_extents` provides the size (in the
517 : /// number of blocks) of the initial aligned edifice to construct. The
518 : /// `block_indices_to_exclude` parameter is used in refining the shape of
519 : /// the edifice from a cube to sometime more non-trivial, such as an L-shape
520 : /// or the net of a tesseract. The `block_demarcations` and
521 : /// `use_equiangular_map` parameters determine the CoordinateMaps to be used.
522 : /// `orientations_of_all_blocks` contains the OrientationMap of the edifice
523 : /// relative to each block.
524 : ///
525 : /// The `identifications` parameter is used when identifying the faces of
526 : /// blocks in an edifice. This is used to identify the 1D boundaries in the 2D
527 : /// net for a 3D cube to construct a domain with topology S2. Note: If the user
528 : /// wishes to rotate the blocks as well as manually identify their faces, the
529 : /// user must provide the PairOfFaces corresponding to the rotated corners.
530 : template <size_t VolumeDim>
531 1 : Domain<VolumeDim> rectilinear_domain(
532 : const Index<VolumeDim>& domain_extents,
533 : const std::array<std::vector<double>, VolumeDim>& block_demarcations,
534 : const std::vector<Index<VolumeDim>>& block_indices_to_exclude = {},
535 : const std::vector<OrientationMap<VolumeDim>>& orientations_of_all_blocks =
536 : {},
537 : const std::array<bool, VolumeDim>& dimension_is_periodic =
538 : make_array<VolumeDim>(false),
539 : const std::vector<PairOfFaces>& identifications = {},
540 : const std::array<std::vector<domain::CoordinateMaps::Distribution>,
541 : VolumeDim>& distributions = {},
542 : const std::array<std::vector<std::optional<double>>, VolumeDim>&
543 : singularity_positions = {},
544 : bool use_equiangular_map = false);
545 :
546 : /// \ingroup ComputationalDomainGroup
547 : /// Iterates over the corners of a VolumeDim-dimensional cube.
548 : template <size_t VolumeDim>
549 1 : class VolumeCornerIterator {
550 : public:
551 0 : VolumeCornerIterator() { setup_from_local_corner_number(); }
552 :
553 0 : explicit VolumeCornerIterator(size_t initial_local_corner_number)
554 : : local_corner_number_(initial_local_corner_number) {
555 : setup_from_local_corner_number();
556 : }
557 0 : VolumeCornerIterator(
558 : // The block index is also global corner
559 : // index of the lowest corner of the block.
560 : Index<VolumeDim> block_index, Index<VolumeDim> global_corner_extents)
561 : : global_corner_number_(
562 : collapsed_index(block_index, global_corner_extents)),
563 : global_corner_index_(block_index),
564 : global_corner_extents_(global_corner_extents) {}
565 :
566 0 : void operator++() {
567 : ++local_corner_number_;
568 : setup_from_local_corner_number();
569 : }
570 :
571 0 : explicit operator bool() const {
572 : return local_corner_number_ < two_to_the(VolumeDim);
573 : }
574 :
575 0 : size_t local_corner_number() const { return local_corner_number_; }
576 :
577 0 : size_t global_corner_number() const {
578 : std::array<size_t, VolumeDim> new_indices{};
579 : for (size_t i = 0; i < VolumeDim; i++) {
580 : gsl::at(new_indices, i) =
581 : global_corner_index_[i] +
582 : (gsl::at(array_sides_, i) == Side::Upper ? 1 : 0);
583 : }
584 : const Index<VolumeDim> interior_multi_index(new_indices);
585 : return collapsed_index(interior_multi_index, global_corner_extents_);
586 : }
587 :
588 0 : const std::array<Side, VolumeDim>& operator()() const { return array_sides_; }
589 :
590 0 : const std::array<Side, VolumeDim>& operator*() const { return array_sides_; }
591 :
592 0 : const std::array<double, VolumeDim>& coords_of_corner() const {
593 : return coords_of_corner_;
594 : }
595 :
596 0 : const std::array<Direction<VolumeDim>, VolumeDim>& directions_of_corner()
597 : const {
598 : return array_directions_;
599 : }
600 :
601 0 : void setup_from_local_corner_number() {
602 : for (size_t i = 0; i < VolumeDim; i++) {
603 : gsl::at(coords_of_corner_, i) =
604 : 2.0 * get_nth_bit(local_corner_number_, i) - 1.0;
605 : gsl::at(array_sides_, i) =
606 : 2 * get_nth_bit(local_corner_number_, i) - 1 == 1 ? Side::Upper
607 : : Side::Lower;
608 : gsl::at(array_directions_, i) =
609 : Direction<VolumeDim>(i, gsl::at(array_sides_, i));
610 : }
611 : }
612 :
613 : private:
614 0 : size_t local_corner_number_ = 0;
615 0 : size_t global_corner_number_{std::numeric_limits<size_t>::max()};
616 0 : Index<VolumeDim> global_corner_index_{};
617 0 : Index<VolumeDim> global_corner_extents_{};
618 0 : std::array<Side, VolumeDim> array_sides_ = make_array<VolumeDim>(Side::Lower);
619 0 : std::array<Direction<VolumeDim>, VolumeDim> array_directions_{};
620 0 : std::array<double, VolumeDim> coords_of_corner_ = make_array<VolumeDim>(-1.0);
621 : };
622 :
623 : /// \ingroup ComputationalDomainGroup
624 : /// Iterates over the 2^(VolumeDim-1) logical corners of the face of a
625 : /// VolumeDim-dimensional cube in the given direction.
626 : template <size_t VolumeDim>
627 1 : class FaceCornerIterator {
628 : public:
629 0 : explicit FaceCornerIterator(Direction<VolumeDim> direction);
630 :
631 0 : void operator++() {
632 : face_index_++;
633 : do {
634 : index_++;
635 : } while (get_nth_bit(index_, direction_.dimension()) ==
636 : (direction_.side() == Side::Upper ? 0 : 1));
637 : for (size_t i = 0; i < VolumeDim; ++i) {
638 : corner_[i] = 2 * static_cast<int>(get_nth_bit(index_, i)) - 1;
639 : }
640 : }
641 :
642 0 : explicit operator bool() const {
643 : return face_index_ < two_to_the(VolumeDim - 1);
644 : }
645 :
646 0 : tnsr::I<double, VolumeDim, Frame::BlockLogical> operator()() const {
647 : return corner_;
648 : }
649 :
650 0 : tnsr::I<double, VolumeDim, Frame::BlockLogical> operator*() const {
651 : return corner_;
652 : }
653 :
654 : // Returns the value used to construct the logical corner.
655 0 : size_t volume_index() const { return index_; }
656 :
657 : // Returns the number of times operator++ has been called.
658 0 : size_t face_index() const { return face_index_; }
659 :
660 : private:
661 0 : const Direction<VolumeDim> direction_;
662 0 : size_t index_;
663 0 : size_t face_index_ = 0;
664 0 : tnsr::I<double, VolumeDim, Frame::BlockLogical> corner_;
665 : };
666 :
667 : template <size_t VolumeDim>
668 : FaceCornerIterator<VolumeDim>::FaceCornerIterator(
669 : Direction<VolumeDim> direction)
670 : : direction_(std::move(direction)),
671 : index_(direction_.side() == Side::Upper
672 : ? two_to_the(direction_.dimension())
673 : : 0) {
674 : for (size_t i = 0; i < VolumeDim; ++i) {
675 : corner_[i] = 2 * static_cast<int>(get_nth_bit(index_, i)) - 1;
676 : }
677 : }
678 :
679 0 : std::ostream& operator<<(std::ostream& os, const ShellWedges& which_wedges);
680 :
681 : template <>
682 0 : struct Options::create_from_yaml<ShellWedges> {
683 : template <typename Metavariables>
684 0 : static ShellWedges create(const Options::Option& options) {
685 : return create<void>(options);
686 : }
687 : };
688 : template <>
689 0 : ShellWedges Options::create_from_yaml<ShellWedges>::create<void>(
690 : const Options::Option& options);
|
