refactor: new trial set algebra

demos/FiniteVolume/CMakeLists.txt

@@ -17,6 +17,7 @@ level_set.cpp:finite-volume-level-set
 level_set_from_scratch.cpp:finite-volume-level-set-from-scratch
 advection_1d.cpp:finite-volume-advection-1d
 advection_2d.cpp:finite-volume-advection-2d
+advection_3d.cpp:finite-volume-advection-3d
 advection_2d_user_bc.cpp:finite-volume-advection-2d-user-bc
 scalar_burgers_2d.cpp:finite-volume-scalar-burgers-2d
 linear_convection_obstacle.cpp:finite-volume-linear-convection-obstacle

demos/FiniteVolume/advection_3d.cpp

@@ -0,0 +1,165 @@
+// Copyright 2018-2025 the samurai's authors
+// SPDX-License-Identifier:  BSD-3-Clause
+
+#include <array>
+
+#include <xtensor/xfixed.hpp>
+
+#include <samurai/algorithm.hpp>
+#include <samurai/bc.hpp>
+#include <samurai/field.hpp>
+#include <samurai/io/hdf5.hpp>
+#include <samurai/io/restart.hpp>
+#include <samurai/mr/adapt.hpp>
+#include <samurai/mr/mesh.hpp>
+#include <samurai/samurai.hpp>
+#include <samurai/stencil_field.hpp>
+#include <samurai/subset/node.hpp>
+
+#include <filesystem>
+namespace fs = std::filesystem;
+
+template <class Field>
+void init(Field& u)
+{
+    auto& mesh = u.mesh();
+    u.resize();
+
+    samurai::for_each_cell(mesh,
+                           [&](auto& cell)
+                           {
+                               auto center           = cell.center();
+                               const double radius   = .2;
+                               const double x_center = 0.3;
+                               const double y_center = 0.3;
+                               const double z_center = 0.3;
+                               if (((center[0] - x_center) * (center[0] - x_center) + (center[1] - y_center) * (center[1] - y_center)
+                                    + (center[2] - z_center) * (center[2] - z_center))
+                                   <= radius * radius)
+                               {
+                                   u[cell] = 1;
+                               }
+                               else
+                               {
+                                   u[cell] = 0;
+                               }
+                           });
+}
+
+template <class Field>
+void save(const fs::path& path, const std::string& filename, const Field& u, const std::string& suffix = "")
+{
+    auto& mesh = u.mesh();
+
+#ifdef SAMURAI_WITH_MPI
+    mpi::communicator world;
+    samurai::save(path, fmt::format("{}_size_{}{}", filename, world.size(), suffix), mesh, u);
+#else
+    samurai::save(path, fmt::format("{}{}", filename, suffix), mesh, u);
+    samurai::dump(path, fmt::format("{}_restart{}", filename, suffix), mesh, u);
+#endif
+}
+
+int main(int argc, char* argv[])
+{
+    auto& app = samurai::initialize("Finite volume example for the advection equation in 3d using multiresolution", argc, argv);
+
+    constexpr std::size_t dim = 3;
+    using Config              = samurai::MRConfig<dim>;
+
+    // Simulation parameters
+    xt::xtensor_fixed<double, xt::xshape<dim>> min_corner = {0., 0., 0.};
+    xt::xtensor_fixed<double, xt::xshape<dim>> max_corner = {1., 1., 1.};
+    std::array<double, dim> a{
+        {1, 1, 1}
+    };
+    double Tf  = .1;
+    double cfl = 0.25;
+    double t   = 0.;
+    std::string restart_file;
+
+    // Multiresolution parameters
+    std::size_t min_level = 4;
+    std::size_t max_level = 10;
+
+    // Output parameters
+    fs::path path        = fs::current_path();
+    std::string filename = "FV_advection_3d";
+    std::size_t nfiles   = 1;
+
+    app.add_option("--min-corner", min_corner, "The min corner of the box")->capture_default_str()->group("Simulation parameters");
+    app.add_option("--max-corner", max_corner, "The max corner of the box")->capture_default_str()->group("Simulation parameters");
+    app.add_option("--velocity", a, "The velocity of the advection equation")->capture_default_str()->group("Simulation parameters");
+    app.add_option("--cfl", cfl, "The CFL")->capture_default_str()->group("Simulation parameters");
+    app.add_option("--Ti", t, "Initial time")->capture_default_str()->group("Simulation parameters");
+    app.add_option("--Tf", Tf, "Final time")->capture_default_str()->group("Simulation parameters");
+    app.add_option("--restart-file", restart_file, "Restart file")->capture_default_str()->group("Simulation parameters");
+    app.add_option("--min-level", min_level, "Minimum level of the multiresolution")->capture_default_str()->group("Multiresolution");
+    app.add_option("--max-level", max_level, "Maximum level of the multiresolution")->capture_default_str()->group("Multiresolution");
+    app.add_option("--path", path, "Output path")->capture_default_str()->group("Output");
+    app.add_option("--filename", filename, "File name prefix")->capture_default_str()->group("Output");
+    app.add_option("--nfiles", nfiles, "Number of output files")->capture_default_str()->group("Output");
+
+    SAMURAI_PARSE(argc, argv);
+
+    const samurai::Box<double, dim> box(min_corner, max_corner);
+    samurai::MRMesh<Config> mesh;
+    auto u = samurai::make_scalar_field<double>("u", mesh);
+
+    if (restart_file.empty())
+    {
+        mesh = {
+            box,
+            min_level,
+            max_level,
+            {true, true, true}
+        };
+        init(u);
+    }
+    else
+    {
+        samurai::load(restart_file, mesh, u);
+    }
+    // samurai::make_bc<samurai::Dirichlet<1>>(u, 0.);
+
+    double dt            = cfl * mesh.cell_length(max_level);
+    const double dt_save = Tf / static_cast<double>(nfiles);
+
+    auto unp1 = samurai::make_scalar_field<double>("unp1", mesh);
+
+    auto MRadaptation = samurai::make_MRAdapt(u);
+    auto mra_config   = samurai::mra_config().epsilon(2e-4);
+    MRadaptation(mra_config);
+    save(path, filename, u, "_init");
+
+    std::size_t nsave = 1;
+    std::size_t nt    = 0;
+
+    while (t != Tf)
+    {
+        MRadaptation(mra_config);
+
+        t += dt;
+        if (t > Tf)
+        {
+            dt += Tf - t;
+            t = Tf;
+        }
+
+        std::cout << fmt::format("iteration {}: t = {}, dt = {}", nt++, t, dt) << std::endl;
+
+        samurai::update_ghost_mr(u);
+        unp1.resize();
+        unp1 = u - dt * samurai::upwind(a, u);
+
+        std::swap(u.array(), unp1.array());
+
+        if (t >= static_cast<double>(nsave + 1) * dt_save || t == Tf)
+        {
+            const std::string suffix = (nfiles != 1) ? fmt::format("_ite_{}", nsave++) : "";
+            save(path, filename, u, suffix);
+        }
+    }
+    samurai::finalize();
+    return 0;
+}

demos/FiniteVolume/linear_convection_obstacle.cpp

@@ -136,7 +136,8 @@ int main(int argc, char* argv[])
     }
 
     double t = 0;
-    while (t != Tf)
+    //~ while (t != Tf)
+    for (int ite = 0; ite != 23; ++ite)
     {
         // Move to next timestep
         t += dt;

demos/tutorial/set_operator.cpp

@@ -98,7 +98,7 @@ int main()
                                u[cell] = cell.indices[0];
                            });
 
-    auto subset1 = samurai::intersection(ca[0], samurai::contraction(ca[1], 1));
+    auto subset1 = samurai::intersection(ca[0], samurai::contract(ca[1], 1));
     subset1.on(0)(
         [&](const auto& i, auto)
         {

include/samurai/algorithm.hpp

@@ -14,6 +14,9 @@
 #include "cell.hpp"
 #include "mesh_holder.hpp"
 
+#include "concepts.hpp"
+#include "subset/node.hpp"
+
 using namespace xt::placeholders;
 
 namespace samurai
@@ -111,18 +114,28 @@ namespace samurai
         }
     }
 
-    template <class Mesh, class Func>
+    template <IsMesh Mesh, class Func>
     inline void for_each_interval(const Mesh& mesh, Func&& f)
     {
         using mesh_id_t = typename Mesh::config::mesh_id_t;
         for_each_interval(mesh[mesh_id_t::cells], std::forward<Func>(f));
     }
 
-    template <class Op, class StartEndOp, class... S>
-    class Subset;
-
-    template <class Func, class Op, class StartEndOp, class... S>
-    inline void for_each_interval(Subset<Op, StartEndOp, S...>& set, Func&& f)
+    //~ template <class Op, class StartEndOp, class... S>
+    //~ class Subset;
+
+    //~ template <class Func, class Op, class StartEndOp, class... S>
+    //~ inline void for_each_interval(Subset<Op, StartEndOp, S...>& set, Func&& f)
+    //~ {
+    //~     set(
+    //~         [&](const auto& i, const auto& index)
+    //~         {
+    //~             f(set.level(), i, index);
+    //~         });
+    //~ }
+
+    template <class Func, class Set>
+    inline void for_each_interval(const SetBase<Set>& set, Func&& f)
     {
         set(
             [&](const auto& i, const auto& index)

include/samurai/boundary.hpp

@@ -9,7 +9,7 @@ namespace samurai
     {
         using mesh_id_t = typename Mesh::mesh_id_t;
 
-        auto& cells = mesh[mesh_id_t::cells][level];
+        const auto& cells = mesh[mesh_id_t::cells][level];
 
         return difference(cells, translate(self(domain).on(level), -layer_width * direction));
     }
@@ -43,7 +43,7 @@ namespace samurai
     {
         using mesh_id_t = typename Mesh::mesh_id_t;
 
-        auto& cells = mesh[mesh_id_t::cells][level];
+        const auto& cells = mesh[mesh_id_t::cells][level];
 
         return difference(cells, contract(self(mesh.domain()).on(level), 1));
     }

include/samurai/interval.hpp

@@ -362,7 +362,8 @@ namespace samurai
     template <class value_t, class index_t>
     inline bool operator==(const Interval<value_t, index_t>& i1, const Interval<value_t, index_t>& i2)
     {
-        return !(i1.start != i2.start || i1.end != i2.end || i1.step != i2.step || i1.index != i2.index);
+        //~ return !(i1.start != i2.start || i1.end != i2.end || i1.step != i2.step || i1.index != i2.index);
+        return !(i1.start != i2.start || i1.end != i2.end || i1.step != i2.step);
     }
 
     template <class value_t, class index_t>
@@ -376,6 +377,17 @@ namespace samurai
     {
         return i1.start < i2.start;
     }
+
+    template <typename T>
+    struct IsInterval : std::false_type
+    {
+    };
+
+    template <class value_t, class index_t>
+    struct IsInterval<Interval<value_t, index_t>> : std::true_type
+    {
+    };
+
 } // namespace samurai
 
 template <class TValue, class TIndex>

include/samurai/level_cell_array.hpp

@@ -94,8 +94,11 @@ namespace samurai
         LevelCellArray() = default;
         LevelCellArray(const LevelCellList<Dim, TInterval>& lcl);
 
-        template <class Op, class StartEndOp, class... S>
-        LevelCellArray(Subset<Op, StartEndOp, S...> set);
+        template <class Set>
+        explicit LevelCellArray(const SetBase<Set>& set);
+
+        template <class Set>
+        LevelCellArray(const SetBase<Set>& set, const coords_t& origin_point, const double scaling_factor);
 
         LevelCellArray(std::size_t level, const Box<value_t, dim>& box);
         LevelCellArray(std::size_t level,
@@ -344,9 +347,21 @@ namespace samurai
         }
     }
 
+    //~ template <std::size_t Dim, class TInterval>
+    //~ template <class Op, class StartEndOp, class... S>
+    //~ inline LevelCellArray<Dim, TInterval>::LevelCellArray(Subset<Op, StartEndOp, S...> set)
+    //~     : m_level(set.level())
+    //~ {
+    //~     set(
+    //~         [this](const auto& i, const auto& index)
+    //~         {
+    //~             add_interval_back(i, index);
+    //~         });
+    //~ }
+
     template <std::size_t Dim, class TInterval>
-    template <class Op, class StartEndOp, class... S>
-    inline LevelCellArray<Dim, TInterval>::LevelCellArray(Subset<Op, StartEndOp, S...> set)
+    template <class Set>
+    inline LevelCellArray<Dim, TInterval>::LevelCellArray(const SetBase<Set>& set)
         : m_level(set.level())
     {
         set(
@@ -356,6 +371,32 @@ namespace samurai
             });
     }
 
+    template <std::size_t Dim, class TInterval>
+    template <class Set>
+    inline LevelCellArray<Dim, TInterval>::LevelCellArray(const SetBase<Set>& set, const coords_t& origin_point, const double scaling_factor)
+        : m_level(set.level())
+        , m_origin_point(origin_point)
+        , m_scaling_factor(scaling_factor)
+    {
+        set(
+            [this](const auto& i, const auto& index)
+            {
+                add_interval_back(i, index);
+            });
+    }
+
+    //~ template <std::size_t Dim, class TInterval>
+    //~ template <class Op, class StartEndOp, class... S>
+    //~ inline LevelCellArray<Dim, TInterval>::LevelCellArray(Subset<Op, StartEndOp, S...> set)
+    //~     : m_level(set.level())
+    //~ {
+    //~     set(
+    //~         [this](const auto& i, const auto& index)
+    //~         {
+    //~             add_interval_back(i, index);
+    //~         });
+    //~ }
+
     template <std::size_t Dim, class TInterval>
     inline LevelCellArray<Dim, TInterval>::LevelCellArray(std::size_t level, const Box<value_t, dim>& box)
         : m_level{level}

include/samurai/list_of_intervals.hpp

@@ -62,6 +62,11 @@ namespace samurai
 
         void add_point(value_t point);
         void add_interval(const interval_t& interval);
+
+        void clear()
+        {
+            std::forward_list<Interval<TValue, TIndex>>::clear();
+        }
     };
 
     ////////////////////////////////////

include/samurai/mesh.hpp

@@ -815,16 +815,18 @@ namespace samurai
             {
                 if constexpr (dim == 2)
                 {
-                    m_corners.push_back(difference(
-                        m_domain,
-                        union_(translate(m_domain, direction_t{-direction[0], 0}), translate(m_domain, direction_t{0, -direction[1]}))));
+                    m_corners.push_back(difference(m_domain,
+                                                   union_(translate(m_domain, direction_t{-direction[0], 0}),
+                                                          translate(m_domain, direction_t{0, -direction[1]})))
+                                            .to_lca());
                 }
                 else if constexpr (dim == 3)
                 {
                     m_corners.push_back(difference(m_domain,
                                                    union_(translate(m_domain, direction_t{-direction[0], 0, 0}),
                                                           translate(m_domain, direction_t{0, -direction[1], 0}),
-                                                          translate(m_domain, direction_t{0, 0, -direction[2]}))));
+                                                          translate(m_domain, direction_t{0, 0, -direction[2]})))
+                                            .to_lca());
                 }
             });
     }

include/samurai/samurai_config.hpp

@@ -20,6 +20,9 @@ namespace samurai
         static constexpr std::size_t graduation_width = 1;
         static constexpr std::size_t prediction_order = 1;
 
+        static constexpr bool projection_with_list_of_intervals = true;
+        static constexpr bool use_native_expand                 = true;
+
         using index_t    = signed long long int;
         using value_t    = int;
         using interval_t = Interval<value_t, index_t>;