diff --git a/src/Domain/Block.cpp b/src/Domain/Block.cpp index 1cc92abc78be..5bc37ca41ab5 100644 --- a/src/Domain/Block.cpp +++ b/src/Domain/Block.cpp @@ -27,6 +27,31 @@ Block::Block( std::string name) : stationary_map_(std::move(stationary_map)), id_(id), + geometry_(domain::BlockGeometry::Cube), + name_(std::move(name)) { + for (auto& [direction, neighbor] : neighbors) { + neighbors_[direction].insert(std::move(neighbor)); + } + // Loop over Directions to search which Directions were not set to neighbors_, + // set these Directions to external_boundaries_. + for (const auto& direction : Direction::all_directions()) { + if (neighbors_.find(direction) == neighbors_.end()) { + external_boundaries_.emplace(direction); + } + } +} + +template +Block::Block( + std::unique_ptr>&& stationary_map, + const size_t id, + DirectionMap>> + neighbors, + std::string name) + : stationary_map_(std::move(stationary_map)), + id_(id), + geometry_(domain::BlockGeometry::SphericalShell), neighbors_(std::move(neighbors)), name_(std::move(name)) { // Loop over Directions to search which Directions were not set to neighbors_, @@ -118,7 +143,7 @@ void Block::inject_time_dependent_map( template void Block::pup(PUP::er& p) { - size_t version = 1; + size_t version = 2; p | version; // Remember to increment the version number when making changes to this // function. Retain support for unpacking data written by previous versions @@ -130,7 +155,18 @@ void Block::pup(PUP::er& p) { p | moving_mesh_grid_to_distorted_map_; p | moving_mesh_distorted_to_inertial_map_; p | id_; - p | neighbors_; + if (version < 2) { + geometry_ = domain::BlockGeometry::Cube; + DirectionMap> neighbors; + p | neighbors; + neighbors_.clear(); + for (auto& [direction, neighbor] : neighbors) { + neighbors_[direction].insert(std::move(neighbor)); + } + } else { + p | geometry_; + p | neighbors_; + } p | external_boundaries_; } if (version >= 1) { @@ -141,6 +177,7 @@ void Block::pup(PUP::er& p) { template std::ostream& operator<<(std::ostream& os, const Block& block) { os << "Block " << block.id() << " (" << block.name() << "):\n"; + os << "Geometry: " << block.geometry() << '\n'; os << "Neighbors: " << block.neighbors() << '\n'; os << "External boundaries: " << block.external_boundaries() << '\n'; os << "Is time dependent: " << std::boolalpha << block.is_time_dependent(); @@ -150,7 +187,8 @@ std::ostream& operator<<(std::ostream& os, const Block& block) { template bool operator==(const Block& lhs, const Block& rhs) { bool blocks_are_equal = - (lhs.id() == rhs.id() and lhs.neighbors() == rhs.neighbors() and + (lhs.id() == rhs.id() and lhs.geometry() == rhs.geometry() and + lhs.neighbors() == rhs.neighbors() and lhs.external_boundaries() == rhs.external_boundaries() and lhs.name() == rhs.name() and lhs.is_time_dependent() == rhs.is_time_dependent() and diff --git a/src/Domain/Block.hpp b/src/Domain/Block.hpp index 5ccad0069a04..efab9b2a0690 100644 --- a/src/Domain/Block.hpp +++ b/src/Domain/Block.hpp @@ -13,6 +13,7 @@ #include #include "Domain/CoordinateMaps/CoordinateMap.hpp" +#include "Domain/Structure/BlockGeometry.hpp" #include "Domain/Structure/BlockNeighbor.hpp" #include "Domain/Structure/Direction.hpp" #include "Domain/Structure/DirectionMap.hpp" @@ -33,7 +34,10 @@ class er; /// Elements that cover a region of the computational domain. /// /// Each codimension 1 boundary of a Block is either an external -/// boundary or identical to a boundary of one other Block. +/// boundary or an internal boundary to one or more neighboring blocks. The only +/// currently supported case where blocks can have more than one neighbor are +/// interfaces between a spherical shell and wedges. In all other cases the +/// internal boundaries between neighboring blocks must be identical. /// /// A Block has logical coordinates that go from -1 to +1 in each /// dimension. The global coordinates are obtained from the logical @@ -44,6 +48,8 @@ class er; template class Block { public: + /// Block with one neighbor per direction. Currently always a deformed cube. + /// /// \param stationary_map the CoordinateMap. /// \param id a unique ID. /// \param neighbors info about the Blocks that share a codimension 1 @@ -54,6 +60,15 @@ class Block { size_t id, DirectionMap> neighbors, std::string name = ""); + /// Block with multiple neighbors per direction. Currently only supports + /// spherical shells. + Block(std::unique_ptr>&& stationary_map, + size_t id, + DirectionMap>> + neighbors, + std::string name = ""); + Block() = default; ~Block() = default; Block(const Block&) = delete; @@ -61,6 +76,8 @@ class Block { Block& operator=(const Block&) = delete; Block& operator=(Block&&) = default; + domain::BlockGeometry geometry() const { return geometry_; } + /// \brief The map used when the coordinate map is time-independent. /// /// \see is_time_dependent() @@ -148,7 +165,8 @@ class Block { size_t id() const { return id_; } /// Information about the neighboring Blocks. - const DirectionMap>& neighbors() const { + const DirectionMap>>& + neighbors() const { return neighbors_; } @@ -186,7 +204,9 @@ class Block { moving_mesh_distorted_to_inertial_map_{nullptr}; size_t id_{0}; - DirectionMap> neighbors_; + domain::BlockGeometry geometry_{domain::BlockGeometry::Cube}; + DirectionMap>> + neighbors_; std::unordered_set> external_boundaries_; std::string name_; }; diff --git a/src/Domain/CMakeLists.txt b/src/Domain/CMakeLists.txt index 39833400bf6a..693bd08fc694 100644 --- a/src/Domain/CMakeLists.txt +++ b/src/Domain/CMakeLists.txt @@ -40,6 +40,7 @@ spectre_target_headers( CreateInitialElement.hpp Domain.hpp DomainHelpers.hpp + DomainHelpers.tpp ElementDistribution.hpp ElementLogicalCoordinates.hpp ElementMap.hpp diff --git a/src/Domain/CoordinateMaps/CMakeLists.txt b/src/Domain/CoordinateMaps/CMakeLists.txt index 340cc7d147c2..085dabc7167e 100644 --- a/src/Domain/CoordinateMaps/CMakeLists.txt +++ b/src/Domain/CoordinateMaps/CMakeLists.txt @@ -31,6 +31,7 @@ spectre_target_sources( Interval.cpp KerrHorizonConforming.cpp Rotation.cpp + ShellType.cpp SpecialMobius.cpp SphericalToCartesianPfaffian.cpp UniformCylindricalEndcap.cpp @@ -72,6 +73,7 @@ spectre_target_headers( ProductMaps.hpp ProductMaps.tpp Rotation.hpp + ShellType.hpp SpecialMobius.hpp SphericalToCartesianPfaffian.hpp Tags.hpp diff --git a/src/Domain/CoordinateMaps/ShellType.cpp b/src/Domain/CoordinateMaps/ShellType.cpp new file mode 100644 index 000000000000..7da241165741 --- /dev/null +++ b/src/Domain/CoordinateMaps/ShellType.cpp @@ -0,0 +1,39 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Domain/CoordinateMaps/ShellType.hpp" + +#include +#include + +#include "Options/Options.hpp" +#include "Options/ParseOptions.hpp" +#include "Utilities/ErrorHandling/Error.hpp" + +namespace domain::CoordinateMaps { + +std::ostream& operator<<(std::ostream& os, const ShellType shell_type) { + switch (shell_type) { + case ShellType::Cubed: + return os << "Cubed"; + case ShellType::Spherical: + return os << "Spherical"; + default: + ERROR("Unknown domain::CoordinateMaps::ShellType"); + } +} + +} // namespace domain::CoordinateMaps + +template <> +domain::CoordinateMaps::ShellType +Options::create_from_yaml::create( + const Options::Option& options) { + const auto shell_type = options.parse_as(); + if (shell_type == "Cubed") { + return domain::CoordinateMaps::ShellType::Cubed; + } else if (shell_type == "Spherical") { + return domain::CoordinateMaps::ShellType::Spherical; + } + PARSE_ERROR(options.context(), "ShellType must be 'Cubed' or 'Spherical'."); +} diff --git a/src/Domain/CoordinateMaps/ShellType.hpp b/src/Domain/CoordinateMaps/ShellType.hpp new file mode 100644 index 000000000000..a8f42112b75e --- /dev/null +++ b/src/Domain/CoordinateMaps/ShellType.hpp @@ -0,0 +1,44 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include + +/// \cond +namespace Options { +class Option; +template +struct create_from_yaml; +} // namespace Options +/// \endcond + +namespace domain::CoordinateMaps { + +/*! + * \brief Type of spherical shell: built from four (2D) or six (3D) wedges + * ("cubed") or from a single spherical shell ("spherical"). + * + * Used to select the shell type in the input file. + */ +enum class ShellType { + Cubed, + Spherical, +}; + +std::ostream& operator<<(std::ostream& os, ShellType shell_type); + +} // namespace domain::CoordinateMaps + +template <> +struct Options::create_from_yaml { + template + static domain::CoordinateMaps::ShellType create( + const Options::Option& options) { + return create(options); + } +}; +template <> +domain::CoordinateMaps::ShellType +Options::create_from_yaml::create( + const Options::Option& options); diff --git a/src/Domain/CreateInitialElement.cpp b/src/Domain/CreateInitialElement.cpp index 384b75175f32..c784c82ec823 100644 --- a/src/Domain/CreateInitialElement.cpp +++ b/src/Domain/CreateInitialElement.cpp @@ -12,6 +12,7 @@ #include "Domain/Structure/Direction.hpp" #include "Domain/Structure/Element.hpp" #include "Domain/Structure/ElementId.hpp" +#include "Domain/Structure/InitialElementIds.hpp" #include "Domain/Structure/Neighbors.hpp" #include "Domain/Structure/OrientationMap.hpp" #include "Domain/Structure/SegmentId.hpp" @@ -25,18 +26,82 @@ namespace domain::Initialization { template Element create_initial_element( - const ElementId& element_id, const Block& block, + const ElementId& element_id, + const std::vector>& blocks, const std::vector>& initial_refinement_levels) { + const auto& block = blocks[element_id.block_id()]; const auto& neighbors_of_block = block.neighbors(); const auto segment_ids = element_id.segment_ids(); // Declare two helper lambdas for setting the neighbors of an element const auto compute_element_neighbor_in_other_block = - [&block, &initial_refinement_levels, &neighbors_of_block, &segment_ids, - grid_index = - element_id.grid_index()](const Direction& direction) { - const auto& block_neighbor = neighbors_of_block.at(direction); + [&block, &blocks, &initial_refinement_levels, &neighbors_of_block, + &segment_ids, grid_index = element_id.grid_index()]( + const Direction& direction) { + const auto& block_neighbors = neighbors_of_block.at(direction); + + // Handle spherical shells. We assume that spherical shells have no + // angular h-refinement. That makes the logic here quite simple: + // The single-element shell is the neighbor of all radial neighbors. + if (block.geometry() == domain::BlockGeometry::SphericalShell or + (block_neighbors.size() == 1 and + blocks[block_neighbors.begin()->id()].geometry() == + domain::BlockGeometry::SphericalShell)) { +#ifdef SPECTRE_DEBUG + const auto check_angular_refinement = + [&initial_refinement_levels, + &direction](const Block& test_block) { + if (test_block.geometry() != + domain::BlockGeometry::SphericalShell) { + return; + } + const auto& refinement_levels = + initial_refinement_levels[test_block.id()]; + for (size_t d = 0; d < VolumeDim; ++d) { + if (d == direction.dimension()) { + continue; + } + ASSERT(refinement_levels[d] == 0, + "Spherical shells are assumed here to have no angular " + "refinement in angular directions."); + } + }; + check_angular_refinement(block); +#endif // SPECTRE_DEBUG + std::unordered_set> neighbor_ids; + for (const auto& block_neighbor : block_neighbors) { +#ifdef SPECTRE_DEBUG + check_angular_refinement(blocks[block_neighbor.id()]); +#endif // SPECTRE_DEBUG + const auto& orientation = block_neighbor.orientation(); + const auto direction_from_neighbor = + orientation(direction.opposite()); + const auto& refinement_of_neighbor = + initial_refinement_levels[block_neighbor.id()]; + const size_t index_of_neighbor = + direction_from_neighbor.side() == Side::Lower + ? 0 + : (two_to_the(refinement_of_neighbor[direction_from_neighbor + .dimension()]) - + 1); + for (const auto& neighbor_id : initial_element_ids( + block_neighbor.id(), refinement_of_neighbor, grid_index)) { + if (neighbor_id.segment_id(direction_from_neighbor.dimension()) + .index() == index_of_neighbor) { + neighbor_ids.insert(neighbor_id); + } + } + } + return std::make_pair( + direction, Neighbors( + std::move(neighbor_ids), + // TODO: may have to set this orientation + OrientationMap::create_aligned(), + blocks[block_neighbors.begin()->id()].geometry())); + } + + const auto& block_neighbor = *block_neighbors.begin(); const auto& orientation = block_neighbor.orientation(); const auto direction_in_neighbor = orientation(direction); @@ -131,29 +196,28 @@ Element create_initial_element( } return std::make_pair( direction, - Neighbors(std::move(neighbor_ids), orientation)); + Neighbors(std::move(neighbor_ids), orientation, + blocks[block_neighbor.id()].geometry())); }; - const auto compute_element_neighbor_in_same_block = [&element_id, - &segment_ids]( - const Direction< - VolumeDim>& - direction) { - auto segment_ids_of_neighbor = segment_ids; - auto& perpendicular_segment_id = - gsl::at(segment_ids_of_neighbor, direction.dimension()); - const auto index = perpendicular_segment_id.index(); - perpendicular_segment_id = - SegmentId(perpendicular_segment_id.refinement_level(), - direction.side() == Side::Upper ? index + 1 : index - 1); - return std::make_pair( - direction, - Neighbors( - {{ElementId{element_id.block_id(), - std::move(segment_ids_of_neighbor), - element_id.grid_index()}}}, - OrientationMap::create_aligned())); - }; + const auto compute_element_neighbor_in_same_block = + [&element_id, &segment_ids, + geometry = block.geometry()](const Direction& direction) { + auto segment_ids_of_neighbor = segment_ids; + auto& perpendicular_segment_id = + gsl::at(segment_ids_of_neighbor, direction.dimension()); + const auto index = perpendicular_segment_id.index(); + perpendicular_segment_id = + SegmentId(perpendicular_segment_id.refinement_level(), + direction.side() == Side::Upper ? index + 1 : index - 1); + return std::make_pair( + direction, + Neighbors( + {{ElementId{element_id.block_id(), + std::move(segment_ids_of_neighbor), + element_id.grid_index()}}}, + OrientationMap::create_aligned(), geometry)); + }; typename Element::Neighbors_t neighbors_of_element; for (size_t d = 0; d < VolumeDim; ++d) { @@ -186,10 +250,10 @@ Element create_initial_element( #define DIM(data) BOOST_PP_TUPLE_ELEM(0, data) -#define INSTANTIATE(_, data) \ - template Element \ - domain::Initialization::create_initial_element( \ - const ElementId&, const Block&, \ +#define INSTANTIATE(_, data) \ + template Element \ + domain::Initialization::create_initial_element( \ + const ElementId&, const std::vector>&, \ const std::vector>&); GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3)) diff --git a/src/Domain/CreateInitialElement.hpp b/src/Domain/CreateInitialElement.hpp index e33b0796fac1..d1e3b1f68b86 100644 --- a/src/Domain/CreateInitialElement.hpp +++ b/src/Domain/CreateInitialElement.hpp @@ -29,7 +29,8 @@ namespace Initialization { */ template Element create_initial_element( - const ElementId& element_id, const Block& block, + const ElementId& element_id, + const std::vector>& blocks, const std::vector>& initial_refinement_levels); } // namespace Initialization diff --git a/src/Domain/Creators/BinaryCompactObject.cpp b/src/Domain/Creators/BinaryCompactObject.cpp index cda5a42579d2..edb713cc2bd5 100644 --- a/src/Domain/Creators/BinaryCompactObject.cpp +++ b/src/Domain/Creators/BinaryCompactObject.cpp @@ -36,6 +36,7 @@ #include "Domain/Creators/TimeDependentOptions/BinaryCompactObject.hpp" #include "Domain/Domain.hpp" #include "Domain/DomainHelpers.hpp" +#include "Domain/DomainHelpers.tpp" #include "Domain/ExcisionSphere.hpp" #include "Domain/FunctionsOfTime/FixedSpeedCubic.hpp" #include "Domain/FunctionsOfTime/PiecewisePolynomial.hpp" @@ -459,20 +460,19 @@ Domain<3> BinaryCompactObject::create_domain() const { length_inner_cube_ * 0.5, std::array{{offset_x_coord_a_, 0.0, 0.0}})); Maps maps_center_A = - domain::make_vector_coordinate_map_base( - sph_wedge_coordinate_maps( - object_a.inner_radius, object_a.outer_radius, - inner_sphericity_A, 1.0, use_equiangular_map_, - offset_a_optional, false, {}, object_A_radial_distribution), - translation_A); + spherical_shells_coordinate_maps( + object_a.inner_radius, object_a.outer_radius, inner_sphericity_A, + 1.0, use_equiangular_map_, offset_a_optional, false, {}, + object_A_radial_distribution, + {domain::CoordinateMaps::ShellType::Cubed}, ShellWedges::All, + M_PI_2, translation_A); Maps maps_cube_A = - domain::make_vector_coordinate_map_base( - sph_wedge_coordinate_maps( - object_a.outer_radius, sqrt(3.0) * 0.5 * length_inner_cube_, - 1.0, 0.0, use_equiangular_map_, offset_a_optional), - translation_A); + spherical_shells_coordinate_maps( + object_a.outer_radius, sqrt(3.0) * 0.5 * length_inner_cube_, 1.0, + 0.0, use_equiangular_map_, offset_a_optional, false, {}, + {domain::CoordinateMaps::Distribution::Linear}, + {domain::CoordinateMaps::ShellType::Cubed}, ShellWedges::All, + M_PI_2, translation_A); std::move(maps_center_A.begin(), maps_center_A.end(), std::back_inserter(maps)); std::move(maps_cube_A.begin(), maps_cube_A.end(), std::back_inserter(maps)); @@ -502,20 +502,19 @@ Domain<3> BinaryCompactObject::create_domain() const { length_inner_cube_ * 0.5, std::array{{offset_x_coord_b_, 0.0, 0.0}})); Maps maps_center_B = - domain::make_vector_coordinate_map_base( - sph_wedge_coordinate_maps( - object_b.inner_radius, object_b.outer_radius, - inner_sphericity_B, 1.0, use_equiangular_map_, - offset_b_optional, false, {}, object_B_radial_distribution), - translation_B); + spherical_shells_coordinate_maps( + object_b.inner_radius, object_b.outer_radius, inner_sphericity_B, + 1.0, use_equiangular_map_, offset_b_optional, false, {}, + object_B_radial_distribution, + {domain::CoordinateMaps::ShellType::Cubed}, ShellWedges::All, + M_PI_2, translation_B); Maps maps_cube_B = - domain::make_vector_coordinate_map_base( - sph_wedge_coordinate_maps( - object_b.outer_radius, sqrt(3.0) * 0.5 * length_inner_cube_, - 1.0, 0.0, use_equiangular_map_, offset_b_optional), - translation_B); + spherical_shells_coordinate_maps( + object_b.outer_radius, sqrt(3.0) * 0.5 * length_inner_cube_, 1.0, + 0.0, use_equiangular_map_, offset_b_optional, false, {}, + {domain::CoordinateMaps::Distribution::Linear}, + {domain::CoordinateMaps::ShellType::Cubed}, ShellWedges::All, + M_PI_2, translation_B); std::move(maps_center_B.begin(), maps_center_B.end(), std::back_inserter(maps)); std::move(maps_cube_B.begin(), maps_cube_B.end(), std::back_inserter(maps)); @@ -545,12 +544,11 @@ Domain<3> BinaryCompactObject::create_domain() const { // --- Outer spherical shell (10 blocks) --- Maps maps_outer_shell = - domain::make_vector_coordinate_map_base( - sph_wedge_coordinate_maps(envelope_radius_, outer_radius_, 1.0, 1.0, - use_equiangular_map_, std::nullopt, true, - {}, {radial_distribution_outer_shell_}, - ShellWedges::All, opening_angle_)); + spherical_shells_coordinate_maps( + envelope_radius_, outer_radius_, 1.0, 1.0, use_equiangular_map_, + std::nullopt, true, {}, {radial_distribution_outer_shell_}, + {domain::CoordinateMaps::ShellType::Cubed}, ShellWedges::All, + opening_angle_); std::move(maps_outer_shell.begin(), maps_outer_shell.end(), std::back_inserter(maps)); diff --git a/src/Domain/Creators/Sphere.cpp b/src/Domain/Creators/Sphere.cpp index 668afdc5fb22..4d752e2ca330 100644 --- a/src/Domain/Creators/Sphere.cpp +++ b/src/Domain/Creators/Sphere.cpp @@ -30,6 +30,7 @@ #include "Domain/Creators/TimeDependence/None.hpp" #include "Domain/Domain.hpp" #include "Domain/DomainHelpers.hpp" +#include "Domain/DomainHelpers.tpp" #include "Domain/Structure/BlockNeighbor.hpp" #include "Domain/Structure/Direction.hpp" #include "Domain/Structure/DirectionMap.hpp" @@ -37,6 +38,7 @@ #include "Utilities/ErrorHandling/Assert.hpp" #include "Utilities/GetOutput.hpp" #include "Utilities/MakeArray.hpp" +#include "Utilities/Overloader.hpp" namespace Frame { struct Inertial; @@ -79,7 +81,7 @@ Sphere::Sphere( std::optional equatorial_compression, std::vector radial_partitioning, const typename RadialDistribution::type& radial_distribution, - ShellWedges which_wedges, + const typename ShellType::type& shell_type, ShellWedges which_wedges, std::optional time_dependent_options, std::unique_ptr outer_boundary_condition, @@ -153,22 +155,72 @@ Sphere::Sphere( "Add entries to 'RadialPartitioning' to add outer shells for " "which you can select different radial distributions."); } + shell_types_ = std::visit( + Overloader{ + [&num_shells = num_shells_]( + const domain::CoordinateMaps::ShellType& uniform_shell_type) { + return std::vector( + num_shells, uniform_shell_type); + }, + [](const std::vector& + shell_types) { return shell_types; }}, + shell_type); + if (shell_types_.size() != num_shells_) { + PARSE_ERROR(context, + "Specify a 'ShellType' for every spherical shell. You " + "specified " + << shell_types_.size() << " items, but the domain has " + << num_shells_ << " shells."); + } + if (fill_interior_ and + shell_types_.front() != domain::CoordinateMaps::ShellType::Cubed) { + PARSE_ERROR(context, + "The 'ShellType' must be 'Cubed' for the innermost " + "shell filled with a cube. " + "Add entries to 'RadialPartitioning' to add outer shells for " + "which you can select different shell types."); + } + if (which_wedges_ != ShellWedges::All and + shell_types_ != + std::vector{ + num_shells_, domain::CoordinateMaps::ShellType::Cubed}) { + PARSE_ERROR(context, + "To specify 'ShellWedges' other than 'All', you must specify " + "'ShellType' as 'Cubed' for all shells."); + } // Create grid anchors grid_anchors_["Center"] = tnsr::I{std::array{0.0, 0.0, 0.0}}; // Determine number of blocks - num_blocks_per_shell_ = which_wedges_ == ShellWedges::All ? 6 - : which_wedges_ == ShellWedges::FourOnEquator ? 4 + num_blocks_per_cubed_shell_ = which_wedges_ == ShellWedges::All ? 6 + : which_wedges_ == ShellWedges::FourOnEquator + ? 4 : 1; - num_blocks_ = num_blocks_per_shell_ * num_shells_ + (fill_interior_ ? 1 : 0); + num_blocks_ = alg::accumulate( + shell_types_, 0_st, + [num_blocks_per_cubed_shell = num_blocks_per_cubed_shell_]( + const size_t num_blocks, + const domain::CoordinateMaps::ShellType shell_type) { + if (shell_type == domain::CoordinateMaps::ShellType::Cubed) { + return num_blocks + num_blocks_per_cubed_shell; + } else if (shell_type == domain::CoordinateMaps::ShellType::Spherical) { + return num_blocks + 1; + } else { + ERROR("Unknown ShellType"); + } + }); + if (fill_interior_) { + ++num_blocks_; + } // Create block names and groups static std::array wedge_directions{ "UpperZ", "LowerZ", "UpperY", "LowerY", "UpperX", "LowerX"}; for (size_t shell = 0; shell < num_shells_; ++shell) { std::string shell_prefix = "Shell" + std::to_string(shell); + if (shell_types_[shell] == domain::CoordinateMaps::ShellType::Cubed) { for (size_t direction = which_wedge_index(which_wedges_); direction < 6; ++direction) { const std::string wedge_name = @@ -176,6 +228,11 @@ Sphere::Sphere( block_names_.emplace_back(wedge_name); block_groups_[shell_prefix].insert(wedge_name); block_groups_["Wedges"].insert(wedge_name); + } + } else if (shell_types_[shell] == + domain::CoordinateMaps::ShellType::Spherical) { + block_names_.emplace_back(shell_prefix); + block_groups_["SphericalShells"].insert(shell_prefix); } } if (fill_interior_) { @@ -264,10 +321,9 @@ Sphere::Sphere( } Domain<3> Sphere::create_domain() const { - std::vector> corners = - corners_for_radially_layered_domains(num_shells_, fill_interior_, - {{1, 2, 3, 4, 5, 6, 7, 8}}, - which_wedges_); + std::vector> blocks{}; + blocks.reserve(num_blocks_); + double aspect_ratio = 1.0; size_t index_polar_axis = 2; if (equatorial_compression_.has_value()) { @@ -277,14 +333,48 @@ Domain<3> Sphere::create_domain() const { const domain::CoordinateMaps::EquatorialCompression compression{ aspect_ratio, index_polar_axis}; - auto coord_maps = domain::make_vector_coordinate_map_base( - sph_wedge_coordinate_maps( + auto coord_maps = + spherical_shells_coordinate_maps( inner_radius_, outer_radius_, fill_interior_ ? std::get(interior_).sphericity : 1.0, 1.0, use_equiangular_map_, std::nullopt, false, radial_partitioning_, - radial_distribution_, which_wedges_), - compression); + radial_distribution_, shell_types_, which_wedges_, M_PI, compression); + + // Assuming only the outer shell is a spherical shell + const auto corners_of_wedges = corners_for_radially_layered_domains( + num_shells_ - 1, fill_interior_, {{1, 2, 3, 4, 5, 6, 7, 8}}, + which_wedges_); + std::vector>> neighbors_of_wedges{}; + set_internal_boundaries<3>(make_not_null(&neighbors_of_wedges), + corners_of_wedges); + std::unordered_set> neighbors_of_outer_shell{}; + for (size_t i = 0; i < coord_maps.size() - 1; i++) { + if (fill_interior_ and i < 7) { + // Adjacent to the inner cube + neighbors_of_wedges[i][Direction<3>::lower_zeta()] = BlockNeighbor<3>{ + num_blocks_ - 1, + neighbors_of_wedges[i][Direction<3>::lower_zeta()].orientation()}; + } + if (i >= coord_maps.size() - 7) { + // Adjacent to the outer shell + // Account for different ordering of logical coordinates in wedges [theta, + // phi, r] and spherical shells [r, theta, phi] + const OrientationMap<3> orientation{ + {{Direction<3>::upper_zeta(), Direction<3>::upper_xi(), + Direction<3>::upper_eta()}}, + {{Direction<3>::upper_xi(), Direction<3>::upper_eta(), + Direction<3>::upper_zeta()}}}; + neighbors_of_wedges[i][Direction<3>::upper_zeta()] = + BlockNeighbor<3>{coord_maps.size() - 1, orientation}; + neighbors_of_outer_shell.emplace(i, orientation.inverse_map()); + } + blocks.emplace_back(std::move(coord_maps[i]), i, + std::move(neighbors_of_wedges[i]), block_names_[i]); + } + blocks.emplace_back(Block<3>( + std::move(coord_maps.back()), coord_maps.size() - 1, + {{Direction<3>::lower_xi(), std::move(neighbors_of_outer_shell)}}, + block_names_[coord_maps.size() - 1])); std::unordered_map> excision_spheres{}; @@ -318,6 +408,9 @@ Domain<3> Sphere::create_domain() const { BulgedCube{inner_radius_, inner_cube_sphericity, use_equiangular_map_})); } + blocks.emplace_back(std::move(coord_maps.back()), num_blocks_ - 1, + std::move(neighbors_of_wedges.back()), + block_names_.back()); } else { // Set up excision sphere only for ShellWedges::All // - The first 6 blocks enclose the excised sphere, see @@ -338,14 +431,14 @@ Domain<3> Sphere::create_domain() const { } } - Domain<3> domain{std::move(coord_maps), corners, {}, - std::move(excision_spheres), block_names_, block_groups_}; + Domain<3> domain{std::move(blocks), std::move(excision_spheres), + block_groups_}; ASSERT(domain.blocks().size() == num_blocks_, "Unexpected number of blocks. Expected " << num_blocks_ << " but created " << domain.blocks().size() << "."); - if (time_dependent_options_.has_value()) { + if (time_dependent_options_.has_value()) { // TODO std::vector>> block_maps_grid_to_inertial{num_blocks_}; @@ -362,18 +455,18 @@ Domain<3> Sphere::create_domain() const { time_dependent_options_.value()); const size_t first_block_outer_shell = - (num_shells_ - 1) * num_blocks_per_shell_; + (num_shells_ - 1) * num_blocks_per_cubed_shell_; for (size_t block_id = 0; block_id < num_blocks_; block_id++) { const bool is_outer_shell = block_id >= first_block_outer_shell and - block_id < (first_block_outer_shell + num_blocks_per_shell_); + block_id < (first_block_outer_shell + num_blocks_per_cubed_shell_); const bool is_inner_cube = fill_interior_ and (block_id == num_blocks_ - 1); // Correct for 'which_wedges' option - const size_t shell = block_id / num_blocks_per_shell_; + const size_t shell = block_id / num_blocks_per_cubed_shell_; const size_t block_number = shell * 6 + which_wedge_index(which_wedges_) + - block_id % num_blocks_per_shell_; + block_id % num_blocks_per_cubed_shell_; block_maps_grid_to_distorted[block_id] = hard_coded_options.grid_to_distorted_map(block_number, is_inner_cube); diff --git a/src/Domain/Creators/Sphere.hpp b/src/Domain/Creators/Sphere.hpp index f70a2e2dea5b..0b3d4518e3a9 100644 --- a/src/Domain/Creators/Sphere.hpp +++ b/src/Domain/Creators/Sphere.hpp @@ -15,6 +15,9 @@ #include "Domain/BoundaryConditions/BoundaryCondition.hpp" #include "Domain/BoundaryConditions/GetBoundaryConditionsBase.hpp" #include "Domain/CoordinateMaps/Distribution.hpp" +#include "Domain/CoordinateMaps/Identity.hpp" +#include "Domain/CoordinateMaps/ShellType.hpp" +#include "Domain/CoordinateMaps/SphericalToCartesianPfaffian.hpp" #include "Domain/Creators/DomainCreator.hpp" #include "Domain/Creators/TimeDependence/TimeDependence.hpp" #include "Domain/Creators/TimeDependentOptions/Sphere.hpp" @@ -178,10 +181,17 @@ class Sphere : public DomainCreator<3> { using Equiangular3D = CoordinateMaps::ProductOf3Maps; using BulgedCube = CoordinateMaps::BulgedCube; + using LogicalToSphericalMap = domain::CoordinateMaps::ProductOf3Maps< + domain::CoordinateMaps::Interval, domain::CoordinateMaps::Identity<1>, + domain::CoordinateMaps::Identity<1>>; public: using maps_list = tmpl::append< tmpl::list< + domain::CoordinateMap< + Frame::BlockLogical, Frame::Inertial, LogicalToSphericalMap, + domain::CoordinateMaps::SphericalToCartesianPfaffian, + domain::CoordinateMaps::EquatorialCompression>, // Inner cube domain::CoordinateMap, @@ -311,6 +321,22 @@ class Sphere : public DomainCreator<3> { "partitions."}; }; + struct ShellType { + using type = std::variant>; + static constexpr Options::String help = { + "Select type of each spherical shell. Specify a list of N+1 shell " + "types for N radial partitions, i.e., one for each shell." + "The shell type can be 'Cubed' to create six wedges that compose the " + "shell, or 'Spherical' to create a single spherical shell with " + "spherical harmonics. The 'Spherical' option " + "is not widely supported yet, so it is unlikely to work. " + "If the interior of the sphere is filled with a cube, the innermost " + "shell must have a 'Cubed' shell type. You can also specify just a " + "single shell type (not in a vector) which will use the same shell " + "type for all partitions."}; + }; + struct WhichWedges { using type = ShellWedges; static constexpr Options::String help = { @@ -340,7 +366,7 @@ class Sphere : public DomainCreator<3> { using basic_options = tmpl::list; template @@ -371,6 +397,8 @@ class Sphere : public DomainCreator<3> { std::vector radial_partitioning = {}, const typename RadialDistribution::type& radial_distribution = domain::CoordinateMaps::Distribution::Linear, + const typename ShellType::type& shell_type = + domain::CoordinateMaps::ShellType::Cubed, ShellWedges which_wedges = ShellWedges::All, std::optional time_dependent_options = std::nullopt, std::unique_ptr @@ -428,6 +456,7 @@ class Sphere : public DomainCreator<3> { std::optional equatorial_compression_{}; std::vector radial_partitioning_{}; std::vector radial_distribution_{}; + std::vector shell_types_{}; ShellWedges which_wedges_ = ShellWedges::All; std::optional time_dependent_options_{}; bool use_hard_coded_maps_{false}; @@ -435,7 +464,7 @@ class Sphere : public DomainCreator<3> { outer_boundary_condition_; size_t num_shells_; size_t num_blocks_; - size_t num_blocks_per_shell_; + size_t num_blocks_per_cubed_shell_; std::vector block_names_{}; std::unordered_map> block_groups_{}; diff --git a/src/Domain/Domain.cpp b/src/Domain/Domain.cpp index 1c956f278835..925812b040ee 100644 --- a/src/Domain/Domain.cpp +++ b/src/Domain/Domain.cpp @@ -21,8 +21,14 @@ struct BlockLogical; } // namespace Frame template -Domain::Domain(std::vector> blocks) - : blocks_(std::move(blocks)) {} +Domain::Domain( + std::vector> blocks, + std::unordered_map> excision_spheres, + std::unordered_map> + block_groups) + : blocks_(std::move(blocks)), + excision_spheres_(std::move(excision_spheres)), + block_groups_(std::move(block_groups)) {} template Domain::Domain( diff --git a/src/Domain/Domain.hpp b/src/Domain/Domain.hpp index b2822d62187a..50e48bfc44e5 100644 --- a/src/Domain/Domain.hpp +++ b/src/Domain/Domain.hpp @@ -98,7 +98,12 @@ namespace domain {} template class Domain { public: - explicit Domain(std::vector> blocks); + explicit Domain( + std::vector> blocks, + std::unordered_map> + excision_spheres = {}, + std::unordered_map> + block_groups = {}); /*! * \brief Create a Domain using CoordinateMaps to encode the Orientations. diff --git a/src/Domain/DomainHelpers.cpp b/src/Domain/DomainHelpers.cpp index 8051cf1b281d..e7990043e2c3 100644 --- a/src/Domain/DomainHelpers.cpp +++ b/src/Domain/DomainHelpers.cpp @@ -592,145 +592,6 @@ size_t which_wedge_index(const ShellWedges& which_wedges) { } } -std::vector> sph_wedge_coordinate_maps( - const double inner_radius, const double outer_radius, - const double inner_sphericity, const double outer_sphericity, - const bool use_equiangular_map, - const std::optional>>& - offset_options, - const bool use_half_wedges, const std::vector& radial_partitioning, - const std::vector& - radial_distribution, - const ShellWedges which_wedges, const double opening_angle) { - ASSERT(not use_half_wedges or which_wedges == ShellWedges::All, - "If we are using half wedges we must also be using ShellWedges::All."); - ASSERT(radial_distribution.size() == 1 + radial_partitioning.size(), - "Specify a radial distribution for every spherical shell. You " - "specified " - << radial_distribution.size() << " items, but the domain has " - << 1 + radial_partitioning.size() << " shells."); - - const auto wedge_orientations = orientations_for_sphere_wrappings(); - - using Wedge3DMap = domain::CoordinateMaps::Wedge<3>; - using Halves = Wedge3DMap::WedgeHalves; - std::vector wedges_for_all_layers{}; - - const size_t number_of_layers = 1 + radial_partitioning.size(); - double temp_inner_radius = inner_radius; - double temp_inner_sphericity = inner_sphericity; - if (offset_options.has_value()) { - for (size_t layer_i = 0; layer_i < number_of_layers; layer_i++) { - const auto& radial_distribution_this_layer = - radial_distribution.at(layer_i); - std::optional optional_outer_radius{}; - if (outer_sphericity != 0.0) { - optional_outer_radius = outer_radius; - } else { - optional_outer_radius = std::nullopt; - } - // Generate wedges/half-wedges a layer at a time. - std::vector wedges_for_this_layer{}; - if (not use_half_wedges) { - for (size_t face_j = which_wedge_index(which_wedges); face_j < 6; - face_j++) { - wedges_for_this_layer.emplace_back( - temp_inner_radius, optional_outer_radius, - offset_options.value().first, offset_options.value().second, - gsl::at(wedge_orientations, face_j), use_equiangular_map, - Halves::Both, radial_distribution_this_layer); - } - } else { - for (size_t i = 0; i < 4; i++) { - wedges_for_this_layer.emplace_back( - temp_inner_radius, optional_outer_radius, - offset_options.value().first, offset_options.value().second, - gsl::at(wedge_orientations, i), use_equiangular_map, - Halves::LowerOnly, radial_distribution_this_layer); - wedges_for_this_layer.emplace_back( - temp_inner_radius, optional_outer_radius, - offset_options.value().first, offset_options.value().second, - gsl::at(wedge_orientations, i), use_equiangular_map, - Halves::UpperOnly, radial_distribution_this_layer); - } - wedges_for_this_layer.emplace_back( - temp_inner_radius, optional_outer_radius, - offset_options.value().first, offset_options.value().second, - gsl::at(wedge_orientations, 4), use_equiangular_map, Halves::Both, - radial_distribution_this_layer); - wedges_for_this_layer.emplace_back( - temp_inner_radius, optional_outer_radius, - offset_options.value().first, offset_options.value().second, - gsl::at(wedge_orientations, 5), use_equiangular_map, Halves::Both, - radial_distribution_this_layer); - } - for (const auto& wedge : wedges_for_this_layer) { - wedges_for_all_layers.push_back(wedge); - } - } - } else { - double temp_outer_radius{}; - for (size_t layer_i = 0; layer_i < number_of_layers; layer_i++) { - const auto& radial_distribution_this_layer = - radial_distribution.at(layer_i); - if (layer_i != radial_partitioning.size()) { - temp_outer_radius = radial_partitioning.at(layer_i); - } else { - temp_outer_radius = outer_radius; - } - // Generate wedges/half-wedges a layer at a time. - std::vector wedges_for_this_layer{}; - if (not use_half_wedges) { - for (size_t face_j = which_wedge_index(which_wedges); face_j < 6; - face_j++) { - wedges_for_this_layer.emplace_back( - temp_inner_radius, temp_outer_radius, temp_inner_sphericity, - outer_sphericity, gsl::at(wedge_orientations, face_j), - use_equiangular_map, Halves::Both, radial_distribution_this_layer, - std::array({{M_PI_2, M_PI_2}})); - } - } else { - for (size_t i = 0; i < 4; i++) { - wedges_for_this_layer.emplace_back( - temp_inner_radius, temp_outer_radius, temp_inner_sphericity, - outer_sphericity, gsl::at(wedge_orientations, i), - use_equiangular_map, Halves::LowerOnly, - radial_distribution_this_layer, - std::array({{opening_angle, M_PI_2}})); - wedges_for_this_layer.emplace_back( - temp_inner_radius, temp_outer_radius, temp_inner_sphericity, - outer_sphericity, gsl::at(wedge_orientations, i), - use_equiangular_map, Halves::UpperOnly, - radial_distribution_this_layer, - std::array({{opening_angle, M_PI_2}})); - } - const double endcap_opening_angle = M_PI - opening_angle; - const std::array endcap_opening_angles = { - {endcap_opening_angle, endcap_opening_angle}}; - wedges_for_this_layer.emplace_back( - temp_inner_radius, temp_outer_radius, temp_inner_sphericity, - outer_sphericity, gsl::at(wedge_orientations, 4), - use_equiangular_map, Halves::Both, radial_distribution_this_layer, - endcap_opening_angles, false); - wedges_for_this_layer.emplace_back( - temp_inner_radius, temp_outer_radius, temp_inner_sphericity, - outer_sphericity, gsl::at(wedge_orientations, 5), - use_equiangular_map, Halves::Both, radial_distribution_this_layer, - endcap_opening_angles, false); - } - for (const auto& wedge : wedges_for_this_layer) { - wedges_for_all_layers.push_back(wedge); - } - - if (layer_i != radial_partitioning.size()) { - temp_inner_radius = radial_partitioning.at(layer_i); - temp_inner_sphericity = outer_sphericity; - } - } - } - return wedges_for_all_layers; -} - std::vector frustum_coordinate_maps( const double length_inner_cube, const double length_outer_cube, const bool equiangular_map_at_outer, const bool equiangular_map_at_inner, diff --git a/src/Domain/DomainHelpers.hpp b/src/Domain/DomainHelpers.hpp index 6c9c802d28f9..40d80befd80d 100644 --- a/src/Domain/DomainHelpers.hpp +++ b/src/Domain/DomainHelpers.hpp @@ -16,6 +16,7 @@ #include "DataStructures/Index.hpp" #include "DataStructures/Tensor/Tensor.hpp" #include "Domain/CoordinateMaps/Distribution.hpp" +#include "Domain/CoordinateMaps/ShellType.hpp" #include "Domain/Structure/Direction.hpp" #include "Domain/Structure/Side.hpp" #include "Utilities/ConstantExpressions.hpp" @@ -179,7 +180,10 @@ size_t which_wedge_index(const ShellWedges& which_wedges); * of pi minus this opening angle. This parameter only has an effect if * `use_half_wedges` is set to `true`. */ -std::vector> sph_wedge_coordinate_maps( +template +std::vector< + std::unique_ptr>> +spherical_shells_coordinate_maps( double inner_radius, double outer_radius, double inner_sphericity, double outer_sphericity, bool use_equiangular_map, const std::optional>>& @@ -188,7 +192,10 @@ std::vector> sph_wedge_coordinate_maps( const std::vector& radial_partitioning = {}, const std::vector& radial_distribution = {domain::CoordinateMaps::Distribution::Linear}, - ShellWedges which_wedges = ShellWedges::All, double opening_angle = M_PI_2); + const std::vector& shell_types = + {domain::CoordinateMaps::ShellType::Cubed}, + ShellWedges which_wedges = ShellWedges::All, double opening_angle = M_PI_2, + const AppendMaps&... append_maps); /// \ingroup ComputationalDomainGroup /// These are the ten Frustums used in the DomainCreators for binary compact diff --git a/src/Domain/DomainHelpers.tpp b/src/Domain/DomainHelpers.tpp new file mode 100644 index 000000000000..f38ddce3d2ba --- /dev/null +++ b/src/Domain/DomainHelpers.tpp @@ -0,0 +1,151 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include "Domain/DomainHelpers.hpp" + +#include +#include +#include + +#include "Domain/CoordinateMaps/Affine.hpp" +#include "Domain/CoordinateMaps/CoordinateMap.tpp" +#include "Domain/CoordinateMaps/Distribution.hpp" +#include "Domain/CoordinateMaps/Identity.hpp" +#include "Domain/CoordinateMaps/Interval.hpp" +#include "Domain/CoordinateMaps/ProductMaps.tpp" +#include "Domain/CoordinateMaps/ShellType.hpp" +#include "Domain/CoordinateMaps/SphericalToCartesianPfaffian.hpp" +#include "Domain/CoordinateMaps/Wedge.hpp" +#include "Utilities/ErrorHandling/Assert.hpp" +#include "Utilities/ErrorHandling/Error.hpp" + +template +std::vector< + std::unique_ptr>> +spherical_shells_coordinate_maps( + const double inner_radius, const double outer_radius, + const double inner_sphericity, const double outer_sphericity, + const bool use_equiangular_map, + const std::optional>>& + offset_options, + const bool use_half_wedges, const std::vector& radial_partitioning, + const std::vector& + radial_distribution, + const std::vector& shell_types, + const ShellWedges which_wedges, const double opening_angle, + const AppendMaps&... append_maps) { + ASSERT(not use_half_wedges or which_wedges == ShellWedges::All, + "If we are using half wedges we must also be using ShellWedges::All."); + ASSERT(radial_distribution.size() == 1 + radial_partitioning.size(), + "Specify a radial distribution for every spherical shell. You " + "specified " + << radial_distribution.size() << " items, but the domain has " + << 1 + radial_partitioning.size() << " shells."); + ASSERT(shell_types.size() == 1 + radial_partitioning.size(), + "Specify a type for every spherical shell. You specified " + << shell_types.size() << " items, but the domain has " + << 1 + radial_partitioning.size() << " shells."); + + const auto wedge_orientations = orientations_for_sphere_wrappings(); + + std::vector< + std::unique_ptr>> + maps{}; + + const size_t number_of_layers = 1 + radial_partitioning.size(); + double temp_inner_radius = inner_radius; + double temp_inner_sphericity = inner_sphericity; + if (offset_options.has_value()) { + ERROR("TODO"); + } else { + double temp_outer_radius{}; + for (size_t layer_i = 0; layer_i < number_of_layers; layer_i++) { + // Determine inner and outer radius for this shell + if (layer_i != radial_partitioning.size()) { + temp_outer_radius = radial_partitioning.at(layer_i); + } else { + temp_outer_radius = outer_radius; + } + if (shell_types[layer_i] == domain::CoordinateMaps::ShellType::Cubed) { + // Compose the shell of wedges (deformed cubes) + using Wedge3DMap = domain::CoordinateMaps::Wedge<3>; + using Halves = Wedge3DMap::WedgeHalves; + std::vector wedges{}; + if (not use_half_wedges) { + for (size_t face_j = which_wedge_index(which_wedges); face_j < 6; + face_j++) { + wedges.emplace_back( + temp_inner_radius, temp_outer_radius, temp_inner_sphericity, + outer_sphericity, gsl::at(wedge_orientations, face_j), + use_equiangular_map, Halves::Both, radial_distribution[layer_i], + std::array({{M_PI_2, M_PI_2}})); + } + } else { + for (size_t i = 0; i < 4; i++) { + wedges.emplace_back( + temp_inner_radius, temp_outer_radius, temp_inner_sphericity, + outer_sphericity, gsl::at(wedge_orientations, i), + use_equiangular_map, Halves::LowerOnly, + radial_distribution[layer_i], + std::array({{opening_angle, M_PI_2}})); + wedges.emplace_back( + temp_inner_radius, temp_outer_radius, temp_inner_sphericity, + outer_sphericity, gsl::at(wedge_orientations, i), + use_equiangular_map, Halves::UpperOnly, + radial_distribution[layer_i], + std::array({{opening_angle, M_PI_2}})); + } + const double endcap_opening_angle = M_PI - opening_angle; + const std::array endcap_opening_angles = { + {endcap_opening_angle, endcap_opening_angle}}; + wedges.emplace_back( + temp_inner_radius, temp_outer_radius, temp_inner_sphericity, + outer_sphericity, gsl::at(wedge_orientations, 4), + use_equiangular_map, Halves::Both, radial_distribution[layer_i], + endcap_opening_angles, false); + wedges.emplace_back( + temp_inner_radius, temp_outer_radius, temp_inner_sphericity, + outer_sphericity, gsl::at(wedge_orientations, 5), + use_equiangular_map, Halves::Both, radial_distribution[layer_i], + endcap_opening_angles, false); + } + auto wedge_maps = + domain::make_vector_coordinate_map_base(std::move(wedges), + append_maps...); + maps.insert(maps.end(), std::make_move_iterator(wedge_maps.begin()), + std::make_move_iterator(wedge_maps.end())); + } else if (shell_types[layer_i] == + domain::CoordinateMaps::ShellType::Spherical) { + // Use a single spherical shell (with spherical harmonics in angular + // directions) + ASSERT( + inner_sphericity == 1. and outer_sphericity == 1., + "Spherical shells only support inner and outer sphericity == 1."); + domain::CoordinateMaps::ProductOf3Maps< + domain::CoordinateMaps::Interval, + domain::CoordinateMaps::Identity<1>, + domain::CoordinateMaps::Identity<1>> + logical_to_spherical_map{{-1., 1., temp_inner_radius, outer_radius, + radial_distribution[layer_i], 0.}, + {}, + {}}; + maps.push_back( + domain::make_coordinate_map_base( + std::move(logical_to_spherical_map), + domain::CoordinateMaps::SphericalToCartesianPfaffian{}, + append_maps...)); + } else { + ERROR("Unknown ShellType"); + } + + if (layer_i != radial_partitioning.size()) { + temp_inner_radius = radial_partitioning.at(layer_i); + temp_inner_sphericity = outer_sphericity; + } + } // for layer + } // if offset_options + return maps; +} diff --git a/src/Domain/ElementDistribution.cpp b/src/Domain/ElementDistribution.cpp index 69b39354c4cc..2cc3ddc5e57d 100644 --- a/src/Domain/ElementDistribution.cpp +++ b/src/Domain/ElementDistribution.cpp @@ -57,14 +57,15 @@ namespace { // local time stepping. template double get_num_points_and_grid_spacing_cost( - const ElementId& element_id, const Block& block, + const ElementId& element_id, const std::vector>& blocks, const std::vector>& initial_refinement_levels, const std::vector>& initial_extents, const Spectral::Quadrature quadrature) { + const auto& block = blocks[element_id.block_id()]; Mesh mesh = ::domain::Initialization::create_initial_mesh( - initial_extents, element_id, quadrature); + initial_extents, element_id, block.geometry(), quadrature); Element element = ::domain::Initialization::create_initial_element( - element_id, block, initial_refinement_levels); + element_id, blocks, initial_refinement_levels); ElementMap element_map{element_id, block}; const tnsr::I logical_coords = logical_coordinates(mesh); @@ -122,7 +123,7 @@ std::unordered_map, double> get_element_costs( element_costs.insert( {element_id, get_num_points_and_grid_spacing_cost( - element_id, block, initial_refinement_levels, + element_id, blocks, initial_refinement_levels, initial_extents, quadrature.value())}); } } @@ -323,7 +324,7 @@ size_t BlockZCurveProcDistribution::get_proc_for_element( template class BlockZCurveProcDistribution; \ double get_num_points_and_grid_spacing_cost( \ const ElementId& element_id, \ - const Block& block, \ + const std::vector>& blocks, \ const std::vector>& \ initial_refinement_levels, \ const std::vector>& initial_extents, \ diff --git a/src/Domain/Structure/BlockGeometry.cpp b/src/Domain/Structure/BlockGeometry.cpp new file mode 100644 index 000000000000..97303c991862 --- /dev/null +++ b/src/Domain/Structure/BlockGeometry.cpp @@ -0,0 +1,24 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Domain/Structure/BlockGeometry.hpp" + +#include +#include + +#include "Utilities/ErrorHandling/Error.hpp" + +namespace domain { + +std::ostream& operator<<(std::ostream& os, const BlockGeometry shell_type) { + switch (shell_type) { + case BlockGeometry::Cube: + return os << "Cube"; + case BlockGeometry::SphericalShell: + return os << "SphericalShell"; + default: + ERROR("Unknown domain::BlockGeometry"); + } +} + +} // namespace domain diff --git a/src/Domain/Structure/BlockGeometry.hpp b/src/Domain/Structure/BlockGeometry.hpp new file mode 100644 index 000000000000..9abd5c49e323 --- /dev/null +++ b/src/Domain/Structure/BlockGeometry.hpp @@ -0,0 +1,25 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +namespace domain { + +/*! + * \ingroup ComputationalDomainGroup + * \brief Geometry of a block in the computational domain + */ +enum class BlockGeometry { + /// A logical cube that can be deformed by coordinate maps. In each direction + /// it has zero or one block neighbor. + Cube, + /// A spherical shell. It only has block neighbors in the radial direction. In + /// each radial direction it has either zero block neighbors (at the + /// boundary), one block neighbor (another spherical shell), or multiple block + /// neighbors (cubes deformed to wedges, 4 in 2D or 6 in 3D). + SphericalShell +}; + +std::ostream& operator<<(std::ostream& os, BlockGeometry shell_type); + +} // namespace domain diff --git a/src/Domain/Structure/BlockNeighbor.hpp b/src/Domain/Structure/BlockNeighbor.hpp index 05efc508d7bb..b8e1b588e3c6 100644 --- a/src/Domain/Structure/BlockNeighbor.hpp +++ b/src/Domain/Structure/BlockNeighbor.hpp @@ -54,6 +54,22 @@ class BlockNeighbor { OrientationMap orientation_; }; +template +size_t hash_value(const BlockNeighbor& block_neighbor) { + // Hashing only the Block ID because we assume it uniquely identifies the + // block neighbor; the orientation map is just extra information. + return std::hash{}(block_neighbor.id()); +} + +namespace std { +template +struct hash> { + size_t operator()(const BlockNeighbor& block_neighbor) const { + return hash_value(block_neighbor); + } +}; +} // namespace std + /// Output operator for BlockNeighbor. template std::ostream& operator<<(std::ostream& os, diff --git a/src/Domain/Structure/CMakeLists.txt b/src/Domain/Structure/CMakeLists.txt index db48be9e2fef..cd1b448383b5 100644 --- a/src/Domain/Structure/CMakeLists.txt +++ b/src/Domain/Structure/CMakeLists.txt @@ -8,6 +8,7 @@ add_spectre_library(${LIBRARY}) spectre_target_sources( ${LIBRARY} PRIVATE + BlockGeometry.cpp BlockGroups.cpp BlockNeighbor.cpp ChildSize.cpp @@ -31,6 +32,7 @@ spectre_target_headers( ${LIBRARY} INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src HEADERS + BlockGeometry.hpp BlockGroups.hpp BlockId.hpp BlockNeighbor.hpp diff --git a/src/Domain/Structure/CreateInitialMesh.cpp b/src/Domain/Structure/CreateInitialMesh.cpp index bf2c6c53d6e9..2335b8d4aec9 100644 --- a/src/Domain/Structure/CreateInitialMesh.cpp +++ b/src/Domain/Structure/CreateInitialMesh.cpp @@ -20,25 +20,38 @@ namespace domain::Initialization { template Mesh create_initial_mesh( const std::vector>& initial_extents, - const ElementId& element_id, const Spectral::Quadrature quadrature, + const ElementId& element_id, const domain::BlockGeometry geometry, + const Spectral::Quadrature quadrature, const OrientationMap& orientation) { const auto& unoriented_extents = initial_extents[element_id.block_id()]; Index extents; for (size_t i = 0; i < Dim; ++i) { extents[i] = gsl::at(unoriented_extents, orientation(i)); } - return {extents.indices(), Spectral::Basis::Legendre, quadrature}; + if (geometry == domain::BlockGeometry::Cube) { + return {extents.indices(), Spectral::Basis::Legendre, quadrature}; + } else { + if constexpr (Dim == 3) { + return {extents.indices(), + {{Spectral::Basis::Legendre, Spectral::Basis::SphericalHarmonic, + Spectral::Basis::SphericalHarmonic}}, + {{quadrature, Spectral::Quadrature::Gauss, + Spectral::Quadrature::Equiangular}}}; + } else { + ERROR("not implemented"); + } + } } } // namespace domain::Initialization #define DIM(data) BOOST_PP_TUPLE_ELEM(0, data) -#define INSTANTIATE(_, data) \ - template Mesh \ - domain::Initialization::create_initial_mesh( \ - const std::vector>&, \ - const ElementId&, const Spectral::Quadrature quadrature, \ - const OrientationMap&); +#define INSTANTIATE(_, data) \ + template Mesh \ + domain::Initialization::create_initial_mesh( \ + const std::vector>&, \ + const ElementId&, domain::BlockGeometry geometry, \ + Spectral::Quadrature quadrature, const OrientationMap&); GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3)) diff --git a/src/Domain/Structure/CreateInitialMesh.hpp b/src/Domain/Structure/CreateInitialMesh.hpp index 485fbe379e74..a02848cd9407 100644 --- a/src/Domain/Structure/CreateInitialMesh.hpp +++ b/src/Domain/Structure/CreateInitialMesh.hpp @@ -8,12 +8,14 @@ #include #include +#include "Domain/Structure/BlockGeometry.hpp" +#include "NumericalAlgorithms/Spectral/Mesh.hpp" +#include "NumericalAlgorithms/Spectral/Quadrature.hpp" + /// \cond template struct ElementId; template -class Mesh; -template struct OrientationMap; namespace Spectral { enum class Quadrature : uint8_t; @@ -36,7 +38,9 @@ namespace domain::Initialization { template Mesh create_initial_mesh( const std::vector>& initial_extents, - const ElementId& element_id, Spectral::Quadrature quadrature, + const ElementId& element_id, + domain::BlockGeometry geometry = domain::BlockGeometry::Cube, + Spectral::Quadrature quadrature = Spectral::Quadrature::GaussLobatto, const OrientationMap& orientation = OrientationMap::create_aligned()); } // namespace domain::Initialization diff --git a/src/Domain/Structure/Neighbors.cpp b/src/Domain/Structure/Neighbors.cpp index c81960df484b..8ee54ceed007 100644 --- a/src/Domain/Structure/Neighbors.cpp +++ b/src/Domain/Structure/Neighbors.cpp @@ -15,14 +15,16 @@ template Neighbors::Neighbors(std::unordered_set> ids, - OrientationMap orientation) - : ids_(std::move(ids)), orientation_(std::move(orientation)) { - // Assuming a maximum 2-to-1 refinement between neighboring elements: - ASSERT(ids_.size() <= maximum_number_of_neighbors_per_direction(VolumeDim), - "Can't have " << ids_.size() << " neighbors in " << VolumeDim - << " dimensions"); + OrientationMap orientation, + domain::BlockGeometry geometry) + : ids_(std::move(ids)), + orientation_(std::move(orientation)), + geometry_(std::move(geometry)) { ASSERT(orientation_ != OrientationMap{}, "Cannot use a default-constructed OrientationMap in Neighbors."); + if (geometry_ == domain::BlockGeometry::SphericalShell) { + ASSERT(ids_.size() == 1, "Only one spherical shell neighbor is possible."); + } } template @@ -31,10 +33,9 @@ void Neighbors::add_ids( for (const auto& id : additional_ids) { ids_.insert(id); } - // Assuming a maximum 2-to-1 refinement between neighboring elements: - ASSERT(ids_.size() <= maximum_number_of_neighbors_per_direction(VolumeDim), - "Can't have " << ids_.size() << " neighbors in " << VolumeDim - << " dimensions"); + if (geometry_ == domain::BlockGeometry::SphericalShell) { + ASSERT(ids_.size() == 1, "Only one spherical shell neighbor is possible."); + } } template @@ -46,7 +47,8 @@ std::ostream& operator<<(std::ostream& os, const Neighbors& n) { template bool operator==(const Neighbors& lhs, const Neighbors& rhs) { - return (lhs.ids() == rhs.ids() and lhs.orientation() == rhs.orientation()); + return (lhs.ids() == rhs.ids() and lhs.orientation() == rhs.orientation() and + lhs.geometry() == rhs.geometry()); } template @@ -59,6 +61,7 @@ template void Neighbors::pup(PUP::er& p) { p | ids_; p | orientation_; + p | geometry_; } #define GET_DIM(data) BOOST_PP_TUPLE_ELEM(0, data) diff --git a/src/Domain/Structure/Neighbors.hpp b/src/Domain/Structure/Neighbors.hpp index a089a8abf9ea..73d478a0a304 100644 --- a/src/Domain/Structure/Neighbors.hpp +++ b/src/Domain/Structure/Neighbors.hpp @@ -10,6 +10,7 @@ #include #include +#include "Domain/Structure/BlockGeometry.hpp" #include "Domain/Structure/OrientationMap.hpp" /// \cond @@ -34,8 +35,11 @@ class Neighbors { /// \param orientation This OrientationMap takes objects in the logical /// coordinate frame of the host Element and maps them to the logical /// coordinate frame of the neighbor Element. + /// \param geometry the geometry of the neighboring block. Neighbors(std::unordered_set> ids, - OrientationMap orientation); + OrientationMap orientation, + // todo: remove default + domain::BlockGeometry geometry = domain::BlockGeometry::Cube); /// Default constructor for Charm++ serialization. Neighbors() = default; @@ -49,6 +53,8 @@ class Neighbors { const OrientationMap& orientation() const { return orientation_; } + const domain::BlockGeometry& geometry() const { return geometry_; } + /// Reset the ids of the neighbors. void set_ids_to(const std::unordered_set> new_ids) { ids_ = std::move(new_ids); @@ -96,6 +102,7 @@ class Neighbors { private: std::unordered_set> ids_; OrientationMap orientation_; + domain::BlockGeometry geometry_; }; /// Output operator for Neighborss. diff --git a/src/Elliptic/DiscontinuousGalerkin/Initialization.cpp b/src/Elliptic/DiscontinuousGalerkin/Initialization.cpp index 762fdce403bf..1859e0ea79e3 100644 --- a/src/Elliptic/DiscontinuousGalerkin/Initialization.cpp +++ b/src/Elliptic/DiscontinuousGalerkin/Initialization.cpp @@ -97,18 +97,18 @@ void InitializeGeometry::apply( const Domain& domain, const domain::FunctionsOfTimeMap& functions_of_time, const Spectral::Quadrature quadrature, const ElementId& element_id) { + const auto& block = domain.blocks()[element_id.block_id()]; // Mesh ASSERT(quadrature == Spectral::Quadrature::GaussLobatto or quadrature == Spectral::Quadrature::Gauss, "The elliptic DG scheme supports Gauss and Gauss-Lobatto " "grids, but the chosen quadrature is: " << quadrature); - *mesh = domain::Initialization::create_initial_mesh(initial_extents, - element_id, quadrature); + *mesh = domain::Initialization::create_initial_mesh( + initial_extents, element_id, block.geometry(), quadrature); // Element - const auto& block = domain.blocks()[element_id.block_id()]; - *element = domain::Initialization::create_initial_element(element_id, block, - initial_refinement); + *element = domain::Initialization::create_initial_element( + element_id, domain.blocks(), initial_refinement); // Neighbor meshes for (const auto& [direction, neighbors] : element->neighbors()) { for (const auto& neighbor_id : neighbors) { diff --git a/src/Evolution/DiscontinuousGalerkin/Initialization/Mortars.cpp b/src/Evolution/DiscontinuousGalerkin/Initialization/Mortars.cpp index bcd041d31c0b..376f0b3f043c 100644 --- a/src/Evolution/DiscontinuousGalerkin/Initialization/Mortars.cpp +++ b/src/Evolution/DiscontinuousGalerkin/Initialization/Mortars.cpp @@ -56,8 +56,8 @@ mortars_apply_impl(const std::vector>& initial_extents, mortar_id, ::dg::mortar_mesh(volume_mesh.slice_away(direction.dimension()), ::domain::Initialization::create_initial_mesh( - initial_extents, neighbor, quadrature, - neighbors.orientation()) + initial_extents, neighbor, neighbors.geometry(), + quadrature, neighbors.orientation()) .slice_away(direction.dimension()))); mortar_sizes.emplace( mortar_id, diff --git a/src/Evolution/Initialization/DgDomain.hpp b/src/Evolution/Initialization/DgDomain.hpp index a26563410a46..49cde9fae000 100644 --- a/src/Evolution/Initialization/DgDomain.hpp +++ b/src/Evolution/Initialization/DgDomain.hpp @@ -145,9 +145,9 @@ struct Domain { const ElementId& element_id) { const auto& my_block = domain.blocks()[element_id.block_id()]; *mesh = ::domain::Initialization::create_initial_mesh( - initial_extents, element_id, quadrature); + initial_extents, element_id, my_block.geometry(), quadrature); *element = ::domain::Initialization::create_initial_element( - element_id, my_block, initial_refinement); + element_id, domain.blocks(), initial_refinement); *element_map = ElementMap{element_id, my_block}; if (my_block.is_time_dependent()) { diff --git a/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/InitializeElementFacesGridCoordinates.hpp b/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/InitializeElementFacesGridCoordinates.hpp index 4d28083f1778..d200253f5c4a 100644 --- a/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/InitializeElementFacesGridCoordinates.hpp +++ b/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/InitializeElementFacesGridCoordinates.hpp @@ -82,10 +82,11 @@ struct InitializeElementFacesGridCoordinates { for (const auto element_id : element_ids) { const auto direction = excision_sphere.abutting_direction(element_id); if (direction.has_value()) { + const auto& current_block = blocks.at(block_id); const auto mesh = ::domain::Initialization::create_initial_mesh( - initial_extents, element_id, quadrature); + initial_extents, element_id, current_block.geometry(), + quadrature); const auto face_mesh = mesh.slice_away(direction.value().dimension()); - const auto& current_block = blocks.at(block_id); const ElementMap element_map{ element_id, current_block.is_time_dependent() diff --git a/src/Executables/Examples/RandomAmr/InitializeDomain.hpp b/src/Executables/Examples/RandomAmr/InitializeDomain.hpp index e0bc500da8f1..5f19779da15b 100644 --- a/src/Executables/Examples/RandomAmr/InitializeDomain.hpp +++ b/src/Executables/Examples/RandomAmr/InitializeDomain.hpp @@ -63,11 +63,11 @@ struct Domain { const std::vector>& initial_refinement, const Spectral::Quadrature& quadrature, const ElementId& element_id) { - const auto& my_block = domain.blocks()[element_id.block_id()]; + const auto& block = domain.blocks()[element_id.block_id()]; *mesh = ::domain::Initialization::create_initial_mesh( - initial_extents, element_id, quadrature); + initial_extents, element_id, block.geometry(), quadrature); *element = ::domain::Initialization::create_initial_element( - element_id, my_block, initial_refinement); + element_id, domain.blocks(), initial_refinement); } }; } // namespace amr::Initialization diff --git a/src/Executables/ExportCoordinates/ExportCoordinates.hpp b/src/Executables/ExportCoordinates/ExportCoordinates.hpp index d7853a1e5aae..d2e8f3e3aca3 100644 --- a/src/Executables/ExportCoordinates/ExportCoordinates.hpp +++ b/src/Executables/ExportCoordinates/ExportCoordinates.hpp @@ -130,13 +130,13 @@ struct ExportCoordinates { const double time = get(box); const auto& mesh = get>(box); - const auto& inv_jacobian = - db::get>(box); + // const auto& inv_jacobian = + // db::get>(box); const auto& inertial_coordinates = db::get>(box); - const auto deriv_inertial_coordinates = - partial_derivative(inertial_coordinates, mesh, inv_jacobian); + // const auto deriv_inertial_coordinates = + // partial_derivative(inertial_coordinates, mesh, inv_jacobian); // Collect volume data // Remove tensor types, only storing individual components std::vector components; @@ -148,13 +148,13 @@ struct ExportCoordinates { inertial_coordinates.get(d)); } - for (size_t i = 0; i < deriv_inertial_coordinates.size(); ++i) { - components.emplace_back( - "DerivInertialCoordinates_" + - deriv_inertial_coordinates.component_name( - deriv_inertial_coordinates.get_tensor_index(i)), - deriv_inertial_coordinates[i]); - } + // for (size_t i = 0; i < deriv_inertial_coordinates.size(); ++i) { + // components.emplace_back( + // "DerivInertialCoordinates_" + + // deriv_inertial_coordinates.component_name( + // deriv_inertial_coordinates.get_tensor_index(i)), + // deriv_inertial_coordinates[i]); + // } // Also output the determinant of the inverse jacobian, which measures // the expansion and compression of the grid @@ -169,27 +169,27 @@ struct ExportCoordinates { // Also output the jacobian diagnostic, which compares the analytic // Jacobian (via the CoordinateMap) to the numerical Jacobian // (computed via logical_partial_derivative) - const auto& jacobian = determinant_and_inverse(inv_jacobian).second; - tnsr::i jac_diag{ - mesh.number_of_grid_points(), 0.0}; - domain::jacobian_diagnostic(make_not_null(&jac_diag), jacobian, - inertial_coordinates, mesh); - for (size_t i = 0; i < Dim; ++i) { - components.emplace_back( - "JacobianDiagnostic_" + - jac_diag.component_name(jac_diag.get_tensor_index(i)), - jac_diag.get(i)); - } + // const auto& jacobian = determinant_and_inverse(inv_jacobian).second; + // tnsr::i jac_diag{ + // mesh.number_of_grid_points(), 0.0}; + // domain::jacobian_diagnostic(make_not_null(&jac_diag), jacobian, + // inertial_coordinates, mesh); + // for (size_t i = 0; i < Dim; ++i) { + // components.emplace_back( + // "JacobianDiagnostic_" + + // jac_diag.component_name(jac_diag.get_tensor_index(i)), + // jac_diag.get(i)); + // } // Also output the computation domain metric - const auto& flat_logical_metric = - db::get>(box); - for (size_t i = 0; i < flat_logical_metric.size(); ++i) { - components.emplace_back( - db::tag_name>() + - flat_logical_metric.component_suffix(i), - flat_logical_metric[i]); - } + // const auto& flat_logical_metric = + // db::get>(box); + // for (size_t i = 0; i < flat_logical_metric.size(); ++i) { + // components.emplace_back( + // db::tag_name>() + + // flat_logical_metric.component_suffix(i), + // flat_logical_metric[i]); + // } // Send data to volume observer auto& local_observer = *Parallel::local_branch( diff --git a/tests/InputFiles/ExportCoordinates/Input3D.yaml b/tests/InputFiles/ExportCoordinates/Input3D.yaml index 919a4a658c49..0d83fa565aad 100644 --- a/tests/InputFiles/ExportCoordinates/Input3D.yaml +++ b/tests/InputFiles/ExportCoordinates/Input3D.yaml @@ -30,17 +30,19 @@ ResourceInfo: Singletons: Auto DomainCreator: - Brick: - LowerBound: [-0.5, -0.5, -0.5] - UpperBound: [0.5, 0.5, 0.5] - Distribution: [Linear, Linear, Linear] - IsPeriodicIn: [false, false, false] - InitialRefinement: [0, 0, 0] - InitialGridPoints: [3, 3, 3] - TimeDependence: - UniformTranslation: - InitialTime: 0.0 - Velocity: [0.5, 0.0, 0.0] + Sphere: + InnerRadius: 2. + OuterRadius: 10. + Interior: Excise + InitialRefinement: 0 + InitialGridPoints: 12 + UseEquiangularMap: True + EquatorialCompression: None + WhichWedges: All + RadialPartitioning: [4.] + RadialDistribution: [Linear, Logarithmic] + ShellType: [Cubed, Spherical] + TimeDependentMaps: None SpatialDiscretization: ActiveGrid: Dg diff --git a/tests/Unit/Domain/Structure/CMakeLists.txt b/tests/Unit/Domain/Structure/CMakeLists.txt index c92226c16520..6eac448f8ab5 100644 --- a/tests/Unit/Domain/Structure/CMakeLists.txt +++ b/tests/Unit/Domain/Structure/CMakeLists.txt @@ -4,6 +4,7 @@ set(LIBRARY "Test_DomainStructure") set(LIBRARY_SOURCES + Test_BlockGeometry.cpp Test_BlockGroups.cpp Test_BlockId.cpp Test_BlockNeighbor.cpp diff --git a/tests/Unit/Domain/Structure/Test_BlockGeometry.cpp b/tests/Unit/Domain/Structure/Test_BlockGeometry.cpp new file mode 100644 index 000000000000..1ffdaf799e68 --- /dev/null +++ b/tests/Unit/Domain/Structure/Test_BlockGeometry.cpp @@ -0,0 +1,18 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Framework/TestingFramework.hpp" + +#include + +#include "Domain/Structure/BlockGeometry.hpp" +#include "Utilities/GetOutput.hpp" + +namespace domain { + +SPECTRE_TEST_CASE("Unit.Domain.Structure.BlockGeometry", "[Domain][Unit]") { + CHECK(get_output(domain::BlockGeometry::Cube) == "Cube"); + CHECK(get_output(domain::BlockGeometry::SphericalShell) == "SphericalShell"); +} + +} // namespace domain diff --git a/tests/Unit/Domain/Structure/Test_BlockNeighbor.cpp b/tests/Unit/Domain/Structure/Test_BlockNeighbor.cpp index 5caea2fd66cf..e7235fbd1935 100644 --- a/tests/Unit/Domain/Structure/Test_BlockNeighbor.cpp +++ b/tests/Unit/Domain/Structure/Test_BlockNeighbor.cpp @@ -43,6 +43,11 @@ SPECTRE_TEST_CASE("Unit.Domain.Structure.BlockNeighbor", "[Domain][Unit]") { // Test serialization: test_serialization(custom_neighbor); + // Test hash: + CHECK(hash_value(custom_neighbor) == + hash_value(serialize_and_deserialize(custom_neighbor))); + CHECK(hash_value(custom_neighbor) != hash_value(BlockNeighbor<3>{1, {}})); + // Test semantics: const auto custom_copy = custom_neighbor; test_copy_semantics(custom_neighbor); diff --git a/tests/Unit/Domain/Test_CreateInitialElement.cpp b/tests/Unit/Domain/Test_CreateInitialElement.cpp index d83529b8806a..25f87226f629 100644 --- a/tests/Unit/Domain/Test_CreateInitialElement.cpp +++ b/tests/Unit/Domain/Test_CreateInitialElement.cpp @@ -34,11 +34,11 @@ namespace { void test_create_initial_element( - const ElementId<2>& element_id, const Block<2>& block, + const ElementId<2>& element_id, const std::vector>& blocks, const std::vector>& refinement_levels, const DirectionMap<2, Neighbors<2>>& expected_neighbors) { const auto created_element = domain::Initialization::create_initial_element( - element_id, block, refinement_levels); + element_id, blocks, refinement_levels); const Element<2> expected_element{element_id, expected_neighbors}; CHECK(created_element == expected_element); } @@ -52,16 +52,17 @@ void test_h_refinement() { const std::unordered_set>& expected_neighbors) { CAPTURE(neighbor_orientation); CAPTURE(neighbor_refinement); - const Block<3> self_block( - domain::make_coordinate_map_base( - domain::CoordinateMaps::Identity<3>{}), - 0, {{neighbor_direction, {1, neighbor_orientation}}}); + std::vector> blocks{}; + blocks.push_back({domain::make_coordinate_map_base( + domain::CoordinateMaps::Identity<3>{}), + 0, + {{neighbor_direction, {1, neighbor_orientation}}}}); const std::vector> refinement_levels{ {{1, 1, 1}}, neighbor_refinement}; const auto refined_neighbors = - domain::Initialization::create_initial_element(self_id, self_block, + domain::Initialization::create_initial_element(self_id, blocks, refinement_levels) .neighbors() .at(neighbor_direction) @@ -398,18 +399,18 @@ SPECTRE_TEST_CASE("Unit.Domain.CreateInitialElement", "[Domain][Unit]") { make_array(Direction<2>::upper_xi(), Direction<2>::upper_eta())); OrientationMap<2> unaligned( make_array(Direction<2>::lower_eta(), Direction<2>::upper_xi())); - Block<2> test_block( - domain::make_coordinate_map_base( - domain::CoordinateMaps::Identity<2>{}), - 0, - {{Direction<2>::upper_xi(), BlockNeighbor<2>{1, aligned}}, - {Direction<2>::upper_eta(), BlockNeighbor<2>{2, unaligned}}}); + std::vector> blocks{}; + blocks.push_back( + {domain::make_coordinate_map_base( + domain::CoordinateMaps::Identity<2>{}), + 0, + {{Direction<2>::upper_xi(), BlockNeighbor<2>{1, aligned}}, + {Direction<2>::upper_eta(), BlockNeighbor<2>{2, unaligned}}}}); std::vector> refinement{{{2, 3}}, {{2, 3}}, {{3, 2}}}; // interior element test_create_initial_element( - ElementId<2>{0, {{SegmentId{2, 2}, SegmentId{3, 4}}}}, test_block, - refinement, + ElementId<2>{0, {{SegmentId{2, 2}, SegmentId{3, 4}}}}, blocks, refinement, {{Direction<2>::upper_xi(), Neighbors<2>{{ElementId<2>{0, {{SegmentId{2, 3}, SegmentId{3, 4}}}}}, aligned}}, @@ -425,8 +426,7 @@ SPECTRE_TEST_CASE("Unit.Domain.CreateInitialElement", "[Domain][Unit]") { // element on external boundary test_create_initial_element( - ElementId<2>{0, {{SegmentId{2, 0}, SegmentId{3, 0}}}}, test_block, - refinement, + ElementId<2>{0, {{SegmentId{2, 0}, SegmentId{3, 0}}}}, blocks, refinement, {{Direction<2>::upper_xi(), Neighbors<2>{{ElementId<2>{0, {{SegmentId{2, 1}, SegmentId{3, 0}}}}}, aligned}}, @@ -436,8 +436,7 @@ SPECTRE_TEST_CASE("Unit.Domain.CreateInitialElement", "[Domain][Unit]") { // element bounding aligned neighbor block test_create_initial_element( - ElementId<2>{0, {{SegmentId{2, 3}, SegmentId{3, 4}}}}, test_block, - refinement, + ElementId<2>{0, {{SegmentId{2, 3}, SegmentId{3, 4}}}}, blocks, refinement, {{Direction<2>::upper_xi(), Neighbors<2>{{ElementId<2>{1, {{SegmentId{2, 0}, SegmentId{3, 4}}}}}, aligned}}, @@ -453,8 +452,7 @@ SPECTRE_TEST_CASE("Unit.Domain.CreateInitialElement", "[Domain][Unit]") { // element bounding unaligned neighbor block test_create_initial_element( - ElementId<2>{0, {{SegmentId{2, 2}, SegmentId{3, 7}}}}, test_block, - refinement, + ElementId<2>{0, {{SegmentId{2, 2}, SegmentId{3, 7}}}}, blocks, refinement, {{Direction<2>::upper_xi(), Neighbors<2>{{ElementId<2>{0, {{SegmentId{2, 3}, SegmentId{3, 7}}}}}, aligned}}, @@ -470,8 +468,7 @@ SPECTRE_TEST_CASE("Unit.Domain.CreateInitialElement", "[Domain][Unit]") { // element bounding both neighbor blocks test_create_initial_element( - ElementId<2>{0, {{SegmentId{2, 3}, SegmentId{3, 7}}}}, test_block, - refinement, + ElementId<2>{0, {{SegmentId{2, 3}, SegmentId{3, 7}}}}, blocks, refinement, {{Direction<2>::upper_xi(), Neighbors<2>{{ElementId<2>{1, {{SegmentId{2, 0}, SegmentId{3, 7}}}}}, aligned}}, @@ -490,7 +487,7 @@ SPECTRE_TEST_CASE("Unit.Domain.CreateInitialElement", "[Domain][Unit]") { const size_t grid_index = 3; test_create_initial_element( ElementId<2>{0, {{SegmentId{2, 2}, SegmentId{3, 4}}}, grid_index}, - test_block, refinement, + blocks, refinement, {{Direction<2>::upper_xi(), Neighbors<2>{ {ElementId<2>{ diff --git a/tests/Unit/Elliptic/SubdomainPreconditioners/Test_MinusLaplacian.cpp b/tests/Unit/Elliptic/SubdomainPreconditioners/Test_MinusLaplacian.cpp index ad5267684801..c1a1d44b2536 100644 --- a/tests/Unit/Elliptic/SubdomainPreconditioners/Test_MinusLaplacian.cpp +++ b/tests/Unit/Elliptic/SubdomainPreconditioners/Test_MinusLaplacian.cpp @@ -238,11 +238,11 @@ auto make_databox_with_boundary_conditions() { const ElementId right_element_id{0, {{{2, 1}, {0, 0}}}}; const ElementId bottom_element_id{1, {{{1, 0}, {1, 1}}}}; Element central_element = domain::Initialization::create_initial_element( - central_element_id, domain.blocks()[0], refinement); + central_element_id, domain.blocks(), refinement); Element right_element = domain::Initialization::create_initial_element( - right_element_id, domain.blocks()[0], refinement); + right_element_id, domain.blocks(), refinement); Element bottom_element = domain::Initialization::create_initial_element( - bottom_element_id, domain.blocks()[1], refinement); + bottom_element_id, domain.blocks(), refinement); // Subdomain LinearSolver::Schwarz::OverlapMap> overlap_elements{ {{Direction::upper_xi(), right_element_id}, @@ -270,9 +270,9 @@ auto make_databox_without_boundary_conditions() { const ElementId central_element_id{0, {{{2, 1}, {2, 1}}}}; const ElementId right_element_id{0, {{{2, 2}, {2, 1}}}}; Element central_element = domain::Initialization::create_initial_element( - central_element_id, domain.blocks()[0], refinement); + central_element_id, domain.blocks(), refinement); Element right_element = domain::Initialization::create_initial_element( - right_element_id, domain.blocks()[0], refinement); + right_element_id, domain.blocks(), refinement); LinearSolver::Schwarz::OverlapMap> overlap_elements{ {{Direction::upper_xi(), right_element_id}, std::move(right_element)}}; diff --git a/tests/Unit/Elliptic/Systems/Xcts/Events/Test_ObserveAdmIntegrals.cpp b/tests/Unit/Elliptic/Systems/Xcts/Events/Test_ObserveAdmIntegrals.cpp index 14a0080df979..8ee1fb0e53d4 100644 --- a/tests/Unit/Elliptic/Systems/Xcts/Events/Test_ObserveAdmIntegrals.cpp +++ b/tests/Unit/Elliptic/Systems/Xcts/Events/Test_ObserveAdmIntegrals.cpp @@ -105,7 +105,7 @@ void test_local_adm_integrals(const double& distance, // Get element information. const auto& current_block = blocks.at(element_id.block_id()); const auto current_element = domain::Initialization::create_initial_element( - element_id, current_block, initial_ref_levels); + element_id, blocks, initial_ref_levels); const ElementMap<3, Frame::Inertial> logical_to_inertial_map( element_id, current_block.stationary_map().get_clone()); diff --git a/tests/Unit/Evolution/DgSubcell/Test_BackgroundGrVars.cpp b/tests/Unit/Evolution/DgSubcell/Test_BackgroundGrVars.cpp index 85c1bc98f669..2eec7e9dcd67 100644 --- a/tests/Unit/Evolution/DgSubcell/Test_BackgroundGrVars.cpp +++ b/tests/Unit/Evolution/DgSubcell/Test_BackgroundGrVars.cpp @@ -166,7 +166,7 @@ void test(const gsl::not_null gen, const bool did_rollback) { const auto domain = brick.create_domain(); const auto element_id = ElementId<3>{0}; Element<3> element = domain::Initialization::create_initial_element( - element_id, domain.blocks().at(0), + element_id, domain.blocks(), std::vector>{{0, 0, 0}}); const Mesh<3> dg_mesh{num_dg_pts, Spectral::Basis::Legendre, diff --git a/tests/Unit/Evolution/DiscontinuousGalerkin/Test_BackgroundGrVars.cpp b/tests/Unit/Evolution/DiscontinuousGalerkin/Test_BackgroundGrVars.cpp index afb3193a00cc..eb79092742ad 100644 --- a/tests/Unit/Evolution/DiscontinuousGalerkin/Test_BackgroundGrVars.cpp +++ b/tests/Unit/Evolution/DiscontinuousGalerkin/Test_BackgroundGrVars.cpp @@ -131,7 +131,7 @@ void test(const gsl::not_null gen) { const auto domain = brick.create_domain(); const auto element_id = ElementId<3>{0}; Element<3> element = domain::Initialization::create_initial_element( - element_id, domain.blocks().at(0), + element_id, domain.blocks(), std::vector>{{0, 0, 0}}); const Mesh<3> mesh{num_dg_pts, Spectral::Basis::Legendre, diff --git a/tests/Unit/Evolution/Systems/Burgers/FiniteDifference/Test_BoundaryConditionGhostData.cpp b/tests/Unit/Evolution/Systems/Burgers/FiniteDifference/Test_BoundaryConditionGhostData.cpp index 3fff46dab117..07f053e6985c 100644 --- a/tests/Unit/Evolution/Systems/Burgers/FiniteDifference/Test_BoundaryConditionGhostData.cpp +++ b/tests/Unit/Evolution/Systems/Burgers/FiniteDifference/Test_BoundaryConditionGhostData.cpp @@ -94,7 +94,7 @@ void test(const BoundaryConditionType& boundary_condition) { auto domain = interval.create_domain(); auto boundary_conditions = interval.external_boundary_conditions(); const auto element = domain::Initialization::create_initial_element( - ElementId<1>{0, {SegmentId{0, 0}}}, domain.blocks().at(0), + ElementId<1>{0, {SegmentId{0, 0}}}, domain.blocks(), std::vector>{{refinement_level_x}}); // Mesh and coordinates diff --git a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/Test_Iterations.cpp b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/Test_Iterations.cpp index f8fd8b8c2fc4..19dc1f098c8d 100644 --- a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/Test_Iterations.cpp +++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/Test_Iterations.cpp @@ -213,7 +213,7 @@ void test_iterations(const size_t max_iterations) { for (const auto& element_id : element_ids) { const auto& my_block = blocks.at(element_id.block_id()); auto element = domain::Initialization::create_initial_element( - element_id, my_block, initial_refinements); + element_id, blocks, initial_refinements); auto mesh = domain::Initialization::create_initial_mesh( initial_extents, element_id, quadrature); const ElementMap element_map(element_id, diff --git a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/Test_ReceiveWorldtubeData.cpp b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/Test_ReceiveWorldtubeData.cpp index 0597c325fc7a..b569bdaf474e 100644 --- a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/Test_ReceiveWorldtubeData.cpp +++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/Test_ReceiveWorldtubeData.cpp @@ -174,7 +174,7 @@ SPECTRE_TEST_CASE("Unit.CurvedScalarWave.Worldtube.ReceiveWorldtubeData", const size_t grid_size = mesh.number_of_grid_points(); const auto& my_block = blocks.at(element_id.block_id()); auto element = domain::Initialization::create_initial_element( - element_id, my_block, initial_refinements); + element_id, blocks, initial_refinements); // we set lapse and shift to Minkowski so dt Psi = - Pi Scalar lapse(grid_size, 1.); diff --git a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/Test_SendToWorldtube.cpp b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/Test_SendToWorldtube.cpp index 01ec01ab6840..8617c423a0b5 100644 --- a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/Test_SendToWorldtube.cpp +++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/Test_SendToWorldtube.cpp @@ -216,7 +216,7 @@ SPECTRE_TEST_CASE("Unit.CurvedScalarWave.Worldtube.SendToWorldtube", "[Unit]") { for (const auto& element_id : element_ids) { const auto& my_block = blocks.at(element_id.block_id()); auto element = domain::Initialization::create_initial_element( - element_id, my_block, initial_refinements); + element_id, blocks, initial_refinements); auto mesh = domain::Initialization::create_initial_mesh( initial_extents, element_id, quadrature); const ElementMap element_map(element_id, diff --git a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/Test_Tags.cpp b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/Test_Tags.cpp index 09562b82a3ee..c169cac7fea6 100644 --- a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/Test_Tags.cpp +++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/Test_Tags.cpp @@ -137,7 +137,7 @@ void test_compute_face_coordinates_grid() { for (const auto& element_id : element_ids) { const auto& my_block = blocks.at(element_id.block_id()); const auto element = domain::Initialization::create_initial_element( - element_id, my_block, initial_refinements); + element_id, blocks, initial_refinements); const ElementMap element_map( element_id, my_block.stationary_map().get_to_grid_frame()); const auto mesh_1 = domain::Initialization::create_initial_mesh( @@ -210,7 +210,7 @@ void test_compute_face_coordinates() { for (const auto& element_id : element_ids) { const auto& my_block = blocks.at(element_id.block_id()); const auto element = domain::Initialization::create_initial_element( - element_id, my_block, initial_refinements); + element_id, blocks, initial_refinements); const auto mesh = domain::Initialization::create_initial_mesh( initial_extents, element_id, quadrature); const ElementMap element_map( @@ -544,7 +544,7 @@ void test_face_quantities_compute() { for (const auto& element_id : element_ids) { const auto& my_block = blocks.at(element_id.block_id()); const auto element = domain::Initialization::create_initial_element( - element_id, my_block, initial_refinements); + element_id, blocks, initial_refinements); const auto mesh = domain::Initialization::create_initial_mesh( initial_extents, element_id, quadrature); const ElementMap element_map(element_id, diff --git a/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Test_BoundaryConditionGhostData.cpp b/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Test_BoundaryConditionGhostData.cpp index 7eb1896a018f..129f52055e5c 100644 --- a/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Test_BoundaryConditionGhostData.cpp +++ b/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Test_BoundaryConditionGhostData.cpp @@ -105,8 +105,7 @@ void test(const BoundaryConditionType& boundary_condition, auto boundary_conditions = brick.external_boundary_conditions(); const auto element = domain::Initialization::create_initial_element( ElementId<3>{0, {SegmentId{0, 0}, SegmentId{0, 0}, SegmentId{0, 0}}}, - domain.blocks().at(0), - std::vector>{{refinement_levels}}); + domain.blocks(), std::vector>{{refinement_levels}}); // Mesh and coordinates const Mesh<3> dg_mesh{num_dg_pts, Spectral::Basis::Legendre, diff --git a/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/Test_NeighborPackagedData.cpp b/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/Test_NeighborPackagedData.cpp index 99d623519db5..a1c3e264e705 100644 --- a/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/Test_NeighborPackagedData.cpp +++ b/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/Test_NeighborPackagedData.cpp @@ -94,17 +94,19 @@ double test(const size_t num_dg_pts) { std::unordered_map> functions_of_time{}; - Block<3> block{ - domain::make_coordinate_map_base( - Affine3D{affine_map, affine_map, affine_map}), - 0, - {}}; + std::vector> blocks{}; + blocks.push_back( + {domain::make_coordinate_map_base( + Affine3D{affine_map, affine_map, affine_map}), + 0, + {}}); + const auto& block = blocks[0]; ElementMap<3, Frame::Grid> element_map{ element_id, block.is_time_dependent() ? block.moving_mesh_logical_to_grid_map().get_clone() : block.stationary_map().get_to_grid_frame()}; const auto element = domain::Initialization::create_initial_element( - element_id, block, + element_id, blocks, std::vector>{std::array{{3, 3, 3}}}); const auto moving_mesh_map = diff --git a/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/Test_TimeDerivative.cpp b/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/Test_TimeDerivative.cpp index 101cad26f9b4..acf4166f0831 100644 --- a/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/Test_TimeDerivative.cpp +++ b/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/Test_TimeDerivative.cpp @@ -131,16 +131,17 @@ double test(const size_t num_dg_pts, std::optional expansion_velocity, ::RelativisticEuler::Solutions::TovStar>>(soln)) .get_clone(); } - Block<3> block{test_non_diagonal_jacobian - ? domain::make_coordinate_map_base( - ::domain::CoordinateMaps::Rotation<3>(0.7, 0, 0.)) - : domain::make_coordinate_map_base( - Affine3D{affine_map, affine_map, affine_map}), - 0, - {}}; - + std::vector> blocks{}; + blocks.push_back({test_non_diagonal_jacobian + ? domain::make_coordinate_map_base( + ::domain::CoordinateMaps::Rotation<3>(0.7, 0, 0.)) + : domain::make_coordinate_map_base( + Affine3D{affine_map, affine_map, affine_map}), + 0, + {}}); + const auto& block = blocks[0]; std::unordered_map> functions_of_time{}; @@ -153,7 +154,7 @@ double test(const size_t num_dg_pts, std::optional expansion_velocity, ::domain::make_coordinate_map_base( ::domain::CoordinateMaps::Identity<3>{}); const auto element = domain::Initialization::create_initial_element( - element_id, block, + element_id, blocks, std::vector>{std::array{{3, 3, 3}}}); const double time = 0.0; diff --git a/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Test_BoundaryConditionGhostData.cpp b/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Test_BoundaryConditionGhostData.cpp index 149c4914ab65..921cd3808f1a 100644 --- a/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Test_BoundaryConditionGhostData.cpp +++ b/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Test_BoundaryConditionGhostData.cpp @@ -105,8 +105,7 @@ void test(const BoundaryConditionType& boundary_condition, auto boundary_conditions = brick.external_boundary_conditions(); const auto element = domain::Initialization::create_initial_element( ElementId<3>{0, {SegmentId{0, 0}, SegmentId{0, 0}, SegmentId{0, 0}}}, - domain.blocks().at(0), - std::vector>{{refinement_levels}}); + domain.blocks(), std::vector>{{refinement_levels}}); // Mesh and coordinates const Mesh<3> dg_mesh{num_dg_pts, Spectral::Basis::Legendre, diff --git a/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/Test_NeighborPackagedData.cpp b/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/Test_NeighborPackagedData.cpp index 7519a9c9707a..e084633897cb 100644 --- a/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/Test_NeighborPackagedData.cpp +++ b/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/Test_NeighborPackagedData.cpp @@ -104,17 +104,19 @@ double test(const size_t num_dg_pts) { std::unordered_map> functions_of_time{}; - Block<3> block{ - domain::make_coordinate_map_base( - Affine3D{affine_map, affine_map, affine_map}), - 0, - {}}; + std::vector> blocks{}; + blocks.push_back( + {domain::make_coordinate_map_base( + Affine3D{affine_map, affine_map, affine_map}), + 0, + {}}); + const auto& block = blocks[0]; ElementMap<3, Frame::Grid> element_map{ element_id, block.is_time_dependent() ? block.moving_mesh_logical_to_grid_map().get_clone() : block.stationary_map().get_to_grid_frame()}; const auto element = domain::Initialization::create_initial_element( - element_id, block, + element_id, blocks, std::vector>{std::array{{3, 3, 3}}}); const auto moving_mesh_map = diff --git a/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/Test_TimeDerivative.cpp b/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/Test_TimeDerivative.cpp index 5af78cbb1742..1b4fb7b53ff8 100644 --- a/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/Test_TimeDerivative.cpp +++ b/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/Test_TimeDerivative.cpp @@ -141,11 +141,12 @@ std::array test(const size_t num_dg_pts, const Affine affine_map{-1.0, 1.0, 1.0, 15.0}; const auto element = domain::Initialization::create_initial_element( ElementId<3>{0, {SegmentId{2, 2}, SegmentId{2, 2}, SegmentId{2, 2}}}, - Block<3>{domain::make_coordinate_map_base( - Affine3D{affine_map, affine_map, affine_map}), - 0, - {}}, + std::vector>{ + {domain::make_coordinate_map_base( + Affine3D{affine_map, affine_map, affine_map}), + 0, + {}}}, std::vector>{std::array{{3, 3, 3}}}); const ElementMap<3, Frame::Grid> element_map{ element.id(), diff --git a/tests/Unit/Evolution/Systems/NewtonianEuler/Subcell/Test_NeighborPackagedData.cpp b/tests/Unit/Evolution/Systems/NewtonianEuler/Subcell/Test_NeighborPackagedData.cpp index 67056983b064..3f0ea9c65c94 100644 --- a/tests/Unit/Evolution/Systems/NewtonianEuler/Subcell/Test_NeighborPackagedData.cpp +++ b/tests/Unit/Evolution/Systems/NewtonianEuler/Subcell/Test_NeighborPackagedData.cpp @@ -104,10 +104,11 @@ auto make_element<1>() { Affine affine_map{-1.0, 1.0, 2.0, 3.0}; return domain::Initialization::create_initial_element( ElementId<1>{0, {SegmentId{3, 4}}}, - Block<1>{domain::make_coordinate_map_base(affine_map), - 0, - {}}, + std::vector>{ + {domain::make_coordinate_map_base(affine_map), + 0, + {}}}, std::vector>{std::array{{3}}}); } @@ -116,11 +117,12 @@ auto make_element<2>() { Affine affine_map{-1.0, 1.0, 2.0, 3.0}; return domain::Initialization::create_initial_element( ElementId<2>{0, {SegmentId{3, 4}, SegmentId{3, 4}}}, - Block<2>{domain::make_coordinate_map_base( - Affine2D{affine_map, affine_map}), - 0, - {}}, + std::vector>{ + {domain::make_coordinate_map_base( + Affine2D{affine_map, affine_map}), + 0, + {}}}, std::vector>{std::array{{3, 3}}}); } @@ -129,11 +131,12 @@ auto make_element<3>() { Affine affine_map{-1.0, 1.0, 2.0, 3.0}; return domain::Initialization::create_initial_element( ElementId<3>{0, {SegmentId{3, 4}, SegmentId{3, 4}, SegmentId{3, 4}}}, - Block<3>{domain::make_coordinate_map_base( - Affine3D{affine_map, affine_map, affine_map}), - 0, - {}}, + std::vector>{ + {domain::make_coordinate_map_base( + Affine3D{affine_map, affine_map, affine_map}), + 0, + {}}}, std::vector>{std::array{{3, 3, 3}}}); } diff --git a/tests/Unit/Evolution/Systems/NewtonianEuler/Subcell/Test_TimeDerivative.cpp b/tests/Unit/Evolution/Systems/NewtonianEuler/Subcell/Test_TimeDerivative.cpp index 165cc10f9745..0e2158fcdc0c 100644 --- a/tests/Unit/Evolution/Systems/NewtonianEuler/Subcell/Test_TimeDerivative.cpp +++ b/tests/Unit/Evolution/Systems/NewtonianEuler/Subcell/Test_TimeDerivative.cpp @@ -97,10 +97,11 @@ auto make_element<1>() { Affine affine_map{-1.0, 1.0, 2.0, 3.0}; return domain::Initialization::create_initial_element( ElementId<1>{0, {SegmentId{2, 2}}}, - Block<1>{domain::make_coordinate_map_base(affine_map), - 0, - {}}, + std::vector>{ + {domain::make_coordinate_map_base(affine_map), + 0, + {}}}, std::vector>{std::array{{3}}}); } @@ -109,11 +110,12 @@ auto make_element<2>() { Affine affine_map{-1.0, 1.0, 2.0, 3.0}; return domain::Initialization::create_initial_element( ElementId<2>{0, {SegmentId{2, 2}, SegmentId{2, 2}}}, - Block<2>{domain::make_coordinate_map_base( - Affine2D{affine_map, affine_map}), - 0, - {}}, + std::vector>{ + {domain::make_coordinate_map_base( + Affine2D{affine_map, affine_map}), + 0, + {}}}, std::vector>{std::array{{3, 3}}}); } @@ -122,11 +124,12 @@ auto make_element<3>() { Affine affine_map{-1.0, 1.0, 2.0, 3.0}; return domain::Initialization::create_initial_element( ElementId<3>{0, {SegmentId{2, 2}, SegmentId{2, 2}, SegmentId{2, 2}}}, - Block<3>{domain::make_coordinate_map_base( - Affine3D{affine_map, affine_map, affine_map}), - 0, - {}}, + std::vector>{ + {domain::make_coordinate_map_base( + Affine3D{affine_map, affine_map, affine_map}), + 0, + {}}}, std::vector>{std::array{{3, 3, 3}}}); } diff --git a/tests/Unit/Helpers/Domain/DomainTestHelpers.cpp b/tests/Unit/Helpers/Domain/DomainTestHelpers.cpp index d012f93ec2aa..5e8047b3fb63 100644 --- a/tests/Unit/Helpers/Domain/DomainTestHelpers.cpp +++ b/tests/Unit/Helpers/Domain/DomainTestHelpers.cpp @@ -580,8 +580,7 @@ void test_initial_domain(const Domain& domain, for (const auto& element_id : element_ids) { elements.emplace(element_id, domain::Initialization::create_initial_element( - element_id, blocks[element_id.block_id()], - initial_refinement_levels)); + element_id, blocks, initial_refinement_levels)); } domain::test_domain_connectivity(domain, elements); domain::test_refinement_levels_of_neighbors<1>(elements); diff --git a/tests/Unit/Helpers/ParallelAlgorithms/LinearSolver/DistributedLinearSolverAlgorithmTestHelpers.hpp b/tests/Unit/Helpers/ParallelAlgorithms/LinearSolver/DistributedLinearSolverAlgorithmTestHelpers.hpp index 89032de8d438..81ceca2b93ab 100644 --- a/tests/Unit/Helpers/ParallelAlgorithms/LinearSolver/DistributedLinearSolverAlgorithmTestHelpers.hpp +++ b/tests/Unit/Helpers/ParallelAlgorithms/LinearSolver/DistributedLinearSolverAlgorithmTestHelpers.hpp @@ -288,9 +288,8 @@ struct InitializeElement { auto mesh = domain::Initialization::create_initial_mesh( initial_extents, element_id, Parallel::get(cache)); - const auto& block = domain.blocks()[element_id.block_id()]; auto element = domain::Initialization::create_initial_element( - element_id, block, initial_refinement); + element_id, domain.blocks(), initial_refinement); auto logical_coords = logical_coordinates(mesh); // Element data const size_t element_index = get_index(element_id); diff --git a/tests/Unit/PointwiseFunctions/Xcts/Test_AdmLinearMomentum.cpp b/tests/Unit/PointwiseFunctions/Xcts/Test_AdmLinearMomentum.cpp index 0f455943c6cf..f5d87f908410 100644 --- a/tests/Unit/PointwiseFunctions/Xcts/Test_AdmLinearMomentum.cpp +++ b/tests/Unit/PointwiseFunctions/Xcts/Test_AdmLinearMomentum.cpp @@ -72,7 +72,7 @@ void test_infinite_surface_integral(const double distance, const double mass, // Get element information const auto& current_block = blocks.at(element_id.block_id()); const auto current_element = domain::Initialization::create_initial_element( - element_id, current_block, initial_ref_levels); + element_id, blocks, initial_ref_levels); const ElementMap<3, Frame::Inertial> logical_to_inertial_map( element_id, current_block.stationary_map().get_clone()); @@ -186,7 +186,7 @@ void test_infinite_volume_integral(const double distance, const double mass, // Get element information const auto& current_block = blocks.at(element_id.block_id()); const auto current_element = domain::Initialization::create_initial_element( - element_id, current_block, initial_ref_levels); + element_id, blocks, initial_ref_levels); const ElementMap<3, Frame::Inertial> logical_to_inertial_map( element_id, current_block.stationary_map().get_clone()); diff --git a/tests/Unit/PointwiseFunctions/Xcts/Test_AdmMass.cpp b/tests/Unit/PointwiseFunctions/Xcts/Test_AdmMass.cpp index a7655108991f..6eac573d3b36 100644 --- a/tests/Unit/PointwiseFunctions/Xcts/Test_AdmMass.cpp +++ b/tests/Unit/PointwiseFunctions/Xcts/Test_AdmMass.cpp @@ -69,7 +69,7 @@ void test_infinite_surface_integral(const double distance, const double mass, // Get element information const auto& current_block = blocks.at(element_id.block_id()); const auto current_element = domain::Initialization::create_initial_element( - element_id, current_block, initial_ref_levels); + element_id, blocks, initial_ref_levels); const ElementMap<3, Frame::Inertial> logical_to_inertial_map( element_id, current_block.stationary_map().get_clone()); @@ -189,7 +189,7 @@ void test_infinite_volume_integral(const double distance, const double mass, // Get element information const auto& current_block = blocks.at(element_id.block_id()); const auto current_element = domain::Initialization::create_initial_element( - element_id, current_block, initial_ref_levels); + element_id, blocks, initial_ref_levels); const ElementMap<3, Frame::Inertial> logical_to_inertial_map( element_id, current_block.stationary_map().get_clone()); diff --git a/tests/Unit/PointwiseFunctions/Xcts/Test_CenterOfMass.cpp b/tests/Unit/PointwiseFunctions/Xcts/Test_CenterOfMass.cpp index c00f3ae7b461..32c34e091164 100644 --- a/tests/Unit/PointwiseFunctions/Xcts/Test_CenterOfMass.cpp +++ b/tests/Unit/PointwiseFunctions/Xcts/Test_CenterOfMass.cpp @@ -81,7 +81,7 @@ void test_infinite_surface_integral(const double distance, // Get element information const auto& current_block = blocks.at(element_id.block_id()); const auto current_element = domain::Initialization::create_initial_element( - element_id, current_block, initial_ref_levels); + element_id, blocks, initial_ref_levels); const ElementMap<3, Frame::Inertial> logical_to_inertial_map( element_id, current_block.stationary_map().get_clone()); @@ -169,7 +169,7 @@ void test_infinite_volume_integral(const double distance, // Get element information const auto& current_block = blocks.at(element_id.block_id()); const auto current_element = domain::Initialization::create_initial_element( - element_id, current_block, initial_ref_levels); + element_id, blocks, initial_ref_levels); const ElementMap<3, Frame::Inertial> logical_to_inertial_map( element_id, current_block.stationary_map().get_clone());