diff --git a/src/Domain/Creators/BinaryCompactObject.cpp b/src/Domain/Creators/BinaryCompactObject.cpp index e7d7bc265e92..52ee6dd0216f 100644 --- a/src/Domain/Creators/BinaryCompactObject.cpp +++ b/src/Domain/Creators/BinaryCompactObject.cpp @@ -269,8 +269,8 @@ BinaryCompactObject::BinaryCompactObject( << number_of_blocks_ << ")."); // Expand initial refinement and number of grid points over all blocks - const ExpandOverBlocks expand_over_blocks{block_names_, - block_groups_}; + const ExpandOverBlocks> expand_over_blocks{ + block_names_, block_groups_}; try { initial_refinement_ = std::visit(expand_over_blocks, initial_refinement); diff --git a/src/Domain/Creators/BlockGroups.cpp b/src/Domain/Creators/BlockGroups.cpp new file mode 100644 index 000000000000..4ac352d83ad7 --- /dev/null +++ b/src/Domain/Creators/BlockGroups.cpp @@ -0,0 +1,61 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Domain/Creators/BlockGroups.hpp" + +#include +#include +#include +#include +#include + +#include "Utilities/ErrorHandling/Error.hpp" +#include "Utilities/StdHelpers.hpp" + +namespace domain { + +bool block_is_in_group( + const std::string& block_name, const std::string& block_or_group_name, + const std::unordered_map>& + block_groups) { + if (block_name == block_or_group_name) { + return true; + } + const auto block_group_it = block_groups.find(block_or_group_name); + return block_group_it != block_groups.end() and + block_group_it->second.count(block_name) == 1; +} + +std::unordered_set expand_block_groups_to_block_names( + const std::vector& block_or_group_names, + const std::vector& all_block_names, + const std::unordered_map>& + block_groups) { + // Use an unordered_set to elide duplicates + std::unordered_set expanded_block_names; + for (const auto& block_or_group_name : block_or_group_names) { + if (const auto block_group_it = block_groups.find(block_or_group_name); + block_group_it != block_groups.end()) { + expanded_block_names.insert(block_group_it->second.begin(), + block_group_it->second.end()); + } else if (const auto block_name_it = + std::find(all_block_names.begin(), all_block_names.end(), + block_or_group_name); + block_name_it != all_block_names.end()) { + expanded_block_names.insert(*block_name_it); + } else { + using ::operator<<; + ERROR_AS( + "The block or group '" + << block_or_group_name + << "' is not one of the block names or groups of the domain. " + "The known block groups are " + << keys_of(block_groups) << " and the known block names are " + << all_block_names, + std::invalid_argument); + } + } + return expanded_block_names; +} + +} // namespace domain diff --git a/src/Domain/Creators/BlockGroups.hpp b/src/Domain/Creators/BlockGroups.hpp new file mode 100644 index 000000000000..0778a74ea8cf --- /dev/null +++ b/src/Domain/Creators/BlockGroups.hpp @@ -0,0 +1,45 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +/// \file +/// Utilities to work with block names and groups + +#pragma once + +#include +#include +#include +#include + +namespace domain { + +/*! + * \brief Check if a block is in the given group + * + * \param block_name The name of the block to check + * \param block_or_group_name The group name we're testing. Returns true if + * the block is in this group. This can also be the name of the block itself. + * \param block_groups All block groups. + */ +bool block_is_in_group( + const std::string& block_name, const std::string& block_or_group_name, + const std::unordered_map>& + block_groups); + +/*! + * \brief Expand a list of block or group names into a list of block names + * + * \param block_or_group_names Block or group names to expand + * \param all_block_names All block names in the domain + * \param block_groups Block groups used to expand the names + * \return std::unordered_set List of block names that appear in + * \p all_block_names. If one of the input names was a group, then all block + * names from that group are included. Overlaps between groups are allowed. + */ +std::unordered_set expand_block_groups_to_block_names( + const std::vector& block_or_group_names, + const std::vector& all_block_names, + const std::unordered_map>& + block_groups); + +} // namespace domain diff --git a/src/Domain/Creators/Brick.cpp b/src/Domain/Creators/Brick.cpp index 8d4d70da87e5..98d6de90b0d0 100644 --- a/src/Domain/Creators/Brick.cpp +++ b/src/Domain/Creators/Brick.cpp @@ -30,20 +30,18 @@ struct Inertial; namespace domain::creators { Brick::Brick( - typename LowerBound::type lower_xyz, typename UpperBound::type upper_xyz, - typename InitialRefinement::type initial_refinement_level_xyz, - typename InitialGridPoints::type initial_number_of_grid_points_in_xyz, - typename IsPeriodicIn::type is_periodic_in_xyz, + std::array lower_xyz, std::array upper_xyz, + std::array initial_refinement_level_xyz, + std::array initial_number_of_grid_points_in_xyz, + std::array is_periodic_in_xyz, std::unique_ptr> time_dependence) - // clang-tidy: trivially copyable - : lower_xyz_(std::move(lower_xyz)), // NOLINT - upper_xyz_(std::move(upper_xyz)), // NOLINT - is_periodic_in_xyz_(std::move(is_periodic_in_xyz)), // NOLINT - initial_refinement_level_xyz_( // NOLINT - std::move(initial_refinement_level_xyz)), // NOLINT - initial_number_of_grid_points_in_xyz_( // NOLINT - std::move(initial_number_of_grid_points_in_xyz)), // NOLINT + : lower_xyz_(lower_xyz), + upper_xyz_(upper_xyz), + is_periodic_in_xyz_(is_periodic_in_xyz), + initial_refinement_level_xyz_(initial_refinement_level_xyz), + initial_number_of_grid_points_in_xyz_( + initial_number_of_grid_points_in_xyz), time_dependence_(std::move(time_dependence)), boundary_condition_in_lower_x_(nullptr), boundary_condition_in_upper_x_(nullptr), @@ -58,9 +56,9 @@ Brick::Brick( } Brick::Brick( - typename LowerBound::type lower_xyz, typename UpperBound::type upper_xyz, - typename InitialRefinement::type initial_refinement_level_xyz, - typename InitialGridPoints::type initial_number_of_grid_points_in_xyz, + std::array lower_xyz, std::array upper_xyz, + std::array initial_refinement_level_xyz, + std::array initial_number_of_grid_points_in_xyz, std::unique_ptr boundary_condition_in_lower_x, std::unique_ptr @@ -76,14 +74,12 @@ Brick::Brick( std::unique_ptr> time_dependence, const Options::Context& context) - // clang-tidy: trivially copyable - : lower_xyz_(std::move(lower_xyz)), // NOLINT - upper_xyz_(std::move(upper_xyz)), // NOLINT + : lower_xyz_(lower_xyz), + upper_xyz_(upper_xyz), is_periodic_in_xyz_{{false, false, false}}, - initial_refinement_level_xyz_( // NOLINT - std::move(initial_refinement_level_xyz)), // NOLINT - initial_number_of_grid_points_in_xyz_( // NOLINT - std::move(initial_number_of_grid_points_in_xyz)), // NOLINT + initial_refinement_level_xyz_(initial_refinement_level_xyz), + initial_number_of_grid_points_in_xyz_( + initial_number_of_grid_points_in_xyz), time_dependence_(std::move(time_dependence)), boundary_condition_in_lower_x_(std::move(boundary_condition_in_lower_x)), boundary_condition_in_upper_x_(std::move(boundary_condition_in_upper_x)), diff --git a/src/Domain/Creators/Brick.hpp b/src/Domain/Creators/Brick.hpp index 7852acecf23b..2846def0973d 100644 --- a/src/Domain/Creators/Brick.hpp +++ b/src/Domain/Creators/Brick.hpp @@ -152,18 +152,16 @@ class Brick : public DomainCreator<3> { static constexpr Options::String help{"Creates a 3D brick."}; - Brick(typename LowerBound::type lower_xyz, - typename UpperBound::type upper_xyz, - typename InitialRefinement::type initial_refinement_level_xyz, - typename InitialGridPoints::type initial_number_of_grid_points_in_xyz, - typename IsPeriodicIn::type is_periodic_in_xyz, + Brick(std::array lower_xyz, std::array upper_xyz, + std::array initial_refinement_level_xyz, + std::array initial_number_of_grid_points_in_xyz, + std::array is_periodic_in_xyz, std::unique_ptr> time_dependence = nullptr); - Brick(typename LowerBound::type lower_xyz, - typename UpperBound::type upper_xyz, - typename InitialRefinement::type initial_refinement_level_xyz, - typename InitialGridPoints::type initial_number_of_grid_points_in_xyz, + Brick(std::array lower_xyz, std::array upper_xyz, + std::array initial_refinement_level_xyz, + std::array initial_number_of_grid_points_in_xyz, std::unique_ptr boundary_condition_in_lower_x, std::unique_ptr @@ -181,10 +179,9 @@ class Brick : public DomainCreator<3> { const Options::Context& context = {}); template - Brick(typename LowerBound::type lower_xyz, - typename UpperBound::type upper_xyz, - typename InitialRefinement::type initial_refinement_level_xyz, - typename InitialGridPoints::type initial_number_of_grid_points_in_xyz, + Brick(std::array lower_xyz, std::array upper_xyz, + std::array initial_refinement_level_xyz, + std::array initial_number_of_grid_points_in_xyz, std::variant, LowerUpperBoundaryCondition> boundary_conditions_in_x, @@ -248,9 +245,9 @@ class Brick : public DomainCreator<3> { template Brick::Brick( - typename LowerBound::type lower_xyz, typename UpperBound::type upper_xyz, - typename InitialRefinement::type initial_refinement_level_xyz, - typename InitialGridPoints::type initial_number_of_grid_points_in_xyz, + std::array lower_xyz, std::array upper_xyz, + std::array initial_refinement_level_xyz, + std::array initial_number_of_grid_points_in_xyz, std::variant, LowerUpperBoundaryCondition> boundary_conditions_in_x, diff --git a/src/Domain/Creators/CMakeLists.txt b/src/Domain/Creators/CMakeLists.txt index f46210f21fb3..a6a2faeed39b 100644 --- a/src/Domain/Creators/CMakeLists.txt +++ b/src/Domain/Creators/CMakeLists.txt @@ -11,6 +11,7 @@ spectre_target_sources( AlignedLattice.cpp BinaryCompactObject.cpp BinaryCompactObjectHelpers.cpp + BlockGroups.cpp Brick.cpp Cylinder.cpp CylindricalBinaryCompactObject.cpp @@ -34,12 +35,14 @@ spectre_target_headers( AlignedLattice.hpp BinaryCompactObject.hpp BinaryCompactObjectHelpers.hpp + BlockGroups.hpp Brick.hpp Cylinder.hpp CylindricalBinaryCompactObject.hpp Disk.hpp DomainCreator.hpp ExpandOverBlocks.hpp + ExpandOverBlocks.tpp Factory.hpp Factory1D.hpp Factory2D.hpp diff --git a/src/Domain/Creators/Cylinder.cpp b/src/Domain/Creators/Cylinder.cpp index 59745aac3c6d..c44b36e5d3c3 100644 --- a/src/Domain/Creators/Cylinder.cpp +++ b/src/Domain/Creators/Cylinder.cpp @@ -174,8 +174,8 @@ Cylinder::Cylinder( } // Expand initial refinement and number of grid points over all blocks - const ExpandOverBlocks expand_over_blocks{block_names_, - block_groups_}; + const ExpandOverBlocks> expand_over_blocks{ + block_names_, block_groups_}; try { initial_refinement_ = std::visit(expand_over_blocks, initial_refinement); } catch (const std::exception& error) { diff --git a/src/Domain/Creators/CylindricalBinaryCompactObject.cpp b/src/Domain/Creators/CylindricalBinaryCompactObject.cpp index 12e9b058f4b4..272e650802f0 100644 --- a/src/Domain/Creators/CylindricalBinaryCompactObject.cpp +++ b/src/Domain/Creators/CylindricalBinaryCompactObject.cpp @@ -304,8 +304,8 @@ CylindricalBinaryCompactObject::CylindricalBinaryCompactObject( } // Expand initial refinement over all blocks - const ExpandOverBlocks expand_over_blocks{block_names_, - block_groups_}; + const ExpandOverBlocks> expand_over_blocks{ + block_names_, block_groups_}; try { initial_refinement_ = std::visit(expand_over_blocks, initial_refinement); } catch (const std::exception& error) { diff --git a/src/Domain/Creators/ExpandOverBlocks.cpp b/src/Domain/Creators/ExpandOverBlocks.cpp index b27644b385d1..37bd10396f83 100644 --- a/src/Domain/Creators/ExpandOverBlocks.cpp +++ b/src/Domain/Creators/ExpandOverBlocks.cpp @@ -1,109 +1,15 @@ // Distributed under the MIT License. // See LICENSE.txt for details. -#include "Domain/Creators/ExpandOverBlocks.hpp" +#include "Domain/Creators/ExpandOverBlocks.tpp" #include -#include #include -#include -#include -#include - -#include "Utilities/ErrorHandling/Assert.hpp" -#include "Utilities/GenerateInstantiations.hpp" -#include "Utilities/MakeArray.hpp" namespace domain { -template -ExpandOverBlocks::ExpandOverBlocks(size_t num_blocks) - : num_blocks_(num_blocks) {} - -template -ExpandOverBlocks::ExpandOverBlocks( - std::vector block_names, - std::unordered_map> - block_groups) - : num_blocks_(block_names.size()), - block_names_(std::move(block_names)), - block_groups_(std::move(block_groups)) {} - -template -std::vector> ExpandOverBlocks::operator()( - const T& value) const { - return {num_blocks_, make_array(value)}; -} - -template -std::vector> ExpandOverBlocks::operator()( - const std::array& value) const { - return {num_blocks_, value}; -} - -template -std::vector> ExpandOverBlocks::operator()( - const std::vector> value) const { - if (value.size() != num_blocks_) { - throw std::length_error{"You supplied " + std::to_string(value.size()) + - " values, but the domain creator has " + - std::to_string(num_blocks_) + " blocks."}; - } - return value; -} - -template -std::vector> ExpandOverBlocks::operator()( - const std::unordered_map>& value) const { - ASSERT(num_blocks_ == block_names_.size(), - "Construct 'ExpandOverBlocks' with block names to use the " - "map-over-block-names feature."); - // Expand group names - auto value_per_block = value; - for (const auto& [name, block_value] : value) { - const auto found_group = block_groups_.find(name); - if (found_group != block_groups_.end()) { - for (const auto& expanded_name : found_group->second) { - if (value_per_block.count(expanded_name) == 0) { - value_per_block[expanded_name] = block_value; - } else { - throw std::invalid_argument{ - "Duplicate block name '" + expanded_name + - // NOLINTNEXTLINE(performance-inefficient-string-concatenation) - "' (expanded from '" + name + "')."}; - } - } - value_per_block.erase(name); - } - } - if (value_per_block.size() != num_blocks_) { - throw std::length_error{ - "You supplied " + std::to_string(value_per_block.size()) + - " values, but the domain creator has " + std::to_string(num_blocks_) + - " blocks: " + boost::algorithm::join(block_names_, ", ")}; - } - std::vector> result{}; - for (const auto& block_name : block_names_) { - const auto found_value = value_per_block.find(block_name); - if (found_value != value_per_block.end()) { - result.push_back(found_value->second); - } else { - throw std::out_of_range{"Value for block '" + block_name + - "' is missing. Did you misspell its name?"}; - } - } - return result; -} - -#define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data) -#define DIM(data) BOOST_PP_TUPLE_ELEM(1, data) -#define INSTANTIATE(_, data) \ - template class ExpandOverBlocks; - -GENERATE_INSTANTIATIONS(INSTANTIATE, (size_t), (1, 2, 3)) - -#undef DIM -#undef DTYPE -#undef INSTANTIATE +template class ExpandOverBlocks>; +template class ExpandOverBlocks>; +template class ExpandOverBlocks>; } // namespace domain diff --git a/src/Domain/Creators/ExpandOverBlocks.hpp b/src/Domain/Creators/ExpandOverBlocks.hpp index f886f763ecb2..8afe20643130 100644 --- a/src/Domain/Creators/ExpandOverBlocks.hpp +++ b/src/Domain/Creators/ExpandOverBlocks.hpp @@ -10,23 +10,44 @@ #include #include +#include "Utilities/MakeArray.hpp" + namespace domain { + +namespace ExpandOverBlocks_detail { +template > +struct is_value_type { + static constexpr bool value = false; +}; + +template +struct is_value_type> { + static constexpr bool value = std::is_same_v; +}; +} // namespace ExpandOverBlocks_detail + /*! - * \brief Produce a distribution of type `T` over all blocks and dimensions in - * the domain, based on values `T` of variable isotropy and homogeneity. + * \brief Produce a std::vector over all blocks of the domain * * This class is useful to option-create values for e.g. the initial refinement * level or initial number of grid points for domain creators. It can be used * with `std::visit` and a `std::variant` with (a subset of) these types: * - * - `T`: Repeat over all blocks and dimensions (isotropic and homogeneous). - * - `std::array`: Repeat over all blocks (homogeneous). - * - `std::vector>>`: Only check if the size matches the + * - `T`: Repeat the given value over all blocks (homogeneous). + * - `std::vector`: Only check if the size matches the * number of blocks, throwing a `std::length_error` if it doesn't. - * - `std::unordered_map>`: Map block names, or + * - `std::unordered_map`: Map block names, or * names of block groups, to values. The map must cover all blocks once the * groups are expanded. To use this option you must pass the list of block * names and groups to the constructor. + * - `T::value_type`: Repeat the given value over all blocks and dimensions + * (isotropic and homogeneous). Only works if `T` is a `std::array`. For + * example, if `T` is `std::array`, this will produce a + * `std::vector>` with the same value repeated + * `num_blocks x 3` times. + * + * If `T` is a `std::unique_ptr`, the class will clone the value for each block + * using `T::get_clone()`. * * Note that the call-operators `throw` when they encounter errors, such as * mismatches in the number of blocks. The exceptions can be used to output @@ -41,9 +62,8 @@ namespace domain { * \snippet Test_ExpandOverBlocks.cpp expand_over_blocks_named_example * * \tparam T The type distributed over the domain - * \tparam Dim The number of spatial dimensions */ -template +template struct ExpandOverBlocks { ExpandOverBlocks(size_t num_blocks); ExpandOverBlocks( @@ -51,25 +71,28 @@ struct ExpandOverBlocks { std::unordered_map> block_groups = {}); - /// Repeat over all blocks and dimensions (isotropic and homogeneous) - std::vector> operator()(const T& value) const; - /// Repeat over all blocks (homogeneous) - std::vector> operator()( - const std::array& value) const; + std::vector operator()(const T& value) const; /// Only check if the size matches the number of blocks, throwing a /// `std::length_error` if it doesn't - std::vector> operator()( - std::vector> value) const; + std::vector operator()(const std::vector& value) const; /// Map block names, or names of block groups, to values. The map must cover /// all blocks once the groups are expanded. To use this option you must pass /// the list of block names and groups to the constructor. Here's an example: /// /// \snippet Test_ExpandOverBlocks.cpp expand_over_blocks_named_example - std::vector> operator()( - const std::unordered_map>& value) const; + std::vector operator()( + const std::unordered_map& value) const; + + /// Repeat over all blocks and dimensions (isotropic and homogeneous) + template < + typename U, + Requires::value> = nullptr> + std::vector operator()(const U& value) const { + return {num_blocks_, make_array>(value)}; + } private: size_t num_blocks_; @@ -77,4 +100,5 @@ struct ExpandOverBlocks { std::unordered_map> block_groups_; }; + } // namespace domain diff --git a/src/Domain/Creators/ExpandOverBlocks.tpp b/src/Domain/Creators/ExpandOverBlocks.tpp new file mode 100644 index 000000000000..32760a4ac8d0 --- /dev/null +++ b/src/Domain/Creators/ExpandOverBlocks.tpp @@ -0,0 +1,121 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include "Domain/Creators/ExpandOverBlocks.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "Utilities/CloneUniquePtrs.hpp" +#include "Utilities/ErrorHandling/Assert.hpp" +#include "Utilities/TypeTraits/IsA.hpp" + +namespace domain { + +template +ExpandOverBlocks::ExpandOverBlocks(size_t num_blocks) + : num_blocks_(num_blocks) {} + +template +ExpandOverBlocks::ExpandOverBlocks( + std::vector block_names, + std::unordered_map> + block_groups) + : num_blocks_(block_names.size()), + block_names_(std::move(block_names)), + block_groups_(std::move(block_groups)) {} + +template +std::vector ExpandOverBlocks::operator()(const T& value) const { + if constexpr (tt::is_a_v) { + std::vector expanded(num_blocks_); + for (size_t i = 0; i < num_blocks_; ++i) { + expanded[i] = value->get_clone(); + } + return expanded; + } else { + return {num_blocks_, value}; + } +} + +template +std::vector ExpandOverBlocks::operator()( + const std::vector& value) const { + if (value.size() != num_blocks_) { + throw std::length_error{"You supplied " + std::to_string(value.size()) + + " values, but the domain creator has " + + std::to_string(num_blocks_) + " blocks."}; + } + if constexpr (tt::is_a_v) { + return clone_unique_ptrs(value); + } else { + return value; + } +} + +template +std::vector ExpandOverBlocks::operator()( + const std::unordered_map& value) const { + ASSERT(num_blocks_ == block_names_.size(), + "Construct 'ExpandOverBlocks' with block names to use the " + "map-over-block-names feature."); + // Expand group names + auto value_per_block = [&value]() { + if constexpr (tt::is_a_v) { + return clone_unique_ptrs(value); + } else { + return value; + } + }(); + for (const auto& [name, block_value] : value) { + const auto found_group = block_groups_.find(name); + if (found_group != block_groups_.end()) { + for (const auto& expanded_name : found_group->second) { + if (value_per_block.count(expanded_name) == 0) { + value_per_block[expanded_name] = [&local_block_value = + block_value]() { + if constexpr (tt::is_a_v) { + return local_block_value->get_clone(); + } else { + return local_block_value; + } + }(); + } else { + throw std::invalid_argument{ + "Duplicate block name '" + expanded_name + + // NOLINTNEXTLINE(performance-inefficient-string-concatenation) + "' (expanded from '" + name + "')."}; + } + } + value_per_block.erase(name); + } + } + if (value_per_block.size() != num_blocks_) { + throw std::length_error{ + "You supplied " + std::to_string(value_per_block.size()) + + " values, but the domain creator has " + std::to_string(num_blocks_) + + " blocks: " + boost::algorithm::join(block_names_, ", ")}; + } + std::vector result{}; + result.reserve(num_blocks_); + for (const auto& block_name : block_names_) { + const auto found_value = value_per_block.find(block_name); + if (found_value != value_per_block.end()) { + result.emplace_back(std::move(found_value->second)); + } else { + throw std::out_of_range{"Value for block '" + block_name + + "' is missing. Did you misspell its name?"}; + } + } + return result; +} + +} // namespace domain diff --git a/src/Domain/Creators/Rectangle.cpp b/src/Domain/Creators/Rectangle.cpp index a14859d7d178..118ab7f3c249 100644 --- a/src/Domain/Creators/Rectangle.cpp +++ b/src/Domain/Creators/Rectangle.cpp @@ -34,20 +34,17 @@ struct BlockLogical; namespace domain::creators { Rectangle::Rectangle( - typename LowerBound::type lower_xy, typename UpperBound::type upper_xy, - typename InitialRefinement::type initial_refinement_level_xy, - typename InitialGridPoints::type initial_number_of_grid_points_in_xy, - typename IsPeriodicIn::type is_periodic_in_xy, + std::array lower_xy, std::array upper_xy, + std::array initial_refinement_level_xy, + std::array initial_number_of_grid_points_in_xy, + std::array is_periodic_in_xy, std::unique_ptr> time_dependence) - // clang-tidy: trivially copyable - : lower_xy_(std::move(lower_xy)), // NOLINT - upper_xy_(std::move(upper_xy)), // NOLINT - is_periodic_in_xy_(std::move(is_periodic_in_xy)), // NOLINT - initial_refinement_level_xy_( // NOLINT - std::move(initial_refinement_level_xy)), // NOLINT - initial_number_of_grid_points_in_xy_( // NOLINT - std::move(initial_number_of_grid_points_in_xy)), // NOLINT + : lower_xy_(lower_xy), + upper_xy_(upper_xy), + is_periodic_in_xy_(is_periodic_in_xy), + initial_refinement_level_xy_(initial_refinement_level_xy), + initial_number_of_grid_points_in_xy_(initial_number_of_grid_points_in_xy), time_dependence_(std::move(time_dependence)), boundary_condition_(nullptr) { if (time_dependence_ == nullptr) { @@ -57,9 +54,9 @@ Rectangle::Rectangle( } Rectangle::Rectangle( - typename LowerBound::type lower_xy, typename UpperBound::type upper_xy, - typename InitialRefinement::type initial_refinement_level_xy, - typename InitialGridPoints::type initial_number_of_grid_points_in_xy, + std::array lower_xy, std::array upper_xy, + std::array initial_refinement_level_xy, + std::array initial_number_of_grid_points_in_xy, std::unique_ptr boundary_condition, std::unique_ptr> @@ -105,7 +102,10 @@ Domain<2> Rectangle::create_domain() const { make_vector_coordinate_map_base( Affine2D{Affine{-1., 1., lower_xy_[0], upper_xy_[0]}, Affine{-1., 1., lower_xy_[1], upper_xy_[1]}}), - std::vector>{{{0, 1, 2, 3}}}, identifications}; + std::vector>{{{0, 1, 2, 3}}}, + identifications, + {}, + block_names_}; if (not time_dependence_->is_none()) { domain.inject_time_dependent_map_for_block( diff --git a/src/Domain/Creators/Rectangle.hpp b/src/Domain/Creators/Rectangle.hpp index 532a1f6f11a4..a65f495ade08 100644 --- a/src/Domain/Creators/Rectangle.hpp +++ b/src/Domain/Creators/Rectangle.hpp @@ -9,6 +9,9 @@ #include #include #include +#include +#include +#include #include #include "Domain/BoundaryConditions/BoundaryCondition.hpp" @@ -107,17 +110,17 @@ class Rectangle : public DomainCreator<2> { static constexpr Options::String help{"Creates a 2D rectangle."}; Rectangle( - typename LowerBound::type lower_xy, typename UpperBound::type upper_xy, - typename InitialRefinement::type initial_refinement_level_xy, - typename InitialGridPoints::type initial_number_of_grid_points_in_xy, - typename IsPeriodicIn::type is_periodic_in_xy, + std::array lower_xy, std::array upper_xy, + std::array initial_refinement_level_xy, + std::array initial_number_of_grid_points_in_xy, + std::array is_periodic_in_xy, std::unique_ptr> time_dependence = nullptr); Rectangle( - typename LowerBound::type lower_xy, typename UpperBound::type upper_xy, - typename InitialRefinement::type initial_refinement_level_xy, - typename InitialGridPoints::type initial_number_of_grid_points_in_xy, + std::array lower_xy, std::array upper_xy, + std::array initial_refinement_level_xy, + std::array initial_number_of_grid_points_in_xy, std::unique_ptr boundary_condition, std::unique_ptr> @@ -147,6 +150,8 @@ class Rectangle : public DomainCreator<2> { std::string, std::unique_ptr> override; + std::vector block_names() const override { return block_names_; } + private: typename LowerBound::type lower_xy_{}; typename UpperBound::type upper_xy_{}; @@ -157,6 +162,7 @@ class Rectangle : public DomainCreator<2> { time_dependence_{nullptr}; std::unique_ptr boundary_condition_{nullptr}; + inline static const std::vector block_names_{"Rectangle"}; }; } // namespace creators } // namespace domain diff --git a/src/Domain/Creators/Sphere.cpp b/src/Domain/Creators/Sphere.cpp index 64dc617554e8..b26eceb45086 100644 --- a/src/Domain/Creators/Sphere.cpp +++ b/src/Domain/Creators/Sphere.cpp @@ -187,8 +187,8 @@ Sphere::Sphere( << num_blocks_ << " but is " << block_names_.size() << "."); // Expand initial refinement and number of grid points over all blocks - const ExpandOverBlocks expand_over_blocks{block_names_, - block_groups_}; + const ExpandOverBlocks> expand_over_blocks{ + block_names_, block_groups_}; try { initial_refinement_ = std::visit(expand_over_blocks, initial_refinement); } catch (const std::exception& error) { diff --git a/src/Elliptic/Executables/Elasticity/CMakeLists.txt b/src/Elliptic/Executables/Elasticity/CMakeLists.txt index 0b3a1081c639..9e806a609db0 100644 --- a/src/Elliptic/Executables/Elasticity/CMakeLists.txt +++ b/src/Elliptic/Executables/Elasticity/CMakeLists.txt @@ -9,6 +9,7 @@ set(LIBS_TO_LINK DiscontinuousGalerkin DomainCreators Elasticity + ElasticityActions ElasticityBoundaryConditions ElasticityPointwiseFunctions ElasticitySolutions diff --git a/src/Elliptic/Executables/Elasticity/SolveElasticity.hpp b/src/Elliptic/Executables/Elasticity/SolveElasticity.hpp index 63312741bbf9..f9c88ef9dec3 100644 --- a/src/Elliptic/Executables/Elasticity/SolveElasticity.hpp +++ b/src/Elliptic/Executables/Elasticity/SolveElasticity.hpp @@ -14,6 +14,7 @@ #include "Elliptic/BoundaryConditions/BoundaryCondition.hpp" #include "Elliptic/DiscontinuousGalerkin/DgElementArray.hpp" #include "Elliptic/Executables/Solver.hpp" +#include "Elliptic/Systems/Elasticity/Actions/InitializeConstitutiveRelation.hpp" #include "Elliptic/Systems/Elasticity/BoundaryConditions/Factory.hpp" #include "Elliptic/Systems/Elasticity/FirstOrderSystem.hpp" #include "Elliptic/Systems/Elasticity/Tags.hpp" @@ -78,8 +79,10 @@ struct Metavariables { error_compute>; // Collect all items to store in the cache. - using const_global_cache_tags = - tmpl::list>; + using const_global_cache_tags = tmpl::list<>; + + // Just a name in the input file + struct ObserveNormsPerLayer {}; struct factory_creation : tt::ConformsTo { @@ -98,12 +101,17 @@ struct Metavariables { Elasticity::ConstitutiveRelations::ConstitutiveRelation, Elasticity::ConstitutiveRelations::standard_constitutive_relations< volume_dim>>, - tmpl::pair>>>, + tmpl::pair< + Event, + tmpl::flatten, + // Observation per material layer + ::Events::ObserveNorms>>>, tmpl::pair>, tmpl::pair, Parallel::Actions::TerminatePhase>; using register_actions = diff --git a/src/Elliptic/Systems/Elasticity/Actions/CMakeLists.txt b/src/Elliptic/Systems/Elasticity/Actions/CMakeLists.txt new file mode 100644 index 000000000000..1f1567c39e2c --- /dev/null +++ b/src/Elliptic/Systems/Elasticity/Actions/CMakeLists.txt @@ -0,0 +1,31 @@ +# Distributed under the MIT License. +# See LICENSE.txt for details. + +set(LIBRARY ElasticityActions) + +add_spectre_library(${LIBRARY}) + +spectre_target_sources( + ${LIBRARY} + PRIVATE + InitializeConstitutiveRelation.cpp + ) + +spectre_target_headers( + ${LIBRARY} + INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src + HEADERS + InitializeConstitutiveRelation.hpp + ) + +target_link_libraries( + ${LIBRARY} + PUBLIC + ConstitutiveRelations + DomainCreators + INTERFACE + DataStructures + Domain + Parallel + Utilities + ) diff --git a/src/Elliptic/Systems/Elasticity/Actions/InitializeConstitutiveRelation.cpp b/src/Elliptic/Systems/Elasticity/Actions/InitializeConstitutiveRelation.cpp new file mode 100644 index 000000000000..11dde7e15ec5 --- /dev/null +++ b/src/Elliptic/Systems/Elasticity/Actions/InitializeConstitutiveRelation.cpp @@ -0,0 +1,17 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +// Instantiations of domain::ExpandOverBlocks for Elasticity + +#include +#include + +#include "Domain/Creators/ExpandOverBlocks.tpp" +#include "PointwiseFunctions/Elasticity/ConstitutiveRelations/ConstitutiveRelation.hpp" + +template +using ConstRelPtr = std::unique_ptr< + Elasticity::ConstitutiveRelations::ConstitutiveRelation>; + +template class domain::ExpandOverBlocks>; +template class domain::ExpandOverBlocks>; diff --git a/src/Elliptic/Systems/Elasticity/Actions/InitializeConstitutiveRelation.hpp b/src/Elliptic/Systems/Elasticity/Actions/InitializeConstitutiveRelation.hpp new file mode 100644 index 000000000000..eac8851b25f6 --- /dev/null +++ b/src/Elliptic/Systems/Elasticity/Actions/InitializeConstitutiveRelation.hpp @@ -0,0 +1,201 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "DataStructures/DataBox/DataBox.hpp" +#include "DataStructures/DataBox/Tag.hpp" +#include "Domain/Creators/BlockGroups.hpp" +#include "Domain/Creators/DomainCreator.hpp" +#include "Domain/Creators/ExpandOverBlocks.hpp" +#include "Domain/Creators/OptionTags.hpp" +#include "Domain/Creators/Tags/Domain.hpp" +#include "Domain/Structure/ElementId.hpp" +#include "Domain/Tags.hpp" +#include "Elliptic/Tags.hpp" +#include "IO/Observer/Tags.hpp" +#include "Options/String.hpp" +#include "Parallel/AlgorithmExecution.hpp" +#include "PointwiseFunctions/Elasticity/ConstitutiveRelations/ConstitutiveRelation.hpp" +#include "PointwiseFunctions/Elasticity/ConstitutiveRelations/Tags.hpp" +#include "Utilities/CallWithDynamicType.hpp" +#include "Utilities/ErrorHandling/Error.hpp" +#include "Utilities/TMPL.hpp" + +/// \cond +namespace tuples { +template +struct TaggedTuple; +} // namespace tuples +namespace Parallel { +template +struct GlobalCache; +} // namespace Parallel +/// \endcond + +namespace Elasticity { +namespace Tags { + +/// A constitutive relation in every block of the domain +template +struct ConstitutiveRelationPerBlock : db::SimpleTag, + ConstitutiveRelationPerBlockBase { + using ConstRelPtr = + std::unique_ptr>; + using type = std::vector; + + using option_tags = tmpl::list, + OptionTags::ConstitutiveRelationPerBlock>; + static constexpr bool pass_metavariables = false; + + static type create_from_options( + const std::unique_ptr>& domain_creator, + const std::variant, + std::unordered_map>& + constitutive_relation_per_block) { + const auto block_names = domain_creator->block_names(); + const auto block_groups = domain_creator->block_groups(); + const domain::ExpandOverBlocks expand_over_blocks{ + block_names, block_groups}; + try { + return std::visit(expand_over_blocks, constitutive_relation_per_block); + } catch (const std::exception& error) { + ERROR_NO_TRACE("Invalid 'Material':\n" << error.what()); + } + } +}; + +/// References the constitutive relation for the element's block, which is +/// stored in the global cache +template +struct ConstitutiveRelationReference : ConstitutiveRelation, + db::ReferenceTag { + using base = ConstitutiveRelation; + using argument_tags = + tmpl::list>; + static const ConstitutiveRelations::ConstitutiveRelation& get( + const std::vector< + std::unique_ptr>>& + constitutive_relation_per_block, + const Element& element) { + return *constitutive_relation_per_block.at(element.id().block_id()); + } +}; + +/// Stores the names of the block groups that split the domain into layers with +/// different material properties. Useful to observe quantities in each layer. +template +struct MaterialBlockGroups : db::SimpleTag { + using type = std::unordered_set; + + using option_tags = tmpl::list>; + static constexpr bool pass_metavariables = false; + using ConstRelPtr = + std::unique_ptr>; + + static type create_from_options( + const std::variant, + std::unordered_map>& + constitutive_relation_per_block) { + if (std::holds_alternative>( + constitutive_relation_per_block)) { + const auto& map = std::get>( + constitutive_relation_per_block); + std::unordered_set block_groups; + for (const auto& [block_name, _] : map) { + block_groups.insert(block_name); + } + return block_groups; + } else { + return {}; + } + } +}; + +/// The name of the material layer (name of a block group with some material) +struct MaterialLayerName : db::SimpleTag { + using type = std::optional; +}; + +} // namespace Tags + +/// Actions related to solving Elasticity systems +namespace Actions { + +/*! + * \brief Initialize the constitutive relation describing properties of the + * elastic material + * + * Every block in the domain can have a different constitutive relation, + * allowing for composite materials. All constitutive relations are stored in + * the global cache indexed by block, and elements reference their block's + * constitutive relation in the DataBox. This means an element can retrieve the + * local constitutive relation from the DataBox simply by requesting + * `Elasticity::Tags::ConstitutiveRelation`. + */ +template +struct InitializeConstitutiveRelation { + public: + using const_global_cache_tags = + tmpl::list, + Tags::MaterialBlockGroups>; + using simple_tags = + tmpl::list>; + using compute_tags = tmpl::list>; + + template + static Parallel::iterable_action_return_t apply( + db::DataBox& box, + const tuples::TaggedTuple& /*inboxes*/, + const Parallel::GlobalCache& /*cache*/, + const ElementId& /*array_index*/, const ActionList /*meta*/, + const ParallelComponent* const /*meta*/) { + db::mutate_apply(make_not_null(&box)); + return {Parallel::AlgorithmExecution::Continue, std::nullopt}; + } + + public: + using return_tags = simple_tags; + using argument_tags = + tmpl::list, domain::Tags::Element, + domain::Tags::Domain>; + + static void apply( + const gsl::not_null*> material_layer_name, + const gsl::not_null*> observation_key, + const std::unordered_set& material_block_groups, + const Element& element, const Domain& domain) { + const auto& block = domain.blocks()[element.id().block_id()]; + // Check if this element is in a material layer + *material_layer_name = [&material_block_groups, &domain, + &block]() -> std::optional { + for (const auto& name : material_block_groups) { + if (domain::block_is_in_group(block.name(), name, + domain.block_groups())) { + return name; + } + } + return std::nullopt; + }(); + // Set the corresponding observation key, but only on the finest multigrid + // level. This could be done better by supporting intersections of array + // sections in observation events or something like that. + *observation_key = + element.id().grid_index() == 0 ? *material_layer_name : std::nullopt; + } +}; + +} // namespace Actions +} // namespace Elasticity diff --git a/src/Elliptic/Systems/Elasticity/CMakeLists.txt b/src/Elliptic/Systems/Elasticity/CMakeLists.txt index 395c9fbf2ee3..49febab496ec 100644 --- a/src/Elliptic/Systems/Elasticity/CMakeLists.txt +++ b/src/Elliptic/Systems/Elasticity/CMakeLists.txt @@ -31,4 +31,5 @@ target_link_libraries( LinearOperators ) +add_subdirectory(Actions) add_subdirectory(BoundaryConditions) diff --git a/src/Elliptic/Systems/Elasticity/Equations.cpp b/src/Elliptic/Systems/Elasticity/Equations.cpp index ec16207ca583..45100ca9e54e 100644 --- a/src/Elliptic/Systems/Elasticity/Equations.cpp +++ b/src/Elliptic/Systems/Elasticity/Equations.cpp @@ -5,6 +5,8 @@ #include #include +#include +#include #include "DataStructures/DataBox/PrefixHelpers.hpp" #include "DataStructures/DataBox/Prefixes.hpp" @@ -57,21 +59,24 @@ void add_curved_sources( template void Fluxes::apply( const gsl::not_null*> minus_stress, - const ConstitutiveRelations::ConstitutiveRelation& - constitutive_relation, - const tnsr::I& coordinates, + const std::vector< + std::unique_ptr>>& + constitutive_relation_per_block, + const Element& element, const tnsr::I& coordinates, const tnsr::I& /*displacement*/, const tnsr::iJ& deriv_displacement) { - primal_fluxes(minus_stress, deriv_displacement, constitutive_relation, + primal_fluxes(minus_stress, deriv_displacement, + *constitutive_relation_per_block.at(element.id().block_id()), coordinates); } template void Fluxes::apply( const gsl::not_null*> minus_stress, - const ConstitutiveRelations::ConstitutiveRelation& - constitutive_relation, - const tnsr::I& coordinates, + const std::vector< + std::unique_ptr>>& + constitutive_relation_per_block, + const Element& element, const tnsr::I& coordinates, const tnsr::i& face_normal, const tnsr::I& /*face_normal_vector*/, const tnsr::I& displacement) { @@ -82,6 +87,8 @@ void Fluxes::apply( face_normal.get(j) * displacement.get(i)); } } + const auto& constitutive_relation = + *constitutive_relation_per_block.at(element.id().block_id()); constitutive_relation.stress(minus_stress, strain, coordinates); for (auto& component : *minus_stress) { component *= -1.; diff --git a/src/Elliptic/Systems/Elasticity/Equations.hpp b/src/Elliptic/Systems/Elasticity/Equations.hpp index 6ab8277794b6..0f71b5c6b5db 100644 --- a/src/Elliptic/Systems/Elasticity/Equations.hpp +++ b/src/Elliptic/Systems/Elasticity/Equations.hpp @@ -4,8 +4,11 @@ #pragma once #include +#include +#include #include "DataStructures/Tensor/Tensor.hpp" +#include "Domain/Structure/Element.hpp" #include "Domain/Tags.hpp" #include "Elliptic/Systems/Elasticity/Tags.hpp" #include "PointwiseFunctions/Elasticity/ConstitutiveRelations/Tags.hpp" @@ -63,23 +66,30 @@ void add_curved_sources( template struct Fluxes { using argument_tags = - tmpl::list, + tmpl::list, domain::Tags::Coordinates>; - using volume_tags = tmpl::list>; - using const_global_cache_tags = volume_tags; - static void apply(gsl::not_null*> minus_stress, - const ConstitutiveRelations::ConstitutiveRelation& - constitutive_relation, - const tnsr::I& coordinates, - const tnsr::I& displacement, - const tnsr::iJ& deriv_displacement); - static void apply(gsl::not_null*> minus_stress, - const ConstitutiveRelations::ConstitutiveRelation& - constitutive_relation, - const tnsr::I& coordinates, - const tnsr::i& face_normal, - const tnsr::I& face_normal_vector, - const tnsr::I& displacement); + using volume_tags = tmpl::list>; + using const_global_cache_tags = + tmpl::list; + static void apply( + gsl::not_null*> minus_stress, + const std::vector< + std::unique_ptr>>& + constitutive_relation_per_block, + const Element& element, const tnsr::I& coordinates, + const tnsr::I& displacement, + const tnsr::iJ& deriv_displacement); + static void apply( + gsl::not_null*> minus_stress, + const std::vector< + std::unique_ptr>>& + constitutive_relation_per_block, + const Element& element, const tnsr::I& coordinates, + const tnsr::i& face_normal, + const tnsr::I& face_normal_vector, + const tnsr::I& displacement); }; } // namespace Elasticity diff --git a/src/Elliptic/Utilities/ApplyAt.hpp b/src/Elliptic/Utilities/ApplyAt.hpp index 4f92f9844267..8f3ce616c7bd 100644 --- a/src/Elliptic/Utilities/ApplyAt.hpp +++ b/src/Elliptic/Utilities/ApplyAt.hpp @@ -16,6 +16,7 @@ #include "Utilities/TaggedTuple.hpp" #include "Utilities/TupleSlice.hpp" #include "Utilities/TypeTraits/CreateIsCallable.hpp" +#include "Utilities/TypeTraits/IsMaplike.hpp" namespace elliptic::util { namespace detail { @@ -60,18 +61,24 @@ struct unmap_arg { template static constexpr decltype(auto) apply( const T& arg, const std::tuple& map_keys) { - return unmap_arg::apply( - arg.at(std::get<0>(map_keys)), - tuple_tail(map_keys)); + return unmap_arg<((sizeof...(MapKeys) == 0) or + (not tt::is_maplike_v)), + MapKeys...>::apply(arg.at(std::get<0>(map_keys)), + tuple_tail( + map_keys)); } template static constexpr decltype(auto) apply( const gsl::not_null arg, const std::tuple& map_keys) { - return unmap_arg::apply( - make_not_null(&arg->operator[](std::get<0>(map_keys))), - tuple_tail(map_keys)); + return unmap_arg<((sizeof...(MapKeys) == 0) or + (not tt::is_maplike_v)), + MapKeys...>::apply(make_not_null(&arg-> + operator[](std::get<0>( + map_keys))), + tuple_tail( + map_keys)); } }; diff --git a/src/Evolution/DgSubcell/CMakeLists.txt b/src/Evolution/DgSubcell/CMakeLists.txt index 240e6980a235..cdb6fa5bfed9 100644 --- a/src/Evolution/DgSubcell/CMakeLists.txt +++ b/src/Evolution/DgSubcell/CMakeLists.txt @@ -68,6 +68,7 @@ target_link_libraries( PUBLIC DataStructures Domain + DomainCreators Interpolation Parallel Spectral diff --git a/src/Evolution/DgSubcell/SubcellOptions.cpp b/src/Evolution/DgSubcell/SubcellOptions.cpp index 92136127896b..5177ec8db0db 100644 --- a/src/Evolution/DgSubcell/SubcellOptions.cpp +++ b/src/Evolution/DgSubcell/SubcellOptions.cpp @@ -12,13 +12,13 @@ #include #include +#include "Domain/Creators/BlockGroups.hpp" #include "Domain/Creators/DomainCreator.hpp" #include "Options/Options.hpp" #include "Utilities/Algorithm.hpp" #include "Utilities/ErrorHandling/Error.hpp" #include "Utilities/GenerateInstantiations.hpp" #include "Utilities/Serialization/PupStlCpp17.hpp" -#include "Utilities/StdHelpers.hpp" namespace evolution::dg::subcell { SubcellOptions::SubcellOptions( @@ -49,41 +49,17 @@ SubcellOptions::SubcellOptions( const auto& only_dg_block_and_group_names = subcell_options_with_block_names.only_dg_block_and_group_names_; - const auto block_names_vector = domain_creator.block_names(); - const auto group_names = domain_creator.block_groups(); - std::unordered_set block_names{}; - // Add blocks from block groups into block_names. Use an unordered_set - // to elide duplicates. - for (const auto& block_or_group_name : - only_dg_block_and_group_names.value_or(std::vector{})) { - if (const auto block_group_it = group_names.find(block_or_group_name); - block_group_it != group_names.end()) { - block_names.insert(block_group_it->second.begin(), - block_group_it->second.end()); - } else if (const auto block_name_it = - std::find(block_names_vector.begin(), - block_names_vector.end(), block_or_group_name); - block_name_it != block_names_vector.end()) { - block_names.insert(*block_name_it); - } else { - using ::operator<<; - ERROR_NO_TRACE( - "The block or group '" - << block_or_group_name - << "' is not one of the block names or groups of the domain. " - "The known block groups are " - << keys_of(group_names) << " and the known block names are " - << block_names_vector); - } - } + const auto block_names = domain_creator.block_names(); + const auto only_dg_block_names = domain::expand_block_groups_to_block_names( + only_dg_block_and_group_names.value_or(std::vector{}), + block_names, domain_creator.block_groups()); only_dg_block_ids_ = std::vector{}; - only_dg_block_ids_.value().reserve(block_names.size()); + only_dg_block_ids_.value().reserve(only_dg_block_names.size()); // Get the block ID of each block name - for (const auto& block_name : block_names) { - only_dg_block_ids_.value().push_back(static_cast( - std::distance(block_names_vector.begin(), - std::find(block_names_vector.begin(), - block_names_vector.end(), block_name)))); + for (const auto& block_name : only_dg_block_names) { + only_dg_block_ids_.value().push_back(static_cast(std::distance( + block_names.begin(), + std::find(block_names.begin(), block_names.end(), block_name)))); } // Sort the block IDs just so they're easier to deal with. alg::sort(only_dg_block_ids_.value()); diff --git a/src/NumericalAlgorithms/LinearOperators/FilterAction.hpp b/src/NumericalAlgorithms/LinearOperators/FilterAction.hpp index 9e452adbe0ba..3b4dd16682c9 100644 --- a/src/NumericalAlgorithms/LinearOperators/FilterAction.hpp +++ b/src/NumericalAlgorithms/LinearOperators/FilterAction.hpp @@ -13,10 +13,12 @@ #include "DataStructures/DataBox/DataBox.hpp" #include "DataStructures/DataVector.hpp" #include "DataStructures/Matrix.hpp" +#include "Domain/Creators/BlockGroups.hpp" #include "Domain/Tags.hpp" #include "NumericalAlgorithms/Spectral/Mesh.hpp" #include "Parallel/AlgorithmExecution.hpp" #include "Parallel/GlobalCache.hpp" +#include "Utilities/Algorithm.hpp" #include "Utilities/ErrorHandling/Error.hpp" #include "Utilities/Gsl.hpp" #include "Utilities/StdHelpers.hpp" @@ -121,28 +123,13 @@ class Filter> { // Only do this check if filtering is enabled. A `nullopt` means all blocks // are allowed to do filtering if (enable and filter_helper.blocks_to_filter().has_value()) { - const auto& blocks_to_filter = filter_helper.blocks_to_filter().value(); - - // If nothing in the loop sets this to true, then our block isn't in the - // blocks to filter, so we shouldn't filter in this element - bool block_is_in_blocks_to_filter = false; - for (const std::string& block_to_filter : blocks_to_filter) { - // First check if the block to filter is our current block. If it is, - // then we do want to filter. - if (block_to_filter == block_name) { - block_is_in_blocks_to_filter = true; - break; - } else if (block_groups.count(block_to_filter) == 1 and - block_groups.at(block_to_filter).count(block_name) == 1) { - // If the block to filter isn't our current block, then we need to - // check if the block to filter is a block group. If it is, check that - // our block is in the block group. If it is, then we enable filtering - block_is_in_blocks_to_filter = true; - break; - } - } - - enable = block_is_in_blocks_to_filter; + // Enable filtering for this block if it's in any of the listed groups + enable = alg::any_of( + filter_helper.blocks_to_filter().value(), + [&block_name, &block_groups](const std::string& block_to_filter) { + return domain::block_is_in_group(block_name, block_to_filter, + block_groups); + }); } if (not enable) { diff --git a/src/ParallelAlgorithms/Events/ObserveNorms.hpp b/src/ParallelAlgorithms/Events/ObserveNorms.hpp index 62b7663892bc..0f6f18062fa8 100644 --- a/src/ParallelAlgorithms/Events/ObserveNorms.hpp +++ b/src/ParallelAlgorithms/Events/ObserveNorms.hpp @@ -38,6 +38,7 @@ #include "Utilities/ErrorHandling/Error.hpp" #include "Utilities/Functional.hpp" #include "Utilities/OptionalHelpers.hpp" +#include "Utilities/PrettyType.hpp" #include "Utilities/Serialization/CharmPupable.hpp" #include "Utilities/TMPL.hpp" @@ -88,17 +89,24 @@ namespace Events { * of elements. The `observers::Tags::ObservationKey` must be * available in the DataBox. It identifies the section and is used as a suffix * for the path in the output file. + * + * \par Option name + * The `OptionName` template parameter is used to give the event a name in the + * input file. If it is not specified, the name defaults to "ObserveNorms". If + * you have multiple `ObserveNorms` events in the input file, you must specify a + * unique name for each one. This can happen, for example, if you want to + * observe norms the full domain and also over a section of the domain. */ template , - typename ArraySectionIdTag = void> + typename ArraySectionIdTag = void, typename OptionName = void> class ObserveNorms; template + typename ArraySectionIdTag, typename OptionName> class ObserveNorms, - tmpl::list, ArraySectionIdTag> - : public Event { + tmpl::list, ArraySectionIdTag, + OptionName> : public Event { private: struct ObserveTensor { static constexpr Options::String help = { @@ -162,6 +170,14 @@ class ObserveNorms, funcl::ElementWise>>>; public: + static std::string name() { + if constexpr (std::is_same_v) { + return "ObserveNorms"; + } else { + return pretty_type::name(); + } + } + /// The name of the subfile inside the HDF5 file struct SubfileName { using type = std::string; @@ -276,19 +292,19 @@ class ObserveNorms, /// \cond template + typename ArraySectionIdTag, typename OptionName> ObserveNorms, - tmpl::list, - ArraySectionIdTag>::ObserveNorms(CkMigrateMessage* msg) + tmpl::list, ArraySectionIdTag, + OptionName>::ObserveNorms(CkMigrateMessage* msg) : Event(msg) {} template + typename ArraySectionIdTag, typename OptionName> ObserveNorms, - tmpl::list, - ArraySectionIdTag>::ObserveNorms(const std::string& subfile_name, - const std::vector& - observe_tensors) + tmpl::list, ArraySectionIdTag, + OptionName>::ObserveNorms(const std::string& subfile_name, + const std::vector& + observe_tensors) : subfile_path_("/" + subfile_name) { tensor_names_.reserve(observe_tensors.size()); tensor_norm_types_.reserve(observe_tensors.size()); @@ -301,13 +317,14 @@ ObserveNorms, } template -ObserveNorms, - tmpl::list, ArraySectionIdTag>:: - ObserveTensor::ObserveTensor(std::string in_tensor, - std::string in_norm_type, - std::string in_components, - const Options::Context& context) + typename ArraySectionIdTag, typename OptionName> +ObserveNorms< + tmpl::list, tmpl::list, + ArraySectionIdTag, + OptionName>::ObserveTensor::ObserveTensor(std::string in_tensor, + std::string in_norm_type, + std::string in_components, + const Options::Context& context) : tensor(std::move(in_tensor)), norm_type(std::move(in_norm_type)), components(std::move(in_components)) { @@ -349,11 +366,12 @@ void fill_norm_values_and_names( } // namespace ObserveNorms_impl template + typename ArraySectionIdTag, typename OptionName> template void ObserveNorms, - tmpl::list, ArraySectionIdTag>:: + tmpl::list, ArraySectionIdTag, + OptionName>:: observe_norms_impl( const gsl::not_null, } template + typename ArraySectionIdTag, typename OptionName> template void ObserveNorms, - tmpl::list, ArraySectionIdTag>:: + tmpl::list, ArraySectionIdTag, + OptionName>:: operator()(const ObservationBox& box, Parallel::GlobalCache& cache, const ElementId& array_index, @@ -449,10 +468,10 @@ operator()(const ObservationBox& box, } template + typename ArraySectionIdTag, typename OptionName> void ObserveNorms, - tmpl::list, - ArraySectionIdTag>::pup(PUP::er& p) { + tmpl::list, ArraySectionIdTag, + OptionName>::pup(PUP::er& p) { Event::pup(p); p | subfile_path_; p | tensor_names_; @@ -460,10 +479,12 @@ void ObserveNorms, p | tensor_components_; } +// NOLINTBEGIN(cppcoreguidelines-avoid-non-const-global-variables) template + typename ArraySectionIdTag, typename OptionName> PUP::able::PUP_ID ObserveNorms, tmpl::list, - ArraySectionIdTag>::my_PUP_ID = 0; // NOLINT + ArraySectionIdTag, OptionName>::my_PUP_ID = 0; +// NOLINTEND(cppcoreguidelines-avoid-non-const-global-variables) /// \endcond } // namespace Events diff --git a/src/PointwiseFunctions/Elasticity/ConstitutiveRelations/ConstitutiveRelation.hpp b/src/PointwiseFunctions/Elasticity/ConstitutiveRelations/ConstitutiveRelation.hpp index 259a4bc673dc..c2857e86d22e 100644 --- a/src/PointwiseFunctions/Elasticity/ConstitutiveRelations/ConstitutiveRelation.hpp +++ b/src/PointwiseFunctions/Elasticity/ConstitutiveRelations/ConstitutiveRelation.hpp @@ -4,6 +4,7 @@ #pragma once #include +#include #include #include "DataStructures/Tensor/TypeAliases.hpp" @@ -57,6 +58,10 @@ class ConstitutiveRelation : public PUP::able { WRAPPED_PUPable_abstract(ConstitutiveRelation); // NOLINT + /// Returns a `std::unique_ptr` pointing to a copy of the + /// `ConstitutiveRelation`. + virtual std::unique_ptr> get_clone() const = 0; + /// The constitutive relation that characterizes the elastic properties of a /// material virtual void stress(gsl::not_null*> stress, diff --git a/src/PointwiseFunctions/Elasticity/ConstitutiveRelations/CubicCrystal.cpp b/src/PointwiseFunctions/Elasticity/ConstitutiveRelations/CubicCrystal.cpp index f5c5f1ce3d4a..90264a961937 100644 --- a/src/PointwiseFunctions/Elasticity/ConstitutiveRelations/CubicCrystal.cpp +++ b/src/PointwiseFunctions/Elasticity/ConstitutiveRelations/CubicCrystal.cpp @@ -4,11 +4,13 @@ #include "PointwiseFunctions/Elasticity/ConstitutiveRelations/CubicCrystal.hpp" #include -#include // IWYU pragma: keep +#include +#include #include "DataStructures/DataVector.hpp" -#include "DataStructures/Tensor/Tensor.hpp" // IWYU pragma: keep +#include "DataStructures/Tensor/Tensor.hpp" #include "DataStructures/Tensor/TypeAliases.hpp" +#include "PointwiseFunctions/Elasticity/ConstitutiveRelations/ConstitutiveRelation.hpp" #include "Utilities/ErrorHandling/Assert.hpp" #include "Utilities/MakeWithValue.hpp" @@ -27,6 +29,10 @@ CubicCrystal::CubicCrystal(const double c_11, const double c_12, "must be positive and the poisson ratio smaller or equal to 0.5."); } +std::unique_ptr> CubicCrystal::get_clone() const { + return std::make_unique(*this); +} + void CubicCrystal::stress(const gsl::not_null*> stress, const tnsr::ii& strain, const tnsr::I& /*x*/) const { diff --git a/src/PointwiseFunctions/Elasticity/ConstitutiveRelations/CubicCrystal.hpp b/src/PointwiseFunctions/Elasticity/ConstitutiveRelations/CubicCrystal.hpp index 48be1691c7b2..a65385ae87c4 100644 --- a/src/PointwiseFunctions/Elasticity/ConstitutiveRelations/CubicCrystal.hpp +++ b/src/PointwiseFunctions/Elasticity/ConstitutiveRelations/CubicCrystal.hpp @@ -5,11 +5,12 @@ #include #include +#include #include #include "DataStructures/Tensor/TypeAliases.hpp" #include "Options/String.hpp" -#include "PointwiseFunctions/Elasticity/ConstitutiveRelations/ConstitutiveRelation.hpp" // IWYU pragma: keep +#include "PointwiseFunctions/Elasticity/ConstitutiveRelations/ConstitutiveRelation.hpp" #include "Utilities/Serialization/CharmPupable.hpp" #include "Utilities/TMPL.hpp" @@ -111,6 +112,8 @@ class CubicCrystal : public ConstitutiveRelation<3> { CubicCrystal(double c_11, double c_12, double c_44); + std::unique_ptr> get_clone() const override; + /// The constitutive relation that characterizes the elastic properties of a /// material void stress(gsl::not_null*> stress, diff --git a/src/PointwiseFunctions/Elasticity/ConstitutiveRelations/IsotropicHomogeneous.cpp b/src/PointwiseFunctions/Elasticity/ConstitutiveRelations/IsotropicHomogeneous.cpp index fa464b1e03da..0406d5a478b1 100644 --- a/src/PointwiseFunctions/Elasticity/ConstitutiveRelations/IsotropicHomogeneous.cpp +++ b/src/PointwiseFunctions/Elasticity/ConstitutiveRelations/IsotropicHomogeneous.cpp @@ -4,11 +4,13 @@ #include "PointwiseFunctions/Elasticity/ConstitutiveRelations/IsotropicHomogeneous.hpp" #include -#include // IWYU pragma: keep +#include +#include #include "DataStructures/DataVector.hpp" -#include "DataStructures/Tensor/Tensor.hpp" // IWYU pragma: keep +#include "DataStructures/Tensor/Tensor.hpp" #include "DataStructures/Tensor/TypeAliases.hpp" +#include "PointwiseFunctions/Elasticity/ConstitutiveRelations/ConstitutiveRelation.hpp" #include "Utilities/GenerateInstantiations.hpp" #include "Utilities/Gsl.hpp" #include "Utilities/MakeWithValue.hpp" @@ -20,6 +22,12 @@ IsotropicHomogeneous::IsotropicHomogeneous(double bulk_modulus, double shear_modulus) : bulk_modulus_(bulk_modulus), shear_modulus_(shear_modulus) {} +template +std::unique_ptr> +IsotropicHomogeneous::get_clone() const { + return std::make_unique(*this); +} + template <> void IsotropicHomogeneous<3>::stress( const gsl::not_null*> stress, diff --git a/src/PointwiseFunctions/Elasticity/ConstitutiveRelations/IsotropicHomogeneous.hpp b/src/PointwiseFunctions/Elasticity/ConstitutiveRelations/IsotropicHomogeneous.hpp index ac6ec5c63bd3..45031db56784 100644 --- a/src/PointwiseFunctions/Elasticity/ConstitutiveRelations/IsotropicHomogeneous.hpp +++ b/src/PointwiseFunctions/Elasticity/ConstitutiveRelations/IsotropicHomogeneous.hpp @@ -5,11 +5,12 @@ #include #include +#include #include #include "DataStructures/Tensor/TypeAliases.hpp" #include "Options/String.hpp" -#include "PointwiseFunctions/Elasticity/ConstitutiveRelations/ConstitutiveRelation.hpp" // IWYU pragma: keep +#include "PointwiseFunctions/Elasticity/ConstitutiveRelations/ConstitutiveRelation.hpp" #include "Utilities/Gsl.hpp" #include "Utilities/Serialization/CharmPupable.hpp" #include "Utilities/TMPL.hpp" @@ -127,6 +128,8 @@ class IsotropicHomogeneous : public ConstitutiveRelation { IsotropicHomogeneous(double bulk_modulus, double shear_modulus); + std::unique_ptr> get_clone() const override; + /// The constitutive relation that characterizes the elastic properties of a /// material void stress(gsl::not_null*> stress, diff --git a/src/PointwiseFunctions/Elasticity/ConstitutiveRelations/Tags.hpp b/src/PointwiseFunctions/Elasticity/ConstitutiveRelations/Tags.hpp index be5ff0116891..382c2c540a4c 100644 --- a/src/PointwiseFunctions/Elasticity/ConstitutiveRelations/Tags.hpp +++ b/src/PointwiseFunctions/Elasticity/ConstitutiveRelations/Tags.hpp @@ -6,6 +6,9 @@ #include #include #include +#include +#include +#include #include "DataStructures/DataBox/Tag.hpp" #include "Options/String.hpp" @@ -24,6 +27,17 @@ struct ConstitutiveRelation : db::SimpleTag { using type = std::unique_ptr>; }; + +template +struct ConstitutiveRelationPerBlock { + static std::string name() { return "Material"; } + using ConstRelPtr = + std::unique_ptr>; + using type = std::variant, + std::unordered_map>; + static constexpr Options::String help = + "A constitutive relation in every block of the domain."; +}; } // namespace OptionTags namespace Tags { @@ -45,5 +59,8 @@ struct ConstitutiveRelation : db::SimpleTag { } }; +/// A constitutive relation in every block of the domain +struct ConstitutiveRelationPerBlockBase : db::BaseTag {}; + } // namespace Tags } // namespace Elasticity diff --git a/tests/InputFiles/Elasticity/HalfSpaceMirror.yaml b/tests/InputFiles/Elasticity/HalfSpaceMirror.yaml index 7671cb70337f..4fd88d6d7069 100644 --- a/tests/InputFiles/Elasticity/HalfSpaceMirror.yaml +++ b/tests/InputFiles/Elasticity/HalfSpaceMirror.yaml @@ -36,7 +36,10 @@ Background: &solution RelativeTolerance: 1e-10 Material: - IsotropicHomogeneous: *fused_silica + Layer0: + IsotropicHomogeneous: *fused_silica + Layer1: + IsotropicHomogeneous: *fused_silica InitialGuess: Zero @@ -52,9 +55,9 @@ DomainCreator: InitialGridPoints: [2, 3, 3] UseEquiangularMap: True RadialPartitioning: [] - PartitioningInZ: [] + PartitioningInZ: [0.1] RadialDistribution: [Linear] - DistributionInZ: [Linear] + DistributionInZ: [Linear, Linear] BoundaryConditions: LowerZ: AnalyticSolution: @@ -83,23 +86,23 @@ Observers: LinearSolver: Gmres: ConvergenceCriteria: - MaxIterations: 1 + MaxIterations: 2 RelativeResidual: 1.e-4 AbsoluteResidual: 1.e-12 - Verbosity: Verbose + Verbosity: Quiet Multigrid: Iterations: 1 MaxLevels: Auto PreSmoothing: True PostSmoothingAtBottom: True - Verbosity: Quiet + Verbosity: Silent OutputVolumeData: False SchwarzSmoother: Iterations: 3 MaxOverlap: 2 - Verbosity: Verbose + Verbosity: Silent SubdomainSolver: Gmres: ConvergenceCriteria: @@ -119,18 +122,20 @@ LinearSolver: EventsAndTriggers: - Trigger: HasConverged Events: - - ObserveNorms: + - ObserveNorms: &error_norms SubfileName: ErrorNorms TensorsToObserve: - Name: Error(Displacement) NormType: L2Norm Components: Sum - - ObserveNorms: + - ObserveNorms: &volume_integrals SubfileName: VolumeIntegrals TensorsToObserve: - Name: PotentialEnergyDensity NormType: VolumeIntegral Components: Individual + - ObserveNormsPerLayer: *error_norms + - ObserveNormsPerLayer: *volume_integrals - ObserveFields: SubfileName: VolumeData VariablesToObserve: diff --git a/tests/Unit/Domain/Creators/CMakeLists.txt b/tests/Unit/Domain/Creators/CMakeLists.txt index af910feff237..503033444f6d 100644 --- a/tests/Unit/Domain/Creators/CMakeLists.txt +++ b/tests/Unit/Domain/Creators/CMakeLists.txt @@ -10,6 +10,7 @@ set(LIBRARY_SOURCES Test_AlignedLattice.cpp Test_BinaryCompactObject.cpp Test_BinaryCompactObjectHelpers.cpp + Test_BlockGroups.cpp Test_Brick.cpp Test_Cylinder.cpp Test_CylindricalBinaryCompactObject.cpp diff --git a/tests/Unit/Domain/Creators/Test_BlockGroups.cpp b/tests/Unit/Domain/Creators/Test_BlockGroups.cpp new file mode 100644 index 000000000000..7b1b58f56a76 --- /dev/null +++ b/tests/Unit/Domain/Creators/Test_BlockGroups.cpp @@ -0,0 +1,41 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Framework/TestingFramework.hpp" + +#include +#include +#include +#include + +#include "Domain/Creators/BlockGroups.hpp" + +namespace domain { + +SPECTRE_TEST_CASE("Unit.Domain.BlockGroups", "[Domain][Unit]") { + const std::unordered_map> + block_groups{ + {"Group1", {"Block1", "Block2"}}, + {"Group2", {"Block3", "Block4"}}, + }; + CHECK(block_is_in_group("Block1", "Group1", block_groups)); + CHECK(block_is_in_group("Block2", "Group1", block_groups)); + CHECK(not block_is_in_group("Block3", "Group1", block_groups)); + CHECK(not block_is_in_group("Block1", "Group2", block_groups)); + + const std::vector all_block_names{"Block1", "Block2", "Block3", + "Block4"}; + const std::unordered_set expanded_block_names{"Block1", "Block3", + "Block4"}; + CHECK(expand_block_groups_to_block_names({"Block1", "Group2"}, + all_block_names, block_groups) == + expanded_block_names); + CHECK_THROWS_WITH( + expand_block_groups_to_block_names({"NoBlock", "Group2"}, all_block_names, + block_groups), + Catch::Matchers::ContainsSubstring( + "The block or group 'NoBlock' is not one of the block names or " + "groups of the domain.")); +} + +} // namespace domain diff --git a/tests/Unit/Domain/Creators/Test_ExpandOverBlocks.cpp b/tests/Unit/Domain/Creators/Test_ExpandOverBlocks.cpp index 5e1866fcd51e..62670426bdc8 100644 --- a/tests/Unit/Domain/Creators/Test_ExpandOverBlocks.cpp +++ b/tests/Unit/Domain/Creators/Test_ExpandOverBlocks.cpp @@ -5,10 +5,11 @@ #include #include +#include #include #include -#include "Domain/Creators/ExpandOverBlocks.hpp" +#include "Domain/Creators/ExpandOverBlocks.tpp" #include "Utilities/ErrorHandling/Error.hpp" namespace domain { @@ -18,10 +19,23 @@ void test_expand_over_blocks( input_value, const size_t num_blocks, const std::vector>& expected_expanded_value) { - CHECK(std::visit(ExpandOverBlocks{num_blocks}, input_value) == - expected_expanded_value); + CHECK(std::visit(ExpandOverBlocks>{num_blocks}, + input_value) == expected_expanded_value); } +struct Base { + virtual ~Base() = default; + virtual std::unique_ptr get_clone() const = 0; +}; + +struct Derived : Base { + explicit Derived(const size_t local_id) : id{local_id} {} + std::unique_ptr get_clone() const override { + return std::make_unique(*this); + } + size_t id; +}; + SPECTRE_TEST_CASE("Unit.Domain.ExpandOverBlocks", "[Domain][Unit]") { test_expand_over_blocks(size_t{2}, 3, {3, {2}}); test_expand_over_blocks(std::array{2}, 3, {3, {2}}); @@ -40,7 +54,7 @@ SPECTRE_TEST_CASE("Unit.Domain.ExpandOverBlocks", "[Domain][Unit]") { std::vector>{{2, 3, 4}, {3, 4, 5}, {4, 5, 6}}, 3, {{2, 3, 4}, {3, 4, 5}, {4, 5, 6}}); CHECK_THROWS_WITH( - (ExpandOverBlocks{3}( + (ExpandOverBlocks>{3}( std::vector>{{2}, {3}})), Catch::Matchers::ContainsSubstring( "You supplied 2 values, but the domain creator has 3 blocks.")); @@ -59,7 +73,7 @@ SPECTRE_TEST_CASE("Unit.Domain.ExpandOverBlocks", "[Domain][Unit]") { try { // Invoke `ExpandOverBlocks`: const auto initial_refinement = - std::visit(ExpandOverBlocks{num_blocks}, + std::visit(ExpandOverBlocks>{num_blocks}, initial_refinement_from_options); // Since a single number was specified, we expect the vector over blocks // is homogeneously and isotropically filled with that number: @@ -83,7 +97,7 @@ SPECTRE_TEST_CASE("Unit.Domain.ExpandOverBlocks", "[Domain][Unit]") { block_groups{{"Wedges", {"East", "North", "West", "South"}}}; // Now we can expand values over blocks by giving their names. This can also // be used with a std::variant like in the other example. - ExpandOverBlocks expand{block_names, block_groups}; + ExpandOverBlocks> expand{block_names, block_groups}; CHECK(expand({{"West", {{3, 4}}}, {"InnerCube", {{2, 3}}}, {"South", {{3, 4}}}, @@ -113,6 +127,50 @@ SPECTRE_TEST_CASE("Unit.Domain.ExpandOverBlocks", "[Domain][Unit]") { "Value for block 'InnerCube' is missing. Did " "you misspell its name?")); } + { + INFO("Expand non-array"); + ExpandOverBlocks expand{3}; + CHECK(expand(std::string{"A"}) == std::vector{3, "A"}); + } + { + INFO("Expand unique_ptr"); + ExpandOverBlocks> expand{{"A", "B"}}; + { + const auto expanded = expand(std::make_unique(size_t{1})); + CHECK(expanded.size() == 2); + for (const auto& v : expanded) { + const auto v_derived = dynamic_cast(v.get()); + CHECK(v_derived != nullptr); + CHECK(v_derived->id == 1); + } + } + { + std::vector> original{}; + original.emplace_back(std::make_unique(size_t{1})); + original.emplace_back(std::make_unique(size_t{2})); + const auto expanded = expand(original); + CHECK(expanded.size() == 2); + auto v_derived = dynamic_cast(expanded[0].get()); + CHECK(v_derived != nullptr); + CHECK(v_derived->id == 1); + v_derived = dynamic_cast(expanded[1].get()); + CHECK(v_derived != nullptr); + CHECK(v_derived->id == 2); + } + { + std::unordered_map> original{}; + original.emplace("A", std::make_unique(size_t{1})); + original.emplace("B", std::make_unique(size_t{2})); + const auto expanded = expand(original); + CHECK(expanded.size() == 2); + auto v_derived = dynamic_cast(expanded[0].get()); + CHECK(v_derived != nullptr); + CHECK(v_derived->id == 1); + v_derived = dynamic_cast(expanded[1].get()); + CHECK(v_derived != nullptr); + CHECK(v_derived->id == 2); + } + } } } // namespace domain diff --git a/tests/Unit/Domain/Test_Domain.cpp b/tests/Unit/Domain/Test_Domain.cpp index 922b833dc13d..41f0404615a9 100644 --- a/tests/Unit/Domain/Test_Domain.cpp +++ b/tests/Unit/Domain/Test_Domain.cpp @@ -519,8 +519,8 @@ Domain<3> create_serialized_domain() { add_outer_region("OuterShell"); // 10 blocks // Expand initial refinement and number of grid points over all blocks - const ExpandOverBlocks expand_over_blocks{block_names, - block_groups}; + const ExpandOverBlocks> expand_over_blocks{ + block_names, block_groups}; using BCO = domain::creators::BinaryCompactObject; diff --git a/tests/Unit/Elliptic/DiscontinuousGalerkin/SubdomainOperator/CMakeLists.txt b/tests/Unit/Elliptic/DiscontinuousGalerkin/SubdomainOperator/CMakeLists.txt index d38b73805c12..f7c7ab90297a 100644 --- a/tests/Unit/Elliptic/DiscontinuousGalerkin/SubdomainOperator/CMakeLists.txt +++ b/tests/Unit/Elliptic/DiscontinuousGalerkin/SubdomainOperator/CMakeLists.txt @@ -19,6 +19,7 @@ target_link_libraries( DomainCreators DomainStructure Elasticity + ElasticityActions ElasticityBoundaryConditions Elliptic EllipticDg diff --git a/tests/Unit/Elliptic/DiscontinuousGalerkin/SubdomainOperator/Test_SubdomainOperator.cpp b/tests/Unit/Elliptic/DiscontinuousGalerkin/SubdomainOperator/Test_SubdomainOperator.cpp index 7ed43e94e1f8..e522bbee5947 100644 --- a/tests/Unit/Elliptic/DiscontinuousGalerkin/SubdomainOperator/Test_SubdomainOperator.cpp +++ b/tests/Unit/Elliptic/DiscontinuousGalerkin/SubdomainOperator/Test_SubdomainOperator.cpp @@ -24,6 +24,7 @@ #include "Domain/Creators/Cylinder.hpp" #include "Domain/Creators/Disk.hpp" #include "Domain/Creators/DomainCreator.hpp" +#include "Domain/Creators/ExpandOverBlocks.hpp" #include "Domain/Creators/Interval.hpp" #include "Domain/Creators/Rectangle.hpp" #include "Domain/Creators/RegisterDerivedWithCharm.hpp" @@ -47,6 +48,7 @@ #include "Elliptic/DiscontinuousGalerkin/SubdomainOperator/InitializeSubdomain.hpp" #include "Elliptic/DiscontinuousGalerkin/SubdomainOperator/SubdomainOperator.hpp" #include "Elliptic/DiscontinuousGalerkin/SubdomainOperator/Tags.hpp" +#include "Elliptic/Systems/Elasticity/Actions/InitializeConstitutiveRelation.hpp" #include "Elliptic/Systems/Elasticity/FirstOrderSystem.hpp" #include "Elliptic/Systems/Poisson/FirstOrderSystem.hpp" #include "Elliptic/Systems/Poisson/Geometry.hpp" @@ -801,7 +803,8 @@ template struct InitializeConstitutiveRelation : tt::ConformsTo<::amr::protocols::Projector> { public: // Iterative action - using simple_tags = tmpl::list>; + using simple_tags = + tmpl::list>; template & /*cache*/, const ArrayIndex& /*array_index*/, const ActionList /*meta*/, const ParallelComponent* const /*meta*/) { - ::Initialization::mutate_assign( - make_not_null(&box), - std::make_unique< - Elasticity::ConstitutiveRelations::IsotropicHomogeneous>(1., - 2.)); + db::mutate_apply(make_not_null(&box)); return {Parallel::AlgorithmExecution::Continue, std::nullopt}; } public: // amr::protocols::Projector - using return_tags = tmpl::list>; - using argument_tags = tmpl::list<>; + using return_tags = simple_tags; + using argument_tags = tmpl::list>; template static void apply( - const gsl::not_null>*> - material, - const AmrData&... /*amr_data*/) { - *material = std::make_unique< - Elasticity::ConstitutiveRelations::IsotropicHomogeneous>(1., 2.); + const gsl::not_null>>*> + constitutive_relations, + const Domain& domain, const AmrData&... /*amr_data*/) { + const domain::ExpandOverBlocks>> + expand_over_blocks{domain.blocks().size()}; + *constitutive_relations = expand_over_blocks( + std::make_unique< + Elasticity::ConstitutiveRelations::IsotropicHomogeneous>(1., + 2.)); } }; diff --git a/tests/Unit/Elliptic/Systems/Elasticity/Actions/CMakeLists.txt b/tests/Unit/Elliptic/Systems/Elasticity/Actions/CMakeLists.txt new file mode 100644 index 000000000000..f149e43b11a4 --- /dev/null +++ b/tests/Unit/Elliptic/Systems/Elasticity/Actions/CMakeLists.txt @@ -0,0 +1,22 @@ +# Distributed under the MIT License. +# See LICENSE.txt for details. + +set(LIBRARY "Test_ElasticityActions") + +set(LIBRARY_SOURCES + Test_InitializeConstitutiveRelation.cpp + ) + +add_test_library(${LIBRARY} "${LIBRARY_SOURCES}") + +target_link_libraries( + ${LIBRARY} + PRIVATE + ConstitutiveRelations + DataStructures + Domain + ElasticityActions + Elliptic + Elasticity + Utilities + ) diff --git a/tests/Unit/Elliptic/Systems/Elasticity/Actions/Test_InitializeConstitutiveRelation.cpp b/tests/Unit/Elliptic/Systems/Elasticity/Actions/Test_InitializeConstitutiveRelation.cpp new file mode 100644 index 000000000000..559ce91542aa --- /dev/null +++ b/tests/Unit/Elliptic/Systems/Elasticity/Actions/Test_InitializeConstitutiveRelation.cpp @@ -0,0 +1,134 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Framework/TestingFramework.hpp" + +#include +#include +#include +#include +#include +#include + +#include "DataStructures/DataBox/DataBox.hpp" +#include "Domain/Creators/Cylinder.hpp" +#include "Domain/Creators/RegisterDerivedWithCharm.hpp" +#include "Domain/Structure/Element.hpp" +#include "Domain/Structure/ElementId.hpp" +#include "Domain/Tags.hpp" +#include "Elliptic/Systems/Elasticity/Actions/InitializeConstitutiveRelation.hpp" +#include "Framework/ActionTesting.hpp" +#include "Options/Protocols/FactoryCreation.hpp" +#include "Parallel/Phase.hpp" +#include "PointwiseFunctions/Elasticity/ConstitutiveRelations/Factory.hpp" +#include "PointwiseFunctions/Elasticity/ConstitutiveRelations/Tags.hpp" +#include "Utilities/ProtocolHelpers.hpp" +#include "Utilities/Serialization/RegisterDerivedClassesWithCharm.hpp" +#include "Utilities/TMPL.hpp" + +namespace Elasticity { +namespace { +constexpr size_t Dim = 3; + +template +struct ElementArray { + using metavariables = Metavariables; + using chare_type = ActionTesting::MockArrayChare; + using array_index = ElementId; + using const_global_cache_tags = tmpl::list<>; + using phase_dependent_action_list = tmpl::list< + Parallel::PhaseActions>>>>, + Parallel::PhaseActions< + Parallel::Phase::Testing, + tmpl::list>>>; +}; + +struct Metavariables { + using component_list = tmpl::list>; + using const_global_cache_tags = tmpl::list>; + struct factory_creation + : tt::ConformsTo { + using factory_classes = tmpl::map, + ConstitutiveRelations::standard_constitutive_relations>>; + }; +}; + +} // namespace + +SPECTRE_TEST_CASE("Unit.Elasticity.Actions.InitializeConstitutiveRelation", + "[Unit][Elliptic][Actions]") { + domain::creators::register_derived_with_charm(); + register_factory_classes_with_charm(); + // A cylinder with two layers in z-direction + const std::unique_ptr> domain_creator = + std::make_unique( + 1., 3., 0., 10., false, 1_st, 3_st, false, std::vector{}, + std::vector{2.}, + std::vector{ + domain::CoordinateMaps::Distribution::Linear}, + std::vector{ + domain::CoordinateMaps::Distribution::Linear, + domain::CoordinateMaps::Distribution::Linear}); + const ElementId element_id_layer1{0, {{{1, 0}, {1, 0}, {1, 0}}}}; + const ElementId element_id_layer2{5, {{{1, 0}, {1, 0}, {1, 0}}}}; + // A different material in each layer + using ConstRelPtr = + std::unique_ptr>; + std::unordered_map material_layers{}; + material_layers["Layer0"] = + std::make_unique>(1., + 2.); + material_layers["Layer1"] = + std::make_unique>(3., + 4.); + const typename OptionTags::ConstitutiveRelationPerBlock::type + material_layers_variant(std::move(material_layers)); + + auto material_per_block = + Tags::ConstitutiveRelationPerBlock::create_from_options( + domain_creator, material_layers_variant); + auto material_block_groups = + Tags::MaterialBlockGroups::create_from_options( + material_layers_variant); + REQUIRE(material_block_groups == + std::unordered_set{"Layer0", "Layer1"}); + + using element_array = ElementArray; + ActionTesting::MockRuntimeSystem runner{ + tuples::TaggedTuple, + Tags::ConstitutiveRelationPerBlock, + Tags::MaterialBlockGroups>{ + domain_creator->create_domain(), std::move(material_per_block), + std::move(material_block_groups)}}; + for (const auto& element_id : {element_id_layer1, element_id_layer2}) { + ActionTesting::emplace_component_and_initialize( + &runner, element_id, {Element{element_id, {}}}); + } + ActionTesting::set_phase(make_not_null(&runner), Parallel::Phase::Testing); + for (const auto& element_id : {element_id_layer1, element_id_layer2}) { + ActionTesting::next_action(make_not_null(&runner), + element_id); + } + + const auto get_tag = + [&runner](auto tag_v, const ElementId& element_id) -> const auto& { + using tag = std::decay_t; + return ActionTesting::get_databox_tag(runner, + element_id); + }; + + const auto check_material = [&get_tag](const ElementId& element_id, + const double expected_bulk_modulus) { + const auto& material = + dynamic_cast&>( + get_tag(Tags::ConstitutiveRelation{}, element_id)); + CHECK(material.bulk_modulus() == expected_bulk_modulus); + }; + + check_material(element_id_layer1, 1.); + check_material(element_id_layer2, 3.); +} +} // namespace Elasticity diff --git a/tests/Unit/Elliptic/Systems/Elasticity/CMakeLists.txt b/tests/Unit/Elliptic/Systems/Elasticity/CMakeLists.txt index 1ce9801a9b86..8bbf7311186d 100644 --- a/tests/Unit/Elliptic/Systems/Elasticity/CMakeLists.txt +++ b/tests/Unit/Elliptic/Systems/Elasticity/CMakeLists.txt @@ -22,4 +22,5 @@ target_link_libraries( Utilities ) +add_subdirectory(Actions) add_subdirectory(BoundaryConditions) diff --git a/tests/Unit/Elliptic/Systems/Elasticity/Test_Equations.cpp b/tests/Unit/Elliptic/Systems/Elasticity/Test_Equations.cpp index db92f3ef5e09..7795357710d2 100644 --- a/tests/Unit/Elliptic/Systems/Elasticity/Test_Equations.cpp +++ b/tests/Unit/Elliptic/Systems/Elasticity/Test_Equations.cpp @@ -13,6 +13,7 @@ #include "DataStructures/Tensor/Tensor.hpp" #include "Domain/Tags.hpp" #include "Elliptic/Protocols/FirstOrderSystem.hpp" +#include "Elliptic/Systems/Elasticity/Actions/InitializeConstitutiveRelation.hpp" #include "Elliptic/Systems/Elasticity/Equations.hpp" #include "Elliptic/Systems/Elasticity/FirstOrderSystem.hpp" #include "Elliptic/Systems/Elasticity/Tags.hpp" @@ -64,24 +65,26 @@ void test_computers(const DataVector& used_for_size) { using deriv_field_tag = ::Tags::deriv, tmpl::size_t, Frame::Inertial>; using field_flux_tag = Elasticity::Tags::MinusStress; - using constitutive_relation_tag = Elasticity::Tags::ConstitutiveRelation; + using constitutive_relation_tag = + Elasticity::Tags::ConstitutiveRelationPerBlock; using coordinates_tag = domain::Tags::Coordinates; const size_t num_points = used_for_size.size(); { INFO("Fluxes" << Dim << "D"); - std::unique_ptr< - Elasticity::ConstitutiveRelations::ConstitutiveRelation> - constitutive_relation = std::make_unique< - Elasticity::ConstitutiveRelations::IsotropicHomogeneous>(1., - 2.); - auto box = db::create< - db::AddSimpleTags>( + const Elasticity::ConstitutiveRelations::IsotropicHomogeneous + constitutive_relation{1., 2.}; + std::vector>> + constitutive_relations{}; + constitutive_relations.push_back(constitutive_relation.get_clone()); + auto box = db::create, coordinates_tag>>( tnsr::I{num_points, 1.}, tnsr::iJ{num_points, 2.}, tnsr::II{num_points, std::numeric_limits::signaling_NaN()}, - std::move(constitutive_relation), + std::move(constitutive_relations), Element{ElementId{0}, {}}, tnsr::I{num_points, 6.}); const Elasticity::Fluxes fluxes_computer{}; @@ -92,9 +95,9 @@ void test_computers(const DataVector& used_for_size) { get(box)); auto expected_field_flux = tnsr::II{ num_points, std::numeric_limits::signaling_NaN()}; - Elasticity::primal_fluxes( - make_not_null(&expected_field_flux), get(box), - get(box), get(box)); + Elasticity::primal_fluxes(make_not_null(&expected_field_flux), + get(box), constitutive_relation, + get(box)); CHECK(get(box) == expected_field_flux); } } diff --git a/tests/Unit/NumericalAlgorithms/LinearOperators/CMakeLists.txt b/tests/Unit/NumericalAlgorithms/LinearOperators/CMakeLists.txt index df42789485e8..e0a90e16a134 100644 --- a/tests/Unit/NumericalAlgorithms/LinearOperators/CMakeLists.txt +++ b/tests/Unit/NumericalAlgorithms/LinearOperators/CMakeLists.txt @@ -24,6 +24,7 @@ target_link_libraries( DataStructures DiscontinuousGalerkin Domain + DomainCreators ErrorHandling LinearOperators MathFunctions diff --git a/tests/Unit/ParallelAlgorithms/Events/Test_ObserveNorms.cpp b/tests/Unit/ParallelAlgorithms/Events/Test_ObserveNorms.cpp index ee7f30885444..6b1a4c5de1ff 100644 --- a/tests/Unit/ParallelAlgorithms/Events/Test_ObserveNorms.cpp +++ b/tests/Unit/ParallelAlgorithms/Events/Test_ObserveNorms.cpp @@ -81,6 +81,9 @@ struct Var0TimesThreeCompute : db::ComputeTag, struct TestSectionIdTag {}; +// Name for option parsing +struct ObserveMyNorms {}; + struct MockContributeReductionData { struct Results { observers::ObservationId observation_id; @@ -151,12 +154,12 @@ struct MockObserverComponent { Parallel::PhaseActions>>; }; -template +template using ObserveNormsEvent = Events::ObserveNorms< tmpl::list, - tmpl::list, ArraySectionIdTag>; + tmpl::list, ArraySectionIdTag, OptionName>; -template +template struct Metavariables { static constexpr size_t volume_dim = Dim; using component_list = tmpl::list, @@ -164,8 +167,8 @@ struct Metavariables { struct factory_creation : tt::ConformsTo { - using factory_classes = tmpl::map< - tmpl::pair>>>; + using factory_classes = tmpl::map>>>; }; }; @@ -300,7 +303,6 @@ void test(const std::unique_ptr observe, CHECK(results.volume_integral_values[1] == approx(41.0)); CHECK(results.volume_integral_values[2] == approx(68.0)); CHECK(results.volume_integral_values[3] == approx(95.0)); - } } // namespace @@ -374,6 +376,21 @@ SPECTRE_TEST_CASE("Unit.Evolution.ObserveNorms", "[Unit][Evolution]") { test(std::move(serialized_event), Spectral::Basis::Legendre, Spectral::Quadrature::GaussLobatto, std::nullopt); + { + INFO("Test option name"); + // Test option name + TestHelpers::test_creation, + Metavariables<3, void, ObserveMyNorms>>( + R"( + ObserveMyNorms: + SubfileName: reduction0 + TensorsToObserve: + - Name: Var0 + NormType: Max + Components: Individual + )"); + } + test(std::make_unique>(ObserveNormsEvent{ "reduction0", {{"Var0", "Max", "Individual"}, diff --git a/tests/Unit/PointwiseFunctions/AnalyticSolutions/Elasticity/Test_BentBeam.cpp b/tests/Unit/PointwiseFunctions/AnalyticSolutions/Elasticity/Test_BentBeam.cpp index dfa4c76a4279..6f5f87002481 100644 --- a/tests/Unit/PointwiseFunctions/AnalyticSolutions/Elasticity/Test_BentBeam.cpp +++ b/tests/Unit/PointwiseFunctions/AnalyticSolutions/Elasticity/Test_BentBeam.cpp @@ -152,10 +152,16 @@ SPECTRE_TEST_CASE( INFO("Test elasticity system with bent beam"); using system = Elasticity::FirstOrderSystem<2>; // Verify that the solution numerically solves the system + std::vector>> + constitutive_relations{}; + constitutive_relations.push_back(constitutive_relation.get_clone()); + const Element<2> element{ElementId<2>{0}, {}}; FirstOrderEllipticSolutionsTestHelpers::verify_solution( solution, mesh, coord_map, std::numeric_limits::epsilon() * 100., - std::make_tuple(constitutive_relation, inertial_coords)); + std::make_tuple(std::move(constitutive_relations), element, + inertial_coords)); }; { diff --git a/tests/Unit/PointwiseFunctions/AnalyticSolutions/Elasticity/Test_HalfSpaceMirror.cpp b/tests/Unit/PointwiseFunctions/AnalyticSolutions/Elasticity/Test_HalfSpaceMirror.cpp index bdb1f90df516..0e53e4a41f54 100644 --- a/tests/Unit/PointwiseFunctions/AnalyticSolutions/Elasticity/Test_HalfSpaceMirror.cpp +++ b/tests/Unit/PointwiseFunctions/AnalyticSolutions/Elasticity/Test_HalfSpaceMirror.cpp @@ -145,9 +145,15 @@ SPECTRE_TEST_CASE( Spectral::Quadrature::GaussLobatto}; const auto logical_coords = logical_coordinates(mesh); const auto inertial_coords = coord_map(logical_coords); + std::vector>> + constitutive_relations{}; + constitutive_relations.push_back(constitutive_relation.get_clone()); + const Element<3> element{ElementId<3>{0}, {}}; FirstOrderEllipticSolutionsTestHelpers::verify_solution( solution, mesh, coord_map, 0.05, - std::make_tuple(constitutive_relation, inertial_coords)); + std::make_tuple(std::move(constitutive_relations), element, + inertial_coords)); }; CHECK_THROWS_WITH( diff --git a/tests/Unit/PointwiseFunctions/AnalyticSolutions/Elasticity/Test_Zero.cpp b/tests/Unit/PointwiseFunctions/AnalyticSolutions/Elasticity/Test_Zero.cpp index cfdf9007de2d..c3c7f83c3989 100644 --- a/tests/Unit/PointwiseFunctions/AnalyticSolutions/Elasticity/Test_Zero.cpp +++ b/tests/Unit/PointwiseFunctions/AnalyticSolutions/Elasticity/Test_Zero.cpp @@ -59,14 +59,19 @@ void test_solution() { test_copy_semantics(solution); using system = Elasticity::FirstOrderSystem; - Elasticity::ConstitutiveRelations::IsotropicHomogeneous - constitutive_relation{1., 1.}; + auto constitutive_relation = std::make_unique< + Elasticity::ConstitutiveRelations::IsotropicHomogeneous>(1., 1.); + std::vector>> + constitutive_relations{}; + constitutive_relations.push_back(std::move(constitutive_relation)); const Mesh mesh{12, Spectral::Basis::Legendre, Spectral::Quadrature::GaussLobatto}; + const Element element{ElementId{0}, {}}; const auto coord_map = make_coord_map(); FirstOrderEllipticSolutionsTestHelpers::verify_solution( solution, mesh, coord_map, 1.e-14, - std::make_tuple(constitutive_relation, + std::make_tuple(std::move(constitutive_relations), element, coord_map(logical_coordinates(mesh)))); } diff --git a/tests/Unit/PointwiseFunctions/Elasticity/ConstitutiveRelations/Test_CubicCrystal.cpp b/tests/Unit/PointwiseFunctions/Elasticity/ConstitutiveRelations/Test_CubicCrystal.cpp index bfcf7331ccd2..be62a52f5708 100644 --- a/tests/Unit/PointwiseFunctions/Elasticity/ConstitutiveRelations/Test_CubicCrystal.cpp +++ b/tests/Unit/PointwiseFunctions/Elasticity/ConstitutiveRelations/Test_CubicCrystal.cpp @@ -19,28 +19,30 @@ #include "PointwiseFunctions/Elasticity/ConstitutiveRelations/IsotropicHomogeneous.hpp" #include "Utilities/Gsl.hpp" #include "Utilities/MakeWithValue.hpp" -// IWYU pragma: no_forward_declare Tensor +namespace Elasticity::ConstitutiveRelations { namespace { void test_semantics() { INFO("Test type traits"); - const Elasticity::ConstitutiveRelations::CubicCrystal relation{3., 2., 1.}; - CHECK(relation == - Elasticity::ConstitutiveRelations::CubicCrystal{3., 2., 1.}); - CHECK(relation != - Elasticity::ConstitutiveRelations::CubicCrystal{2., 2., 2.}); - CHECK(relation != - Elasticity::ConstitutiveRelations::CubicCrystal{3., 0.2, 1.}); + const CubicCrystal relation{3., 2., 1.}; + CHECK(relation == CubicCrystal{3., 2., 1.}); + CHECK(relation != CubicCrystal{2., 2., 2.}); + CHECK(relation != CubicCrystal{3., 0.2, 1.}); test_serialization(relation); test_copy_semantics(relation); - const auto created_relation = TestHelpers::test_creation< - Elasticity::ConstitutiveRelations::CubicCrystal>( - "C_11: 3.\n" - "C_12: 2.\n" - "C_44: 1.\n"); + const auto created_ptr = + TestHelpers::test_factory_creation, CubicCrystal>( + "CubicCrystal:\n" + " C_11: 3.\n" + " C_12: 2.\n" + " C_44: 1.\n") + ->get_clone(); + REQUIRE(dynamic_cast(created_ptr.get()) != nullptr); + const auto& created_relation = + dynamic_cast(*created_ptr); CHECK(created_relation == relation); - Elasticity::ConstitutiveRelations::CubicCrystal moved_relation{3., 2., 1.}; + CubicCrystal moved_relation{3., 2., 1.}; test_move_semantics(std::move(moved_relation), relation); } @@ -53,17 +55,15 @@ void test_consistency(const tnsr::ii& random_strain, youngs_modulus / (1. - 2. * poisson_ratio) / 3.; const double isotropic_shear_modulus = youngs_modulus / (1. + poisson_ratio) / 2.; - const Elasticity::ConstitutiveRelations::IsotropicHomogeneous<3> - isotropic_homogeneous_relation{isotropic_bulk_modulus, - isotropic_shear_modulus}; - const double c_11 = youngs_modulus * (1.- poisson_ratio) / + const IsotropicHomogeneous<3> isotropic_homogeneous_relation{ + isotropic_bulk_modulus, isotropic_shear_modulus}; + const double c_11 = youngs_modulus * (1. - poisson_ratio) / ((1. + poisson_ratio) * (1. - 2. * poisson_ratio)); const double c_12 = youngs_modulus * poisson_ratio / ((1. + poisson_ratio) * (1. - 2. * poisson_ratio)); const double c_44 = isotropic_shear_modulus; // this should be isotropic homogeneous - const Elasticity::ConstitutiveRelations::CubicCrystal - cubic_crystalline_relation{c_11, c_12, c_44}; + const CubicCrystal cubic_crystalline_relation{c_11, c_12, c_44}; tnsr::II cubic_crystalline_stress{ random_strain.begin()->size()}; cubic_crystalline_relation.stress(make_not_null(&cubic_crystalline_stress), @@ -80,7 +80,7 @@ void test_consistency(const tnsr::ii& random_strain, } } // This relation should be the negative identity - const Elasticity::ConstitutiveRelations::CubicCrystal relation{1., 0., 0.5}; + const CubicCrystal relation{1., 0., 0.5}; relation.stress(make_not_null(&cubic_crystalline_stress), random_strain, random_inertial_coords); for (size_t i = 0; i < 3; i++) { @@ -93,13 +93,12 @@ void test_consistency(const tnsr::ii& random_strain, void test_implementation(const double youngs_modulus, const double poisson_ratio, double shear_modulus) { - const Elasticity::ConstitutiveRelations::CubicCrystal relation{ - youngs_modulus, poisson_ratio, shear_modulus}; + const CubicCrystal relation{youngs_modulus, poisson_ratio, shear_modulus}; pypp::check_with_random_values<1>( - static_cast*>, const tnsr::ii&, const tnsr::I&) const>( - &Elasticity::ConstitutiveRelations::CubicCrystal::stress), + &CubicCrystal::stress), relation, "CubicCrystal", {"stress"}, {{{-1., 1.}}}, std::tuple{youngs_modulus, poisson_ratio, shear_modulus}, @@ -161,3 +160,5 @@ SPECTRE_TEST_CASE("Unit.Elasticity.ConstitutiveRelations.CubicCrystal", test_implementation_suite(); } } + +} // namespace Elasticity::ConstitutiveRelations diff --git a/tests/Unit/PointwiseFunctions/Elasticity/ConstitutiveRelations/Test_IsotropicHomogeneous.cpp b/tests/Unit/PointwiseFunctions/Elasticity/ConstitutiveRelations/Test_IsotropicHomogeneous.cpp index b147e52a6772..bc6d00672e3c 100644 --- a/tests/Unit/PointwiseFunctions/Elasticity/ConstitutiveRelations/Test_IsotropicHomogeneous.cpp +++ b/tests/Unit/PointwiseFunctions/Elasticity/ConstitutiveRelations/Test_IsotropicHomogeneous.cpp @@ -18,46 +18,44 @@ #include "PointwiseFunctions/Elasticity/ConstitutiveRelations/IsotropicHomogeneous.hpp" #include "Utilities/Gsl.hpp" #include "Utilities/MakeWithValue.hpp" -// IWYU pragma: no_forward_declare Tensor +namespace Elasticity::ConstitutiveRelations { namespace { template void test_type_traits() { INFO("Test type traits"); - const Elasticity::ConstitutiveRelations::IsotropicHomogeneous relation{ - 1., 2.}; - CHECK(relation == - Elasticity::ConstitutiveRelations::IsotropicHomogeneous{1., 2.}); - CHECK(relation != - Elasticity::ConstitutiveRelations::IsotropicHomogeneous{1., 1.}); - CHECK(relation != - Elasticity::ConstitutiveRelations::IsotropicHomogeneous{2., 2.}); + const IsotropicHomogeneous relation{1., 2.}; + CHECK(relation == IsotropicHomogeneous{1., 2.}); + CHECK(relation != IsotropicHomogeneous{1., 1.}); + CHECK(relation != IsotropicHomogeneous{2., 2.}); test_serialization(relation); test_copy_semantics(relation); - const auto created_relation = TestHelpers::test_creation< - Elasticity::ConstitutiveRelations::IsotropicHomogeneous>( - "BulkModulus: 1.\n" - "ShearModulus: 2.\n"); + const auto created_ptr = + TestHelpers::test_factory_creation, + IsotropicHomogeneous>( + "IsotropicHomogeneous:\n" + " BulkModulus: 1.\n" + " ShearModulus: 2.\n") + ->get_clone(); + REQUIRE(dynamic_cast*>(created_ptr.get()) != + nullptr); + const auto& created_relation = + dynamic_cast&>(*created_ptr); CHECK(created_relation == relation); - Elasticity::ConstitutiveRelations::IsotropicHomogeneous moved_relation{ - 1., 2.}; + IsotropicHomogeneous moved_relation{1., 2.}; test_move_semantics(std::move(moved_relation), relation); } template void test_implementation(const double incompressibility, const double rigidity) { - const Elasticity::ConstitutiveRelations::IsotropicHomogeneous relation{ - incompressibility, rigidity}; + const IsotropicHomogeneous relation{incompressibility, rigidity}; pypp::check_with_random_values<1>( - static_cast::*)( + static_cast::*)( gsl::not_null*>, const tnsr::ii&, const tnsr::I&) - const>( - &Elasticity::ConstitutiveRelations::IsotropicHomogeneous< - Dim>::stress), + const>(&IsotropicHomogeneous::stress), relation, "IsotropicHomogeneous", {"stress"}, {{{-10.0, 10.0}}}, std::tuple{incompressibility, rigidity}, DataVector(5)); } @@ -87,8 +85,7 @@ void test_identity(const tnsr::ii& random_strain, const tnsr::I& random_inertial_coords) { INFO("Identity"); // This relation should be the negative identity - const Elasticity::ConstitutiveRelations::IsotropicHomogeneous relation{ - 1. / 3., 1. / 2.}; + const IsotropicHomogeneous relation{1. / 3., 1. / 2.}; tnsr::II stress{random_strain.begin()->size()}; relation.stress(make_not_null(&stress), random_strain, random_inertial_coords); @@ -109,8 +106,7 @@ void test_trace<3>(const tnsr::ii& random_strain, INFO("Trace"); // This relation should result in a stress trace that is equal the negative // trace of the strain. The shear modulus should be irrelevant here. - const Elasticity::ConstitutiveRelations::IsotropicHomogeneous<3> relation{ - 1. / 3., 10.}; + const IsotropicHomogeneous<3> relation{1. / 3., 10.}; tnsr::II stress{random_strain.begin()->size()}; relation.stress(make_not_null(&stress), random_strain, random_inertial_coords); @@ -127,8 +123,7 @@ template <> void test_trace<2>(const tnsr::ii& random_strain, const tnsr::I& random_inertial_coords) { INFO("Trace"); - const Elasticity::ConstitutiveRelations::IsotropicHomogeneous<2> relation{ - 1. / 3., 2.}; + const IsotropicHomogeneous<2> relation{1. / 3., 2.}; tnsr::II stress{random_strain.begin()->size()}; relation.stress(make_not_null(&stress), random_strain, random_inertial_coords); @@ -148,8 +143,7 @@ void test_traceless(const tnsr::ii& random_strain, INFO("Traceless"); // This relation should result in a traceless stress. // The shear modulus should be irrelevant here. - const Elasticity::ConstitutiveRelations::IsotropicHomogeneous relation{ - 0., 10.}; + const IsotropicHomogeneous relation{0., 10.}; tnsr::II stress{random_strain.begin()->size()}; relation.stress(make_not_null(&stress), random_strain, random_inertial_coords); @@ -201,3 +195,5 @@ SPECTRE_TEST_CASE("Unit.Elasticity.ConstitutiveRelations.IsotropicHomogeneous", test_analytically<2>(); } } + +} // namespace Elasticity::ConstitutiveRelations