diff --git a/src/DataStructures/CMakeLists.txt b/src/DataStructures/CMakeLists.txt index 1396ffa167a7..4cd42f9e3af9 100644 --- a/src/DataStructures/CMakeLists.txt +++ b/src/DataStructures/CMakeLists.txt @@ -41,6 +41,7 @@ spectre_target_headers( DynamicBuffer.hpp DynamicMatrix.hpp DynamicVector.hpp + ExtractPoint.hpp FixedHashMap.hpp FloatingPointType.hpp IdPair.hpp diff --git a/src/DataStructures/ExtractPoint.hpp b/src/DataStructures/ExtractPoint.hpp new file mode 100644 index 000000000000..1f004c5273d3 --- /dev/null +++ b/src/DataStructures/ExtractPoint.hpp @@ -0,0 +1,112 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include + +#include "DataStructures/DataVector.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "DataStructures/Variables.hpp" +#include "Utilities/ErrorHandling/Assert.hpp" +#include "Utilities/Gsl.hpp" +#include "Utilities/TMPL.hpp" + +/// \ingroup DataStructuresGroup +/// Copy a given index of each component of a `Tensor` or +/// `Variables` into a `Tensor`, single point +/// `Tensor`, or single-point `Variables`. +/// +/// \note There is no by-value overload extracting to a +/// `Tensor`. This is both for the practical reason that +/// it would be ambiguous with the `Tensor` overload and +/// because allocating multiple `DataVector`s for the return type +/// would usually be very inefficient. +/// +/// \see overwrite_point +/// @{ +template +void extract_point( + const gsl::not_null*> destination, + const Tensor& source, const size_t index) { + for (size_t i = 0; i < destination->size(); ++i) { + (*destination)[i] = source[i][index]; + } +} + +template +Tensor extract_point( + const Tensor& tensor, const size_t index) { + Tensor result; + extract_point(make_not_null(&result), tensor, index); + return result; +} + +template +void extract_point( + const gsl::not_null*> destination, + const Tensor& source, const size_t index) { + ASSERT(destination->begin()->size() == 1, + "Output tensor components have wrong size: " + << destination->begin()->size()); + for (size_t i = 0; i < destination->size(); ++i) { + (*destination)[i][0] = source[i][index]; + } +} + +template +void extract_point( + const gsl::not_null>*> result, + const Variables>& variables, const size_t index) { + result->initialize(1); + expand_pack((extract_point( + make_not_null(&get(*result)), get(variables), index), 0)...); +} + +template +Variables> extract_point( + const Variables>& variables, const size_t index) { + Variables> result(1); + extract_point(make_not_null(&result), variables, index); + return result; +} +/// @} + +/// \ingroup DataStructuresGroup +/// Copy a `Tensor`, single point `Tensor`, or +/// single-point `Variables` into the given index of each +/// component of a `Tensor` or `Variables`. +/// +/// \see extract_point +/// @{ +template +void overwrite_point( + const gsl::not_null*> destination, + const Tensor& source, const size_t index) { + for (size_t i = 0; i < destination->size(); ++i) { + (*destination)[i][index] = source[i]; + } +} + +template +void overwrite_point( + const gsl::not_null*> destination, + const Tensor& source, const size_t index) { + ASSERT(source.begin()->size() == 1, + "Cannot overwrite with " << source.begin()->size() << " points."); + for (size_t i = 0; i < destination->size(); ++i) { + (*destination)[i][index] = source[i][0]; + } +} + +template +void overwrite_point( + const gsl::not_null>*> destination, + const Variables>& source, const size_t index) { + ASSERT(source.number_of_grid_points() == 1, + "Must overwrite with a single point."); + expand_pack((overwrite_point(make_not_null(&get(*destination)), + extract_point(get(source), 0), index), + 0)...); +} +/// @} diff --git a/src/DataStructures/Tensor/CMakeLists.txt b/src/DataStructures/Tensor/CMakeLists.txt index a549370ebb54..a80b062d22cc 100644 --- a/src/DataStructures/Tensor/CMakeLists.txt +++ b/src/DataStructures/Tensor/CMakeLists.txt @@ -5,7 +5,6 @@ spectre_target_headers( ${LIBRARY} INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src HEADERS - ExtractPoint.hpp Identity.hpp IndexType.hpp Metafunctions.hpp diff --git a/src/DataStructures/Tensor/ExtractPoint.hpp b/src/DataStructures/Tensor/ExtractPoint.hpp deleted file mode 100644 index 470786c1e8b2..000000000000 --- a/src/DataStructures/Tensor/ExtractPoint.hpp +++ /dev/null @@ -1,22 +0,0 @@ -// Distributed under the MIT License. -// See LICENSE.txt for details. - -#pragma once - -#include - -#include "DataStructures/DataVector.hpp" -#include "DataStructures/Tensor/Tensor.hpp" - -/// \ingroup DataStructuresGroup -/// Copy a given index of each component of a `Tensor` -/// into a `Tensor`. -template -Tensor extract_point( - const Tensor& tensor, const size_t index) { - Tensor result; - for (size_t i = 0; i < result.size(); ++i) { - result[i] = tensor[i][index]; - } - return result; -} diff --git a/src/Evolution/DgSubcell/SetInterpolators.hpp b/src/Evolution/DgSubcell/SetInterpolators.hpp index 6125c6947010..fd514325633d 100644 --- a/src/Evolution/DgSubcell/SetInterpolators.hpp +++ b/src/Evolution/DgSubcell/SetInterpolators.hpp @@ -8,8 +8,8 @@ #include #include +#include "DataStructures/ExtractPoint.hpp" #include "DataStructures/FixedHashMap.hpp" -#include "DataStructures/Tensor/ExtractPoint.hpp" #include "DataStructures/Tensor/Tensor.hpp" #include "Domain/Creators/Tags/Domain.hpp" #include "Domain/Domain.hpp" diff --git a/src/Evolution/Systems/GrMhd/ValenciaDivClean/FixConservatives.cpp b/src/Evolution/Systems/GrMhd/ValenciaDivClean/FixConservatives.cpp index 1d23f71360ce..d59364e9e3b7 100644 --- a/src/Evolution/Systems/GrMhd/ValenciaDivClean/FixConservatives.cpp +++ b/src/Evolution/Systems/GrMhd/ValenciaDivClean/FixConservatives.cpp @@ -10,9 +10,9 @@ #include // IWYU pragma: keep #include "DataStructures/DataVector.hpp" +#include "DataStructures/ExtractPoint.hpp" #include "DataStructures/Tags/TempTensor.hpp" #include "DataStructures/Tensor/EagerMath/DotProduct.hpp" -#include "DataStructures/Tensor/ExtractPoint.hpp" #include "DataStructures/Tensor/Tensor.hpp" #include "DataStructures/Variables.hpp" #include "NumericalAlgorithms/RootFinding/TOMS748.hpp" diff --git a/tests/Unit/DataStructures/CMakeLists.txt b/tests/Unit/DataStructures/CMakeLists.txt index 695767b6e2dd..5d387be44b8e 100644 --- a/tests/Unit/DataStructures/CMakeLists.txt +++ b/tests/Unit/DataStructures/CMakeLists.txt @@ -16,6 +16,7 @@ set(LIBRARY_SOURCES Test_DataVectorBinaryOperations.cpp Test_DiagonalModalOperator.cpp Test_DynamicBuffer.cpp + Test_ExtractPoint.cpp Test_FixedHashMap.cpp Test_FloatingPointType.cpp Test_IdPair.cpp diff --git a/tests/Unit/DataStructures/Tensor/CMakeLists.txt b/tests/Unit/DataStructures/Tensor/CMakeLists.txt index 4fc0e087a5ab..6e152a01037e 100644 --- a/tests/Unit/DataStructures/Tensor/CMakeLists.txt +++ b/tests/Unit/DataStructures/Tensor/CMakeLists.txt @@ -4,7 +4,6 @@ set(LIBRARY "Test_Tensor") set(LIBRARY_SOURCES - Test_ExtractPoint.cpp Test_Identity.cpp Test_Metafunctions.cpp Test_Slice.cpp diff --git a/tests/Unit/DataStructures/Tensor/Test_ExtractPoint.cpp b/tests/Unit/DataStructures/Tensor/Test_ExtractPoint.cpp deleted file mode 100644 index daaf189bf37d..000000000000 --- a/tests/Unit/DataStructures/Tensor/Test_ExtractPoint.cpp +++ /dev/null @@ -1,35 +0,0 @@ -// Distributed under the MIT License. -// See LICENSE.txt for details. - -#include "Framework/TestingFramework.hpp" - -#include "DataStructures/DataVector.hpp" -#include "DataStructures/Tensor/ExtractPoint.hpp" -#include "DataStructures/Tensor/Tensor.hpp" - -SPECTRE_TEST_CASE("Unit.DataStructures.Tensor.ExtaractPoint", - "[DataStructures][Unit]") { - const Scalar scalar{DataVector{3.0, 4.0, 5.0}}; - CHECK(extract_point(scalar, 0) == Scalar{3.0}); - CHECK(extract_point(scalar, 1) == Scalar{4.0}); - CHECK(extract_point(scalar, 2) == Scalar{5.0}); - - tnsr::ii tensor; - get<0, 0>(tensor) = DataVector{1.0, 2.0}; - get<0, 1>(tensor) = DataVector{3.0, 4.0}; - get<1, 1>(tensor) = DataVector{5.0, 6.0}; - { - tnsr::ii expected; - get<0, 0>(expected) = 1.0; - get<0, 1>(expected) = 3.0; - get<1, 1>(expected) = 5.0; - CHECK(extract_point(tensor, 0) == expected); - } - { - tnsr::ii expected; - get<0, 0>(expected) = 2.0; - get<0, 1>(expected) = 4.0; - get<1, 1>(expected) = 6.0; - CHECK(extract_point(tensor, 1) == expected); - } -} diff --git a/tests/Unit/DataStructures/Test_ExtractPoint.cpp b/tests/Unit/DataStructures/Test_ExtractPoint.cpp new file mode 100644 index 000000000000..c2f6b9545cd1 --- /dev/null +++ b/tests/Unit/DataStructures/Test_ExtractPoint.cpp @@ -0,0 +1,118 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Framework/TestingFramework.hpp" + +#include + +#include "DataStructures/DataBox/Tag.hpp" +#include "DataStructures/DataVector.hpp" +#include "DataStructures/ExtractPoint.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "Utilities/Literals.hpp" +#include "Utilities/TMPL.hpp" + +namespace { +struct ScalarTag : db::SimpleTag { + using type = Scalar; +}; + +struct TensorTag : db::SimpleTag { + using type = tnsr::ii; +}; + +template +Tensor promote_to_data_vectors( + const Tensor& tensor) { + Tensor result(1_st); + for (size_t i = 0; i < result.size(); ++i) { + result[i] = tensor[i]; + } + return result; +} +} // namespace + +SPECTRE_TEST_CASE("Unit.DataStructures.ExtractPoint", + "[DataStructures][Unit]") { + // Test with owning Tensors as well as Variables to make sure there + // are no assumptions about memory layout in the tensor versions. + const Scalar scalar{DataVector{3.0, 4.0}}; + tnsr::ii tensor; + get<0, 0>(tensor) = DataVector{1.0, 2.0}; + get<0, 1>(tensor) = DataVector{3.0, 4.0}; + get<1, 1>(tensor) = DataVector{5.0, 6.0}; + Variables> variables(2); + get(variables) = scalar; + get(variables) = tensor; + + Scalar reconstructed_scalar{DataVector(2)}; + tnsr::ii reconstructed_tensor(DataVector(2)); + Scalar reconstructed_scalar_from_dv{DataVector(2)}; + tnsr::ii reconstructed_tensor_from_dv(DataVector(2)); + Variables> reconstructed_variables(2); + + const auto check_point = [&](const size_t index, + const double scalar_component, + const std::array& tensor_components) { + const Scalar expected_scalar{scalar_component}; + tnsr::ii expected_tensor; + get<0, 0>(expected_tensor) = tensor_components[0]; + get<0, 1>(expected_tensor) = tensor_components[1]; + get<1, 1>(expected_tensor) = tensor_components[2]; + const Scalar expected_scalar_dv = + promote_to_data_vectors(expected_scalar); + const tnsr::ii expected_tensor_dv = + promote_to_data_vectors(expected_tensor); + Variables> expected_variables(1); + get(expected_variables) = expected_scalar_dv; + get(expected_variables) = expected_tensor_dv; + + { + Scalar result{}; + extract_point(make_not_null(&result), scalar, index); + CHECK(result == expected_scalar); + } + { + tnsr::ii result{}; + extract_point(make_not_null(&result), tensor, index); + CHECK(result == expected_tensor); + } + CHECK(extract_point(scalar, index) == expected_scalar); + CHECK(extract_point(tensor, index) == expected_tensor); + { + Scalar result(1_st); + extract_point(make_not_null(&result), scalar, index); + CHECK(result == expected_scalar_dv); + } + { + tnsr::ii result(1_st); + extract_point(make_not_null(&result), tensor, index); + CHECK(result == expected_tensor_dv); + } + { + Variables> result(1); + extract_point(make_not_null(&result), variables, index); + CHECK(result == expected_variables); + } + CHECK(extract_point(variables, index) == expected_variables); + + overwrite_point(make_not_null(&reconstructed_scalar), expected_scalar, + index); + overwrite_point(make_not_null(&reconstructed_tensor), expected_tensor, + index); + overwrite_point(make_not_null(&reconstructed_scalar_from_dv), + expected_scalar_dv, index); + overwrite_point(make_not_null(&reconstructed_tensor_from_dv), + expected_tensor_dv, index); + overwrite_point(make_not_null(&reconstructed_variables), expected_variables, + index); + }; + check_point(0, 3.0, {{1.0, 3.0, 5.0}}); + check_point(1, 4.0, {{2.0, 4.0, 6.0}}); + + CHECK(reconstructed_scalar == scalar); + CHECK(reconstructed_tensor == tensor); + CHECK(reconstructed_scalar_from_dv == scalar); + CHECK(reconstructed_tensor_from_dv == tensor); + CHECK(reconstructed_variables == variables); +}