diff --git a/DiFfRG/include/DiFfRG/discretization/data/csv_output.hh b/DiFfRG/include/DiFfRG/discretization/data/csv_output.hh index 1d146e365..3db0c5093 100644 --- a/DiFfRG/include/DiFfRG/discretization/data/csv_output.hh +++ b/DiFfRG/include/DiFfRG/discretization/data/csv_output.hh @@ -52,7 +52,7 @@ namespace DiFfRG void set_Lambda(const double Lambda); private: - const JSONValue &json; + [[maybe_unused]] const JSONValue &json; const std::string top_folder; const std::string output_name; std::ofstream output_stream; diff --git a/DiFfRG/include/DiFfRG/discretization/mesh/configuration_mesh.hh b/DiFfRG/include/DiFfRG/discretization/mesh/configuration_mesh.hh new file mode 100644 index 000000000..a4a47798d --- /dev/null +++ b/DiFfRG/include/DiFfRG/discretization/mesh/configuration_mesh.hh @@ -0,0 +1,189 @@ +#pragma once +#include "DiFfRG/common/math.hh" +#include +#include +#include +#include + +namespace DiFfRG +{ + + namespace Config + { + + namespace internal + { + inline std::vector string_to_substrings_array(const std::string str) + { + std::vector array; + std::istringstream ss(str); + std::string buf; + while (std::getline(ss, buf, ',')) + array.push_back(buf); + return array; + } + inline std::array string_to_range(const std::string str) + { + std::vector vec; + std::istringstream ss(str); + std::string buf; + while (std::getline(ss, buf, ':')) + vec.push_back(std::stod(buf)); + if (vec.size() != 3) throw std::runtime_error(str + " is not a range"); + + std::array arr; + std::copy_n(vec.begin(), 3, arr.begin()); + return arr; + } + } // namespace internal + + class GridAxis + { + public: + GridAxis() = delete; + GridAxis(std::string json_string) + { + auto subrange = internal::string_to_range(json_string); + min = subrange[0]; + step = subrange[1]; + max = subrange[2]; + + validate(); + }; + + static std::string get_default_range() { return "0:1e-4:7e-3,7e-3:5e-4:9e-3,9e-3:1e-3:1.1e-2"; } + + double min; + double max; + double step; + + /** + * Return a vector of uniform step widths covering [min, max]. + * If the configured step fits the interval (within is_close) it is used; + * otherwise the count is increased by one and the width is adjusted. + */ + std::vector get_stepwiths() const + { + double local_step = step; + auto steps = static_cast((max - min) / local_step); + if (!is_close((min + local_step * steps), max)) { + local_step = (max - min) / (steps + 1); + ++steps; + } + return std::vector(steps, local_step); + } + + private: + void validate() + { + bool condition = (min < max) && (step <= (max - min)); + if (!condition) { + throw std::runtime_error(std::format( + "invalid range: min={}, step={}, max={} (require min < max and 0 < step <= max-min)", min, step, max)); + } + } + }; + namespace internal + { + inline void check_ranges_consistency(const std::vector &ranges) + { + for (uint i = 1; i < ranges.size(); ++i) + if (!is_close(ranges[i].min, ranges[i - 1].max)) + throw std::runtime_error("Your range definition is inconsistent!"); + } + + template void append_range(T &to_append, const J &range) + { + for (const auto &element : range) { + to_append.push_back(element); + } + } + } // namespace internal + + template class ConfigurationMesh + { + public: + ConfigurationMesh() = delete; + ConfigurationMesh(const DiFfRG::JSONValue &json) + { + static_assert(dim >= 1 && dim <= 3, "ConfigurationMesh only supports dim = 1, 2, 3"); + + const std::array grid_names = {"/discretization/grid/x_grid", "/discretization/grid/y_grid", + "/discretization/grid/z_grid"}; + + for (int i = 0; i < dim; ++i) { + auto subranges = internal::string_to_substrings_array(json.get_string(grid_names[i])); + for (const auto &subrange_str : subranges) + grids[i].emplace_back(subrange_str); + internal::check_ranges_consistency(grids[i]); + } + + refine = json.get_uint("/discretization/grid/refine", 0); + } + + static std::string get_defaults() + { + std::string json_str = R"({ + "discretization": { + "grid": { + "x_grid": ")"; + json_str += GridAxis::get_default_range(); + json_str += R"(")"; + + if constexpr (dim >= 2) { + json_str += R"(, + "y_grid": ")"; + json_str += GridAxis::get_default_range(); + json_str += R"(")"; + } + + if constexpr (dim >= 3) { + json_str += R"(, + "z_grid": ")"; + json_str += GridAxis::get_default_range(); + json_str += R"(")"; + } + + json_str += R"(, + "refine": 0 + } + } +})"; + return json_str; + } + + inline std::vector> get_step_withs_for_triangulation() const + { + std::vector> step_sizes(dim); + for (int i = 0; i < dim; ++i) { + for (const auto &sub_grid_axis : grids[i]) { + internal::append_range(step_sizes[i], sub_grid_axis.get_stepwiths()); + } + } + return step_sizes; + } + + inline dealii::Point get_lower_left() const + { + dealii::Point lower_left; + for (int i = 0; i < dim; ++i) { + lower_left[i] = grids[i][0].min; + } + return lower_left; + } + + inline dealii::Point get_upper_right() const + { + dealii::Point upper_right; + for (int i = 0; i < dim; ++i) { + upper_right[i] = grids[i].back().max; + } + return upper_right; + } + + std::array, dim> grids; + uint refine; + }; + + } // namespace Config +} // namespace DiFfRG \ No newline at end of file diff --git a/DiFfRG/include/DiFfRG/discretization/mesh/rectangular_mesh.hh b/DiFfRG/include/DiFfRG/discretization/mesh/rectangular_mesh.hh index c8681a9bd..886f435fa 100644 --- a/DiFfRG/include/DiFfRG/discretization/mesh/rectangular_mesh.hh +++ b/DiFfRG/include/DiFfRG/discretization/mesh/rectangular_mesh.hh @@ -5,6 +5,7 @@ // DiFfRG #include +#include namespace DiFfRG { @@ -26,15 +27,18 @@ namespace DiFfRG * @brief Construct a new RectangularMesh object. * @param json JSONValue object containing the parameters for the mesh. */ + [[deprecated("Please use RectangularMesh(const Config::ConfigurationMesh &mesh_config) instead")]] RectangularMesh(const JSONValue &json); + RectangularMesh(const Config::ConfigurationMesh &mesh_config); + Triangulation &get_triangulation() { return triangulation; } const Triangulation &get_triangulation() const { return triangulation; } protected: virtual void make_grid(); - const JSONValue &json; + const Config::ConfigurationMesh &mesh_config; Triangulation triangulation; }; } // namespace DiFfRG \ No newline at end of file diff --git a/DiFfRG/src/DiFfRG/discretization/mesh/rectangular_mesh.cc b/DiFfRG/src/DiFfRG/discretization/mesh/rectangular_mesh.cc index 4403b2656..b8e9671eb 100644 --- a/DiFfRG/src/DiFfRG/discretization/mesh/rectangular_mesh.cc +++ b/DiFfRG/src/DiFfRG/discretization/mesh/rectangular_mesh.cc @@ -1,4 +1,5 @@ // standard library +#include "DiFfRG/discretization/mesh/configuration_mesh.hh" #include // external libraries @@ -7,7 +8,6 @@ #include // DiFfRG -#include #include namespace DiFfRG @@ -22,7 +22,14 @@ namespace DiFfRG void append_range(std::vector> &dest, const std::array &range, const uint d); } // namespace internal - template RectangularMesh::RectangularMesh(const JSONValue &json) : json(json) + template + RectangularMesh::RectangularMesh(const JSONValue &json) + : DiFfRG::RectangularMesh(Config::ConfigurationMesh(json)) + { + } + + template + RectangularMesh::RectangularMesh(const Config::ConfigurationMesh &mesh_config) : mesh_config(mesh_config) { static_assert(dim > 0 && dim < 4); make_grid(); @@ -30,54 +37,12 @@ namespace DiFfRG template void RectangularMesh::make_grid() { - std::vector> step_sizes; - std::vector> x_grid_ranges; - std::vector> y_grid_ranges; - std::vector> z_grid_ranges; - try { - auto x_grid_ranges_str = internal::string_to_substrings_array(json.get_string("/discretization/grid/x_grid")); - for (const auto &str : x_grid_ranges_str) - x_grid_ranges.emplace_back(internal::string_to_range(str)); - internal::check_ranges_consistency(x_grid_ranges); - for (const auto &range : x_grid_ranges) - internal::append_range(step_sizes, range, 0); - if constexpr (dim > 1) { - auto y_grid_ranges_str = internal::string_to_substrings_array(json.get_string("/discretization/grid/y_grid")); - for (const auto &str : y_grid_ranges_str) - y_grid_ranges.emplace_back(internal::string_to_range(str)); - internal::check_ranges_consistency(y_grid_ranges); - for (const auto &range : y_grid_ranges) - internal::append_range(step_sizes, range, 1); - } - if constexpr (dim > 2) { - auto z_grid_ranges_str = internal::string_to_substrings_array(json.get_string("/discretization/grid/z_grid")); - for (const auto &str : z_grid_ranges_str) - z_grid_ranges.emplace_back(internal::string_to_range(str)); - internal::check_ranges_consistency(z_grid_ranges); - for (const auto &range : z_grid_ranges) - internal::append_range(step_sizes, range, 2); - } - } catch (std::exception &e) { - std::cerr << "Error in grid construction: " + std::string(e.what()); - spdlog::get("log")->error("Error in grid construction: {}", e.what()); - } - - if constexpr (dim == 1) { - dealii::Point<1> origin{x_grid_ranges[0][0]}; - dealii::Point<1> extent{x_grid_ranges[x_grid_ranges.size() - 1][2]}; - dealii::GridGenerator::subdivided_hyper_rectangle(triangulation, step_sizes, origin, extent); - } else if constexpr (dim == 2) { - dealii::Point<2> origin{x_grid_ranges[0][0], y_grid_ranges[0][0]}; - dealii::Point<2> extent{x_grid_ranges[x_grid_ranges.size() - 1][2], y_grid_ranges[y_grid_ranges.size() - 1][2]}; - dealii::GridGenerator::subdivided_hyper_rectangle(triangulation, step_sizes, origin, extent); - } else if constexpr (dim == 3) { - dealii::Point<3> origin{x_grid_ranges[0][0], y_grid_ranges[0][0], z_grid_ranges[0][0]}; - dealii::Point<3> extent{x_grid_ranges[x_grid_ranges.size() - 1][2], y_grid_ranges[y_grid_ranges.size() - 1][2], - z_grid_ranges[z_grid_ranges.size() - 1][2]}; - dealii::GridGenerator::subdivided_hyper_rectangle(triangulation, step_sizes, origin, extent); - } + const auto lower_left = mesh_config.get_lower_left(); + const auto upper_right = mesh_config.get_upper_right(); + const auto step_sizes = mesh_config.get_step_withs_for_triangulation(); + dealii::GridGenerator::subdivided_hyper_rectangle(triangulation, step_sizes, lower_left, upper_right); - triangulation.refine_global(json.get_uint("/discretization/grid/refine")); + triangulation.refine_global(mesh_config.refine); if constexpr (dim == 2) { std::ofstream out("grid.svg"); @@ -86,54 +51,6 @@ namespace DiFfRG } } - namespace internal - { - std::array string_to_range(const std::string str) - { - std::vector vec; - std::istringstream ss(str); - std::string buf; - while (std::getline(ss, buf, ':')) - vec.push_back(std::stod(buf)); - if (vec.size() != 3) throw std::runtime_error(str + " is not a range"); - - std::array arr; - std::copy_n(vec.begin(), 3, arr.begin()); - return arr; - } - - std::vector string_to_substrings_array(const std::string str) - { - std::vector array; - std::istringstream ss(str); - std::string buf; - while (std::getline(ss, buf, ',')) - array.push_back(buf); - return array; - } - - void check_ranges_consistency(const std::vector> &ranges) - { - for (uint i = 1; i < ranges.size(); ++i) - if (!is_close(ranges[i][0], ranges[i - 1][2])) - throw std::runtime_error("Your range definition is inconsistent!"); - } - - void append_range(std::vector> &dest, const std::array &range, const uint d) - { - if (dest.size() <= d) dest.resize(d + 1); - - uint steps = uint((range[2] - range[0]) / range[1]); - for (uint i = 0; i < steps; ++i) - dest[d].push_back(range[1]); - - if (!is_close(range[0] + range[1] * steps, range[2])) { - if (range[0] + range[1] * steps > range[2]) throw std::runtime_error("grid construction broke down!"); - dest[d].push_back(range[2] - (range[0] + range[1] * steps)); - } - } - } // namespace internal - template class RectangularMesh<1>; template class RectangularMesh<2>; template class RectangularMesh<3>; diff --git a/DiFfRG/tests/discretization/CMakeLists.txt b/DiFfRG/tests/discretization/CMakeLists.txt index c2c77c36c..caf12d999 100644 --- a/DiFfRG/tests/discretization/CMakeLists.txt +++ b/DiFfRG/tests/discretization/CMakeLists.txt @@ -14,3 +14,4 @@ add_subdirectory(coordinates) setup_test(mesh.cc) setup_test(mesh_adaptor.cc) setup_test(EoM_finding.cc) +setup_test(mesh_configuration_test.cpp) diff --git a/DiFfRG/tests/discretization/mesh_configuration_test.cpp b/DiFfRG/tests/discretization/mesh_configuration_test.cpp new file mode 100644 index 000000000..b8311b849 --- /dev/null +++ b/DiFfRG/tests/discretization/mesh_configuration_test.cpp @@ -0,0 +1,223 @@ +// Tests for Mesh Configuration + +#include "DiFfRG/common/json.hh" +#include "catch2/catch_test_macros.hpp" +#include "catch2/generators/catch_generators.hpp" +#include "catch2/generators/catch_generators_adapters.hpp" +#include +#include +#include +#include + +using namespace DiFfRG::Config; + +// // Test case for default constructor +TEST_CASE("Mesh Configuration exists", "[MeshConfiguration]") +{ + SECTION("initializes correctly the mesh config") + { + + // triangulation.refine_global(json.get_uint("/discretization/grid/refine")); + DiFfRG::JSONValue json = DiFfRG::json::value({ + {"discretization", + { + {"grid", {{"x_grid", "0:0.1:1"}, {"y_grid", "0:0.1:1,1:0.2:2"}, {"z_grid", "0:0.1:1"}}}, + }}, + }); + ConfigurationMesh<3> mesh_config(json); + REQUIRE(size(mesh_config.grids[0]) == 1); + REQUIRE(size(mesh_config.grids[1]) == 2); + REQUIRE(size(mesh_config.grids[2]) == 1); + } + + SECTION("refine parameter works correctly") + { + + DiFfRG::JSONValue json_with_refine = DiFfRG::json::value({ + {"discretization", + { + {"grid", {{"x_grid", "0:0.1:1"}, {"refine", 1}}}, + }}, + }); + ConfigurationMesh<1> mesh_config_with_refine(json_with_refine); + CHECK(mesh_config_with_refine.refine == 1); + DiFfRG::JSONValue json_without_refine = DiFfRG::json::value({ + {"discretization", + { + {"grid", {{"x_grid", "0:0.1:1"}}}, + }}, + }); + ConfigurationMesh<1> mesh_config_without_refine(json_without_refine); + CHECK(mesh_config_without_refine.refine == 0); + } + + SECTION("works only with x_axis") + { + DiFfRG::JSONValue json = DiFfRG::json::value({ + {"discretization", + { + {"grid", {{"x_grid", "0:0.1:1"}}}, + }}, + }); + ConfigurationMesh<1> mesh_config(json); + REQUIRE(size(mesh_config.grids[0]) == 1); + } + SECTION("check range consistency") + { + DiFfRG::JSONValue json = DiFfRG::json::value({ + {"discretization", + { + {"grid", {{"x_grid", "0:0.1:1,1.1:0.2:2"}}}, + }}, + }); + CHECK_THROWS(ConfigurationMesh<1>(json)); + } + SECTION("check a correct step size array") + { + DiFfRG::JSONValue json = DiFfRG::json::value({ + {"discretization", + { + {"grid", {{"x_grid", "0:0.5:1,1:0.6:2"}}}, + }}, + }); + ConfigurationMesh<1> mesh_config(json); + REQUIRE(size(mesh_config.get_step_withs_for_triangulation()[0]) == 4); + REQUIRE(mesh_config.get_step_withs_for_triangulation()[0] == std::vector{0.5, 0.5, 0.5, 0.5}); + } + SECTION("check a correct step size array in 2D") + { + DiFfRG::JSONValue json = DiFfRG::json::value({ + {"discretization", + { + {"grid", {{"x_grid", "0:0.5:1,1:0.6:2"}, {"y_grid", "0:1.0:2,2:1.0:4"}}}, + }}, + }); + ConfigurationMesh<2> mesh_config(json); + REQUIRE(size(mesh_config.get_step_withs_for_triangulation()[0]) == 4); + REQUIRE(mesh_config.get_step_withs_for_triangulation()[0] == std::vector{0.5, 0.5, 0.5, 0.5}); + REQUIRE(size(mesh_config.get_step_withs_for_triangulation()[1]) == 4); + REQUIRE(mesh_config.get_step_withs_for_triangulation()[1] == std::vector(4, 1.0)); + } + SECTION("check a correct step size array in 3D") + { + DiFfRG::JSONValue json = DiFfRG::json::value({ + {"discretization", + { + {"grid", {{"x_grid", "0:0.5:1,1:0.6:2"}, {"y_grid", "0:1.0:2,2:1.0:4"}, {"z_grid", "0:2.0:2,2:2.0:4"}}}, + }}, + }); + ConfigurationMesh<3> mesh_config(json); + REQUIRE(size(mesh_config.get_step_withs_for_triangulation()[0]) == 4); + REQUIRE(mesh_config.get_step_withs_for_triangulation()[0] == std::vector{0.5, 0.5, 0.5, 0.5}); + REQUIRE(size(mesh_config.get_step_withs_for_triangulation()[1]) == 4); + REQUIRE(mesh_config.get_step_withs_for_triangulation()[1] == std::vector(4, 1.0)); + REQUIRE(size(mesh_config.get_step_withs_for_triangulation()[2]) == 2); + REQUIRE(mesh_config.get_step_withs_for_triangulation()[2] == std::vector(2, 2.0)); + } + + SECTION("Get lower_left and upper_right of grid in 3D", "[MeshConfiguration]") + { + using dealii::Point; + DiFfRG::JSONValue json = DiFfRG::json::value({ + {"discretization", + { + {"grid", + {{"x_grid", "0:0.5:1,1:0.6:2"}, {"y_grid", "0.1:1.0:2,2:1.0:4"}, {"z_grid", "0.2:0.2:2,2:2.0:4"}}}, + }}, + }); + ConfigurationMesh<3> mesh_config(json); + CHECK(mesh_config.get_lower_left() == Point<3, double>(0.0, 0.1, 0.2)); + CHECK(mesh_config.get_upper_right() == Point<3, double>(2.0, 4.0, 4.0)); + + DiFfRG::JSONValue json2 = DiFfRG::json::value({ + {"discretization", + { + {"grid", + {{"x_grid", "0.2:0.5:1,1:0.6:3"}, {"y_grid", "0.3:1.0:2,2:1.0:6"}, {"z_grid", "0.4:0.2:2,2:2.0:8"}}}, + }}, + }); + ConfigurationMesh<3> mesh_config2(json2); + CHECK(mesh_config2.get_lower_left() == Point<3, double>(0.2, 0.3, 0.4)); + CHECK(mesh_config2.get_upper_right() == Point<3, double>(3.0, 6.0, 8.0)); + } +} + +TEST_CASE("GridAxis tests", "[MeshConfiguration]") +{ + SECTION("parses range_string \"0:0.1:1\" correctly") + { + GridAxis grid_axis("0:0.1:1"); + REQUIRE(grid_axis.min == 0.0); + REQUIRE(grid_axis.step == 0.1); + REQUIRE(grid_axis.max == 1.0); + } + + SECTION("parses range_string \"0.1:0.1:2\" correctly") + { + GridAxis grid_axis("0.1:0.1:2"); + REQUIRE(grid_axis.min == 0.1); + REQUIRE(grid_axis.step == 0.1); + REQUIRE(grid_axis.max == 2.0); + } + + SECTION("returns array of steps") + { + GridAxis grid_axis("0.0:0.5:2"); + REQUIRE(size(grid_axis.get_stepwiths()) == 4); + REQUIRE(grid_axis.get_stepwiths() == std::vector{0.5, 0.5, 0.5, 0.5}); + } + + SECTION("returns array of steps if step does not add up exactly to max") + { + GridAxis grid_axis("0.0:0.6:2"); + REQUIRE(size(grid_axis.get_stepwiths()) == 4); + REQUIRE(grid_axis.get_stepwiths() == std::vector{0.5, 0.5, 0.5, 0.5}); + } + + SECTION("returns array of steps if step does not add up exactly to max, with one step") + { + GridAxis grid_axis("0:2.0:2"); + REQUIRE(size(grid_axis.get_stepwiths()) == 1); + REQUIRE(grid_axis.get_stepwiths() == std::vector{2.0}); + } +} + +TEST_CASE("throws when range is invalid", "[MeshConfiguration]") +{ + CHECK_THROWS(GridAxis("0.1:0.01:0.0")); + CHECK_THROWS(GridAxis("0.1:2.0:1.0")); +} + +TEST_CASE("parses default JSON grid configuration", "[MeshConfiguration]") +{ + SECTION("get_defaults returns valid parseable JSON") + { + std::string default_json_str = ConfigurationMesh<3>::get_defaults(); + DiFfRG::JSONValue default_json = DiFfRG::json::parse(default_json_str); + + // Parse the defaults into a ConfigurationMesh + ConfigurationMesh<3> mesh_config(default_json); + + // Verify it matches expected default configuration - all grids use the same default range + REQUIRE(size(mesh_config.grids[0]) == 3); + CHECK(mesh_config.grids[0][0].min == 0.0); + CHECK(mesh_config.grids[0][0].max == 7e-3); + CHECK(mesh_config.grids[0][1].min == 7e-3); + CHECK(mesh_config.grids[0][1].max == 9e-3); + CHECK(mesh_config.grids[0][2].min == 9e-3); + CHECK(mesh_config.grids[0][2].max == 1.1e-2); + + REQUIRE(size(mesh_config.grids[1]) == 3); + CHECK(mesh_config.grids[1][0].min == 0.0); + CHECK(mesh_config.grids[1][0].max == 7e-3); + + REQUIRE(size(mesh_config.grids[2]) == 3); + CHECK(mesh_config.grids[2][0].min == 0.0); + CHECK(mesh_config.grids[2][0].max == 7e-3); + + CHECK(mesh_config.refine == 0); + + CHECK(mesh_config.get_lower_left() == dealii::Point<3, double>(0.0, 0.0, 0.0)); + CHECK(mesh_config.get_upper_right() == dealii::Point<3, double>(1.1e-2, 1.1e-2, 1.1e-2)); + } +} \ No newline at end of file