Refactoring of RectangularMesh and Config Loading

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;

DiFfRG/include/DiFfRG/discretization/mesh/configuration_mesh.hh

@@ -0,0 +1,189 @@
+#pragma once
+#include "DiFfRG/common/math.hh"
+#include <DiFfRG/common/utils.hh>
+#include <stdexcept>
+#include <string>
+#include <vector>
+
+namespace DiFfRG
+{
+
+  namespace Config
+  {
+
+    namespace internal
+    {
+      inline std::vector<std::string> string_to_substrings_array(const std::string str)
+      {
+        std::vector<std::string> array;
+        std::istringstream ss(str);
+        std::string buf;
+        while (std::getline(ss, buf, ','))
+          array.push_back(buf);
+        return array;
+      }
+      inline std::array<double, 3> string_to_range(const std::string str)
+      {
+        std::vector<double> 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<double, 3> 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<double> get_stepwiths() const
+      {
+        double local_step = step;
+        auto steps = static_cast<uint>((max - min) / local_step);
+        if (!is_close((min + local_step * steps), max)) {
+          local_step = (max - min) / (steps + 1);
+          ++steps;
+        }
+        return std::vector<double>(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<GridAxis> &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 <typename T, typename J> void append_range(T &to_append, const J &range)
+      {
+        for (const auto &element : range) {
+          to_append.push_back(element);
+        }
+      }
+    } // namespace internal
+
+    template <int dim> 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<std::string, 3> 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<std::vector<double>> get_step_withs_for_triangulation() const
+      {
+        std::vector<std::vector<double>> 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<dim, double> get_lower_left() const
+      {
+        dealii::Point<dim, double> lower_left;
+        for (int i = 0; i < dim; ++i) {
+          lower_left[i] = grids[i][0].min;
+        }
+        return lower_left;
+      }
+
+      inline dealii::Point<dim, double> get_upper_right() const
+      {
+        dealii::Point<dim, double> upper_right;
+        for (int i = 0; i < dim; ++i) {
+          upper_right[i] = grids[i].back().max;
+        }
+        return upper_right;
+      }
+
+      std::array<std::vector<GridAxis>, dim> grids;
+      uint refine;
+    };
+
+  } // namespace Config
+} // namespace DiFfRG

DiFfRG/include/DiFfRG/discretization/mesh/rectangular_mesh.hh

@@ -5,6 +5,7 @@
 
 // DiFfRG
 #include <DiFfRG/common/utils.hh>
+#include <DiFfRG/discretization/mesh/configuration_mesh.hh>
 
 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<dim> &mesh_config) instead")]]
     RectangularMesh(const JSONValue &json);
 
+    RectangularMesh(const Config::ConfigurationMesh<dim> &mesh_config);
+
     Triangulation<dim> &get_triangulation() { return triangulation; }
     const Triangulation<dim> &get_triangulation() const { return triangulation; }
 
   protected:
     virtual void make_grid();
 
-    const JSONValue &json;
+    const Config::ConfigurationMesh<dim> &mesh_config;
     Triangulation<dim> triangulation;
   };
 } // namespace DiFfRG

DiFfRG/src/DiFfRG/discretization/mesh/rectangular_mesh.cc

@@ -1,4 +1,5 @@
 // standard library
+#include "DiFfRG/discretization/mesh/configuration_mesh.hh"
 #include <fstream>
 
 // external libraries
@@ -7,7 +8,6 @@
 #include <deal.II/grid/grid_refinement.h>
 
 // DiFfRG
-#include <DiFfRG/common/utils.hh>
 #include <DiFfRG/discretization/mesh/rectangular_mesh.hh>
 
 namespace DiFfRG
@@ -22,62 +22,27 @@ namespace DiFfRG
     void append_range(std::vector<std::vector<double>> &dest, const std::array<double, 3> &range, const uint d);
   } // namespace internal
 
-  template <uint dim> RectangularMesh<dim>::RectangularMesh(const JSONValue &json) : json(json)
+  template <uint dim>
+  RectangularMesh<dim>::RectangularMesh(const JSONValue &json)
+      : DiFfRG::RectangularMesh<dim>(Config::ConfigurationMesh<dim>(json))
+  {
+  }
+
+  template <uint dim>
+  RectangularMesh<dim>::RectangularMesh(const Config::ConfigurationMesh<dim> &mesh_config) : mesh_config(mesh_config)
   {
     static_assert(dim > 0 && dim < 4);
     make_grid();
   }
 
   template <uint dim> void RectangularMesh<dim>::make_grid()
   {
-    std::vector<std::vector<double>> step_sizes;
-    std::vector<std::array<double, 3>> x_grid_ranges;
-    std::vector<std::array<double, 3>> y_grid_ranges;
-    std::vector<std::array<double, 3>> 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<double, 3> string_to_range(const std::string str)
-    {
-      std::vector<double> 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<double, 3> arr;
-      std::copy_n(vec.begin(), 3, arr.begin());
-      return arr;
-    }
-
-    std::vector<std::string> string_to_substrings_array(const std::string str)
-    {
-      std::vector<std::string> 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<std::array<double, 3>> &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<std::vector<double>> &dest, const std::array<double, 3> &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>;

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)