diff --git a/demos/FiniteVolume/CMakeLists.txt b/demos/FiniteVolume/CMakeLists.txt index 8ff4b2bf9..f117398c3 100644 --- a/demos/FiniteVolume/CMakeLists.txt +++ b/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 diff --git a/demos/FiniteVolume/advection_3d.cpp b/demos/FiniteVolume/advection_3d.cpp new file mode 100644 index 000000000..c69731471 --- /dev/null +++ b/demos/FiniteVolume/advection_3d.cpp @@ -0,0 +1,165 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +namespace fs = std::filesystem; + +template +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 +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; + + // Simulation parameters + xt::xtensor_fixed> min_corner = {0., 0., 0.}; + xt::xtensor_fixed> max_corner = {1., 1., 1.}; + std::array 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 box(min_corner, max_corner); + samurai::MRMesh mesh; + auto u = samurai::make_scalar_field("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>(u, 0.); + + double dt = cfl * mesh.cell_length(max_level); + const double dt_save = Tf / static_cast(nfiles); + + auto unp1 = samurai::make_scalar_field("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(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; +} diff --git a/demos/FiniteVolume/linear_convection_obstacle.cpp b/demos/FiniteVolume/linear_convection_obstacle.cpp index a39c513e0..f8329ed83 100644 --- a/demos/FiniteVolume/linear_convection_obstacle.cpp +++ b/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; diff --git a/demos/tutorial/set_operator.cpp b/demos/tutorial/set_operator.cpp index 4f46700b1..82a4bd53e 100644 --- a/demos/tutorial/set_operator.cpp +++ b/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) { diff --git a/include/samurai/algorithm.hpp b/include/samurai/algorithm.hpp index bd6c21617..9c262b4a4 100644 --- a/include/samurai/algorithm.hpp +++ b/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 + template 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(f)); } - template - class Subset; - - template - inline void for_each_interval(Subset& set, Func&& f) + //~ template + //~ class Subset; + + //~ template + //~ inline void for_each_interval(Subset& set, Func&& f) + //~ { + //~ set( + //~ [&](const auto& i, const auto& index) + //~ { + //~ f(set.level(), i, index); + //~ }); + //~ } + + template + inline void for_each_interval(const SetBase& set, Func&& f) { set( [&](const auto& i, const auto& index) diff --git a/include/samurai/boundary.hpp b/include/samurai/boundary.hpp index ac505f65e..b22d62ae2 100644 --- a/include/samurai/boundary.hpp +++ b/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)); } diff --git a/include/samurai/interval.hpp b/include/samurai/interval.hpp index 97ef0cd5d..c84a60263 100644 --- a/include/samurai/interval.hpp +++ b/include/samurai/interval.hpp @@ -362,7 +362,8 @@ namespace samurai template inline bool operator==(const Interval& i1, const Interval& 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 diff --git a/include/samurai/level_cell_array.hpp b/include/samurai/level_cell_array.hpp index 48b79c717..18b8c1eba 100644 --- a/include/samurai/level_cell_array.hpp +++ b/include/samurai/level_cell_array.hpp @@ -94,8 +94,11 @@ namespace samurai LevelCellArray() = default; LevelCellArray(const LevelCellList& lcl); - template - LevelCellArray(Subset set); + template + explicit LevelCellArray(const SetBase& set); + + template + LevelCellArray(const SetBase& set, const coords_t& origin_point, const double scaling_factor); LevelCellArray(std::size_t level, const Box& box); LevelCellArray(std::size_t level, @@ -344,9 +347,21 @@ namespace samurai } } + //~ template + //~ template + //~ inline LevelCellArray::LevelCellArray(Subset set) + //~ : m_level(set.level()) + //~ { + //~ set( + //~ [this](const auto& i, const auto& index) + //~ { + //~ add_interval_back(i, index); + //~ }); + //~ } + template - template - inline LevelCellArray::LevelCellArray(Subset set) + template + inline LevelCellArray::LevelCellArray(const SetBase& set) : m_level(set.level()) { set( @@ -356,6 +371,32 @@ namespace samurai }); } + template + template + inline LevelCellArray::LevelCellArray(const SetBase& 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 + //~ template + //~ inline LevelCellArray::LevelCellArray(Subset set) + //~ : m_level(set.level()) + //~ { + //~ set( + //~ [this](const auto& i, const auto& index) + //~ { + //~ add_interval_back(i, index); + //~ }); + //~ } + template inline LevelCellArray::LevelCellArray(std::size_t level, const Box& box) : m_level{level} diff --git a/include/samurai/list_of_intervals.hpp b/include/samurai/list_of_intervals.hpp index d3889f847..ee8d90e12 100644 --- a/include/samurai/list_of_intervals.hpp +++ b/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>::clear(); + } }; //////////////////////////////////// diff --git a/include/samurai/mesh.hpp b/include/samurai/mesh.hpp index 1fa2192f7..ed17be689 100644 --- a/include/samurai/mesh.hpp +++ b/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()); } }); } diff --git a/include/samurai/samurai_config.hpp b/include/samurai/samurai_config.hpp index b4394f61a..fc8f32d87 100644 --- a/include/samurai/samurai_config.hpp +++ b/include/samurai/samurai_config.hpp @@ -20,6 +20,8 @@ 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 = false; + using index_t = signed long long int; using value_t = int; using interval_t = Interval; diff --git a/include/samurai/subset/apply.hpp b/include/samurai/subset/apply.hpp index b8652afe2..aef13b38d 100644 --- a/include/samurai/subset/apply.hpp +++ b/include/samurai/subset/apply.hpp @@ -3,119 +3,97 @@ #pragma once -#include "concepts.hpp" -#include "utils.hpp" +#include "set_base.hpp" namespace samurai { namespace detail { - template - bool apply_impl(Set&& global_set, Func&& func, Container& index) + template + void apply_rec(const SetBase& set, + Func&& func, + typename SetBase::yz_index_t& yz_index, + std::integral_constant d_ic, + typename SetBase::Workspace& workspace) { - auto set = global_set.template get_local_set(global_set.level(), index); - auto start_and_stop = global_set.template get_start_and_stop_function(); - - if constexpr (dim != 1) + using traverser_t = typename Set::template traverser_t; + using current_interval_t = typename traverser_t::current_interval_t; + using interval_t = typename traverser_t::interval_t; + + set.init_workspace(1, d_ic, workspace); + + //~for (traverser_t traverser = set.get_traverser(yz_index, d_ic, workspace); !traverser.is_empty(); traverser.next_interval()) + //~{ + //~ current_interval_t interval = traverser.current_interval(); + //~ + //~ if constexpr (d == 0) + //~ { + //~ func(interval, yz_index); + //~ } + //~ else + //~ { + //~ for (yz_index[d - 1] = interval.start; yz_index[d - 1] != interval.end; ++yz_index[d - 1]) + //~ { + //~ apply_rec(set, std::forward(func), yz_index, std::integral_constant{}, workspace); + //~ } + //~ } + //~} + traverser_t traverser = set.get_traverser(yz_index, d_ic, workspace); + if (!traverser.is_empty()) { - auto func_int = [&](const auto& interval) + interval_t interval = traverser.current_interval(); + + if constexpr (d == 0) + { + func(interval, yz_index); + } + else { - for (auto i = interval.start; i < interval.end; ++i) + for (yz_index[d - 1] = interval.start; yz_index[d - 1] != interval.end; ++yz_index[d - 1]) { - index[dim - 2] = i; - if (apply_impl(std::forward(global_set), std::forward(func), index)) - { - return true; - } + apply_rec(set, std::forward(func), yz_index, std::integral_constant{}, workspace); } - return false; - }; - return apply(set, start_and_stop, func_int); - } - else - { - auto func_int = [&](const auto& interval) - { - return func(interval, index); - }; - return apply(set, start_and_stop, func_int); - } - } - } - - template - void apply(Set&& global_set, Func&& user_func) - { - constexpr std::size_t dim = std::decay_t::dim; - xt::xtensor_fixed> index; + } + traverser.next_interval(); - auto func = [&](const auto& interval, const auto& yz) - { - user_func(interval, yz); - return false; - }; + while (!traverser.is_empty()) + { + interval_t old_interval = interval; - if (global_set.exist()) - { - detail::apply_impl(std::forward(global_set), func, index); - } - } + interval = traverser.current_interval(); - template - bool empty_check(Set&& global_set) - { - constexpr std::size_t dim = std::decay_t::dim; - xt::xtensor_fixed> index; + assert(old_interval < interval); - auto func = [](const auto&, const auto&) - { - return true; - }; - - if (global_set.exist()) - { - return !detail::apply_impl(std::forward(global_set), func, index); + if constexpr (d == 0) + { + func(interval, yz_index); + } + else + { + for (yz_index[d - 1] = interval.start; yz_index[d - 1] != interval.end; ++yz_index[d - 1]) + { + apply_rec(set, std::forward(func), yz_index, std::integral_constant{}, workspace); + } + } + traverser.next_interval(); + } + } } - return true; } - template - requires IsSetOp || IsIntervalListVisitor - bool apply(Set&& set, StartEnd&& start_and_stop, Func&& func) + template + void apply(const SetBase& set, Func&& func) { - using interval_t = typename std::decay_t::interval_t; - using value_t = typename interval_t::value_t; + using Workspace = typename Set::Workspace; + using yz_index_t = typename Set::yz_index_t; - interval_t result; - int r_ipos = 0; - set.next(0, std::forward(start_and_stop)); - auto scan = set.min(); + constexpr std::size_t dim = Set::dim; - while (scan < sentinel && !set.is_empty()) + yz_index_t yz_index; + if (set.exist()) { - bool is_in = set.is_in(scan); - - if (is_in && r_ipos == 0) - { - result.start = scan; - r_ipos = 1; - } - else if (!is_in && r_ipos == 1) - { - result.end = scan; - r_ipos = 0; - - auto true_result = set.shift() >= 0 ? result >> static_cast(set.shift()) - : result << -static_cast(set.shift()); - if (func(true_result)) - { - return true; - } - } - - set.next(scan, std::forward(start_and_stop)); - scan = set.min(); + Workspace workspace; + detail::apply_rec(set, std::forward(func), yz_index, std::integral_constant{}, workspace); } - return false; } } diff --git a/include/samurai/subset/box_view.hpp b/include/samurai/subset/box_view.hpp new file mode 100644 index 000000000..6ed4f10e5 --- /dev/null +++ b/include/samurai/subset/box_view.hpp @@ -0,0 +1,89 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "set_base.hpp" +#include "traversers/box_traverser.hpp" + +namespace samurai +{ + + template + class BoxView; + + template + struct SetTraits> + { + static_assert(std::same_as, B>); + + template + using traverser_t = BoxTraverser; + + struct Workspace + { + }; + + static constexpr std::size_t dim() + { + return B::dim; + } + }; + + template + class BoxView : public SetBase> + { + using Self = BoxView; + + public: + + SAMURAI_SET_TYPEDEFS + + BoxView(const std::size_t level, const B& box) + : m_level(level) + , m_box(box) + { + } + + inline std::size_t level_impl() const + { + return m_level; + } + + inline bool exist_impl() const + { + return m_box.is_valid(); + } + + inline bool empty_impl() const + { + return !exist_impl(); + } + + template + inline traverser_t get_traverser_impl(const yz_index_t& index, std::integral_constant, Workspace) const + { + return (m_box.min_corner()[d + 1] <= index[d] && index[d] < m_box.max_corner()[d + 1]) + ? traverser_t(m_box.min_corner()[d], m_box.max_corner()[d]) + : traverser_t(0, 0); + } + + template + inline traverser_t + get_traverser_unordered_impl(const yz_index_t& index, std::integral_constant d_ic, Workspace) const + { + return get_traverser_impl(index, d_ic, Workspace{}); + } + + template + inline constexpr void init_workspace_impl(const std::size_t, std::integral_constant, Workspace) const + { + } + + private: + + std::size_t m_level; + const B& m_box; + }; + +} // namespace samurai diff --git a/include/samurai/subset/concepts.hpp b/include/samurai/subset/concepts.hpp deleted file mode 100644 index 4fcdc6125..000000000 --- a/include/samurai/subset/concepts.hpp +++ /dev/null @@ -1,44 +0,0 @@ -// Copyright 2018-2025 the samurai's authors -// SPDX-License-Identifier: BSD-3-Clause - -#pragma once - -#include -#include - -namespace samurai -{ - template - class LevelCellArray; - // } - - // namespace samurai::experimental - // { - template - class SetTraverser; - - template - class IntervalListVisitor; - - template - struct is_setop : std::false_type - { - }; - - template - struct is_setop> : std::true_type - { - }; - - template - constexpr bool is_setop_v{is_setop>::value}; - - template - concept IsSetOp = is_setop_v; - - template - concept IsIntervalListVisitor = std::is_base_of_v::container_t>, std::decay_t>; - - template - concept IsLCA = std::same_as, T>; -} diff --git a/include/samurai/subset/contraction.hpp b/include/samurai/subset/contraction.hpp new file mode 100644 index 000000000..83e394293 --- /dev/null +++ b/include/samurai/subset/contraction.hpp @@ -0,0 +1,138 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "set_base.hpp" +#include "traversers/contraction_traverser.hpp" + +namespace samurai +{ + + template + class Contraction; + + template + struct SetTraits> + { + static_assert(IsSet::value); + + template + using traverser_t = ContractionTraverser>; + + struct Workspace + { + typename Set::Workspace child_workspace; + }; + + static constexpr std::size_t dim() + { + return Set::dim; + } + }; + + template + class Contraction : public SetBase> + { + using Self = Contraction; + + public: + + SAMURAI_SET_TYPEDEFS + + using contraction_t = std::array; + using do_contraction_t = std::array; + + Contraction(const Set& set, const contraction_t& contraction) + : m_set(set) + , m_contraction(contraction) + { + assert(std::all_of(m_contraction.cbegin(), + m_contraction.cend(), + [](const contraction_t& c) + { + return c >= 0; + })); + } + + Contraction(const Set& set, const value_t contraction) + : m_set(set) + { + assert(contraction >= 0); + std::fill(m_contraction.begin(), m_contraction.end(), contraction); + } + + Contraction(const Set& set, const value_t contraction, const do_contraction_t& do_contraction) + : m_set(set) + { + for (std::size_t i = 0; i != m_contraction.size(); ++i) + { + m_contraction[i] = contraction * do_contraction[i]; + } + } + + inline std::size_t level_impl() const + { + return m_set.level(); + } + + inline bool exist_impl() const + { + return m_set.exist(); + } + + inline bool empty_impl() const + { + return Base::empty_default_impl(); + } + + template + inline void + init_workspace_impl(const std::size_t n_traversers, std::integral_constant d_ic, Workspace& workspace) const + { + m_set.init_workspace(n_traversers, d_ic, workspace.child_workspace); + } + + template + inline traverser_t + get_traverser_impl(const yz_index_t& index, std::integral_constant d_ic, Workspace& workspace) const + { + return traverser_t(m_set.get_traverser(index, d_ic, workspace.child_workspace), m_contraction[d]); + } + + template + inline traverser_t + get_traverser_unordered_impl(const yz_index_t& index, std::integral_constant d_ic, Workspace& workspace) const + { + return traverser_t(m_set.get_traverser_unordered(index, d_ic, workspace.child_workspace), m_contraction[d]); + } + + private: + + Set m_set; + contraction_t m_contraction; + }; + + template + auto contract(const Set& set, const typename Contraction>::contraction_t& contraction) + { + return Contraction(self(set), contraction); + } + + template + auto contract(const Set& set, const typename Contraction>::value_t& contraction) + { + return Contraction(self(set), contraction); + } + + template + auto contract(const Set& set, + const typename Contraction>::value_t& contraction, + const typename Contraction>::do_contraction_t& do_contraction) // idk how to make this + // more readable, + // perhaps a traits... + { + return Contraction(self(set), contraction, do_contraction); + } + +} // namespace samurai diff --git a/include/samurai/subset/expansion.hpp b/include/samurai/subset/expansion.hpp new file mode 100644 index 000000000..326bb7b63 --- /dev/null +++ b/include/samurai/subset/expansion.hpp @@ -0,0 +1,205 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "set_base.hpp" +#include "traversers/expansion_traverser.hpp" +#include "traversers/last_dim_expansion_traverser.hpp" + +namespace samurai +{ + + namespace detail + { + template + struct ExpansionWorkspace; + + template + struct ExpansionWorkspace> + { + template + using child_traverser_t = typename Set::template traverser_t; + + template + using array_of_child_traverser_t = std::vector>; + + using Type = std::tuple...>; + + static_assert(std::tuple_size::value == Set::dim - 1); + }; + + } // namespace detail + + template + class Expansion; + + template + struct SetTraits> + { + static_assert(IsSet::value); + + template + using traverser_t = std::conditional_t>, + ExpansionTraverser>>; + + struct Workspace + { + typename detail::ExpansionWorkspace>::Type expansion_workspace; + typename Set::Workspace child_workspace; + }; + + static constexpr std::size_t dim() + { + return Set::dim; + } + }; + + template + class Expansion : public SetBase> + { + using Self = Expansion; + + public: + + SAMURAI_SET_TYPEDEFS + + using expansion_t = std::array; + using do_expansion_t = std::array; + + Expansion(const Set& set, const expansion_t& expansions) + : m_set(set) + , m_expansions(expansions) + { + } + + Expansion(const Set& set, const value_t expansion) + : m_set(set) + { + m_expansions.fill(expansion); + } + + Expansion(const Set& set, const value_t expansion, const do_expansion_t& do_expansion) + : m_set(set) + { + for (std::size_t i = 0; i != m_expansions.size(); ++i) + { + m_expansions[i] = expansion * do_expansion[i]; + } + } + + inline std::size_t level_impl() const + { + return m_set.level(); + } + + inline bool exist_impl() const + { + return m_set.exist(); + } + + inline bool empty_impl() const + { + return m_set.empty(); + } + + template + inline void + init_workspace_impl(const std::size_t n_traversers, std::integral_constant d_ic, Workspace& workspace) const + { + if constexpr (d == Base::dim - 1) + { + m_set.init_workspace(n_traversers, d_ic, workspace.child_workspace); + } + else + { + const std::size_t my_work_size = n_traversers * 2 * std::size_t(m_expansions[d + 1] + 1); + + auto& childTraversers = std::get(workspace.expansion_workspace); + childTraversers.clear(); + childTraversers.reserve(my_work_size); + + m_set.init_workspace(my_work_size, d_ic, workspace.child_workspace); + } + } + + template + inline traverser_t + get_traverser_impl(const yz_index_t& index, std::integral_constant d_ic, Workspace& workspace) const + { + if constexpr (d == Base::dim - 1) + { + return traverser_t(m_set.get_traverser(index, d_ic, workspace.child_workspace), m_expansions[d]); + } + else + { + auto& childTraversers = std::get(workspace.expansion_workspace); + + yz_index_t tmp_index(index); + + const auto childTraversers_begin = childTraversers.end(); + + for (value_t width = 0; width != m_expansions[d + 1] + 1; ++width) + { + // it's d and not d+1 because tmp_index represents m_expansions[1,..,dim] + tmp_index[d] = index[d] + width; + childTraversers.push_back(m_set.get_traverser_unordered(tmp_index, d_ic, workspace.child_workspace)); + if (childTraversers.back().is_empty()) + { + childTraversers.pop_back(); + } + // same + tmp_index[d] = index[d] - width; + childTraversers.push_back(m_set.get_traverser_unordered(tmp_index, d_ic, workspace.child_workspace)); + if (childTraversers.back().is_empty()) + { + childTraversers.pop_back(); + } + } + + return traverser_t(childTraversers_begin, childTraversers.end(), m_expansions[d]); + } + } + + template + inline traverser_t + get_traverser_unordered_impl(const yz_index_t& index, std::integral_constant d_ic, Workspace& workspace) const + { + if constexpr (d == Base::dim - 1) + { + return traverser_t(m_set.get_traverser_unordered(index, d_ic, workspace.child_workspace), m_expansions[d]); + } + else + { + return get_traverser_impl(index, d_ic, workspace); + } + } + + private: + + Set m_set; + expansion_t m_expansions; + }; + + template + auto expand(const Set& set, const typename Contraction>::contraction_t& expansions) + { + return Expansion(self(set), expansions); + } + + template + auto expand(const Set& set, const typename Contraction>::value_t& expansion) + { + return Expansion(self(set), expansion); + } + + template + auto expand(const Set& set, + const typename Contraction>::value_t& expansion, + const typename Contraction>::do_expansion_t& do_expansion) + { + return Expansion(self(set), expansion, do_expansion); + } + +} // namespace samurai diff --git a/include/samurai/subset/lca_view.hpp b/include/samurai/subset/lca_view.hpp new file mode 100644 index 000000000..153656f3d --- /dev/null +++ b/include/samurai/subset/lca_view.hpp @@ -0,0 +1,196 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "set_base.hpp" +#include "traversers/lca_traverser.hpp" + +namespace samurai +{ + + template + class LCAView; + + template + struct SetTraits> + { + static_assert(std::same_as, LCA>); + + template + using traverser_t = LCATraverser; + + struct Workspace + { + Workspace() + { + start_offset_guess.fill(0); + } + + // we do not need the offsets for the last dim + // the offset at dimension d will be initialized when calling + // get_traverser_impl + std::array start_offset; + std::array end_offset; + + std::array start_offset_guess; + }; + + static constexpr std::size_t dim() + { + return LCA::dim; + } + }; + + template + class LCAView : public SetBase> + { + using Self = LCAView; + + using const_interval_iterator = typename std::vector::const_iterator; + + public: + + SAMURAI_SET_TYPEDEFS + + explicit LCAView(const LCA& lca) + : m_lca(lca) + { + } + + inline std::size_t level_impl() const + { + return m_lca.level(); + } + + inline bool exist_impl() const + { + return !empty_impl(); + } + + inline bool empty_impl() const + { + return m_lca.empty(); + } + + template + inline constexpr void init_workspace_impl(const std::size_t, std::integral_constant, Workspace&) const + { + } + + template + inline traverser_t get_traverser_impl(const yz_index_t& index, std::integral_constant, Workspace& workspace) const + { + if constexpr (d == Base::dim - 1) + { + return traverser_t(m_lca[d].cbegin(), m_lca[d].cend()); + } + else + { + // In 3d, we would be in the y dimension + // we need to find an interval that contains the prescibed z. + const auto& z_intervals = m_lca[d + 1]; + const auto begin_z_interval = (d == Base::dim - 2) ? z_intervals.cbegin() + : z_intervals.cbegin() + workspace.start_offset[d + 1]; + const auto end_z_interval = (d == Base::dim - 2) ? z_intervals.cend() : z_intervals.cbegin() + workspace.end_offset[d + 1]; + + const auto z = index[d]; + + const auto z_interval_it = std::find_if(begin_z_interval, + end_z_interval, + [z](const auto& z_interval) + { + return z_interval.contains(z); + }); + + const auto& y_intervals = m_lca[d]; + + auto& y_start_offset = workspace.start_offset[d]; + auto& y_end_offset = workspace.end_offset[d]; + + if (z_interval_it == end_z_interval) + { + y_start_offset = y_end_offset; + return traverser_t(y_intervals.cend(), y_intervals.cend()); + } + const auto& y_offsets = m_lca.offsets(d + 1); + const auto y_offset_idx = std::size_t(z_interval_it->index + z); + + y_start_offset = std::ptrdiff_t(y_offsets[y_offset_idx]); + y_end_offset = std::ptrdiff_t(y_offsets[y_offset_idx + 1]); + return traverser_t(y_intervals.cbegin() + y_start_offset, y_intervals.cbegin() + y_end_offset); + } + } + + template + inline traverser_t + get_traverser_unordered_impl(const yz_index_t& index, std::integral_constant, Workspace& workspace) const + { + ptrdiff_t start_offset = 0; + ptrdiff_t end_offset = std::ssize(m_lca[Base::dim - 1]); + + for (std::size_t dim = Base::dim - 1; dim != d; --dim) + { + const auto y = index[dim - 1]; + const auto& y_intervals = m_lca[dim]; + auto& start_offset_guess = workspace.start_offset_guess[dim - 1]; + + const auto begin_y_intervals = y_intervals.cbegin() + start_offset; + const auto end_y_intervals = y_intervals.cbegin() + end_offset; + const auto y_intervals_size = std::distance(begin_y_intervals, end_y_intervals); + // if guess was wrong for the higer dimensions, the guess is wrong an may even be out of [begin_y_intervals, + // end_y_intervals) we thus need to cap start_offset to ensure begin_y_intervals <= begin_y_intervals_guess <= + // end_y_intervals + const auto begin_y_intervals_guess = begin_y_intervals + std::min(start_offset, y_intervals_size); + + assert(begin_y_intervals <= begin_y_intervals_guess and begin_y_intervals_guess <= end_y_intervals); + + // we know the interval that contains y is likely to be in the range [begin_y_intervals_guess, end_y_intervals) + // first we try to find it within [begin_y_intervals_guess, end_y_intervals) + // hopefully, *begin_y_intervals_guess contains y. + auto y_interval_it = std::find_if(begin_y_intervals_guess, + end_y_intervals, + [y](const auto& y_interval) + { + return y_interval.contains(y); + }); + if (y_interval_it == end_y_intervals) + { + // we did not find an interval that contains y in [begin_y_intervals_guess, end_y_intervals) + // we try to find it in [begin_y_intervals, begin_y_intervals_guess) + y_interval_it = std::find_if(begin_y_intervals, + begin_y_intervals_guess, + [y](const auto& y_interval) + { + return y_interval.contains(y); + }); + if (y_interval_it == begin_y_intervals_guess) + { + // there is no interval that contains y + return traverser_t(m_lca[d].cend(), m_lca[d].cend()); + } + } + start_offset_guess = std::distance(begin_y_intervals, y_interval_it); + assert(0 <= start_offset_guess and start_offset_guess < std::distance(begin_y_intervals, end_y_intervals)); + + const auto& y_offsets = m_lca.offsets(dim); + const auto y_offset_idx = std::size_t(y + y_interval_it->index); + + start_offset = ptrdiff_t(y_offsets[y_offset_idx]); + end_offset = ptrdiff_t(y_offsets[y_offset_idx + 1]); + } + return traverser_t(m_lca[d].cbegin() + start_offset, m_lca[d].cbegin() + end_offset); + } + + private: + + const LCA& m_lca; + }; + + template + LCAView> self(const LevelCellArray& lca) + { + return LCAView>(lca); + } + +} // namespace samurai diff --git a/include/samurai/subset/nary_set_operator.hpp b/include/samurai/subset/nary_set_operator.hpp new file mode 100644 index 000000000..b9b99329f --- /dev/null +++ b/include/samurai/subset/nary_set_operator.hpp @@ -0,0 +1,232 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "set_base.hpp" +#include "traversers/difference_id_traverser.hpp" +#include "traversers/difference_traverser.hpp" +#include "traversers/intersection_traverser.hpp" +#include "traversers/union_traverser.hpp" +#include "utils.hpp" + +namespace samurai +{ + enum class SetOperator + { + UNION, + INTERSECTION, + DIFFERENCE + }; + + template + class NArySetOperator; + + template + struct SetTraits> + { + static_assert((IsSet::value and ...)); + + template + using traverser_t = UnionTraverser...>; + + struct Workspace + { + std::tuple children_workspaces; + }; + + static constexpr std::size_t dim() + { + return std::tuple_element_t<0, std::tuple>::dim; + } + }; + + template + struct SetTraits> + { + static_assert((IsSet::value and ...)); + + template + using traverser_t = IntersectionTraverser...>; + + struct Workspace + { + std::tuple children_workspaces; + }; + + static constexpr std::size_t dim() + { + return std::tuple_element_t<0, std::tuple>::dim; + } + }; + + template + struct SetTraits> + { + static_assert((IsSet::value and ...)); + + template + using traverser_t = std::conditional_t...>, + DifferenceIdTraverser...>>; + + struct Workspace + { + std::tuple children_workspaces; + }; + + static constexpr std::size_t dim() + { + return std::tuple_element_t<0, std::tuple>::dim; + } + }; + + template + class NArySetOperator : public SetBase> + { + using Self = NArySetOperator; + + public: + + SAMURAI_SET_TYPEDEFS + + using Childrens = std::tuple; + + static constexpr std::size_t nIntervals = std::tuple_size_v; + + explicit NArySetOperator(const Sets&... sets) + : m_sets(sets...) + { + m_level = std::apply( + [](const auto&... set) -> std::size_t + { + return vmax(set.level()...); + }, + m_sets); + + enumerate_const_items(m_sets, + [this](const auto i, const auto& set) + { + m_shifts[i] = std::size_t(m_level - set.level()); + }); + } + + inline std::size_t level_impl() const + { + return m_level; + } + + inline bool exist_impl() const + { + return std::apply( + [](const auto first_set, const auto&... other_sets) -> std::size_t + { + if constexpr (op == SetOperator::UNION) + { + return first_set.exist() || (other_sets.exist() || ...); + } + else if constexpr (op == SetOperator::INTERSECTION) + { + return first_set.exist() && (other_sets.exist() && ...); + } + else + { + return first_set.exist(); + } + }, + m_sets); + } + + inline bool empty_impl() const + { + return Base::empty_default_impl(); + } + + template + inline void + init_workspace_impl(const std::size_t n_traversers, std::integral_constant d_ic, Workspace& workspace) const + { + static_for<0, nIntervals>::apply( + [this, n_traversers, d_ic, &workspace](const auto i) -> void + { + std::get(m_sets).init_workspace(n_traversers, d_ic, std::get(workspace.children_workspaces)); + }); + } + + template + inline traverser_t + get_traverser_impl(const yz_index_t& index, std::integral_constant d_ic, Workspace& workspace) const + { + return get_traverser_impl_detail(index, d_ic, std::make_index_sequence{}, workspace); + } + + template + inline traverser_t + get_traverser_unordered_impl(const yz_index_t& index, std::integral_constant d_ic, Workspace& workspace) const + { + return get_traverser_unordered_impl_detail(index, d_ic, std::make_index_sequence{}, workspace); + } + + private: + + template + traverser_t get_traverser_impl_detail(const yz_index_t& index, + std::integral_constant d_ic, + std::index_sequence, + Workspace& workspace) const + { + return traverser_t(m_shifts, + std::get(m_sets).get_traverser(utils::powMinus2(index, m_shifts[Is]), + d_ic, + std::get(workspace.children_workspaces))...); + } + + template + traverser_t get_traverser_unordered_impl_detail(const yz_index_t& index, + std::integral_constant d_ic, + std::index_sequence, + Workspace& workspace) const + { + return traverser_t(m_shifts, + std::get(m_sets).get_traverser_unordered(utils::powMinus2(index, m_shifts[Is]), + d_ic, + std::get(workspace.children_workspaces))...); + } + + Childrens m_sets; + std::size_t m_level; + std::array m_shifts; + }; + + //////////////////////////////////////////////////////////////////////// + //// functions + //////////////////////////////////////////////////////////////////////// + + template + auto union_(const Sets&... sets) + requires(sizeof...(Sets) >= 1) + { + using Union = NArySetOperator...>; + + return Union(self(sets)...); + } + + template + auto intersection(const Sets&... sets) + requires(sizeof...(Sets) >= 2) + { + using Intersection = NArySetOperator...>; + + return Intersection(self(sets)...); + } + + template + auto difference(const Sets&... sets) + requires(sizeof...(Sets) >= 2) + { + using Difference = NArySetOperator...>; + + return Difference(self(sets)...); + } + +} // namespace samurai diff --git a/include/samurai/subset/node.hpp b/include/samurai/subset/node.hpp index 6114bf0f4..d36aa8009 100644 --- a/include/samurai/subset/node.hpp +++ b/include/samurai/subset/node.hpp @@ -3,687 +3,13 @@ #pragma once -#include -#include -#include -#include -#include -#include - -#include - -#include "../algorithm.hpp" #include "apply.hpp" -#include "concepts.hpp" -#include "samurai/list_of_intervals.hpp" -#include "start_end_fct.hpp" -#include "visitor.hpp" - -namespace samurai -{ - - template - void apply(Set&& global_set, Func&& func); - - template - class Subset - { - public: - - static constexpr std::size_t dim = get_set_dim_v; - using set_type = std::tuple; - using interval_t = get_interval_t; - - Subset(Op&& op, StartEndOp&& start_end_op, S&&... s) - : m_operator(std::forward(op)) - , m_start_end_op(std::forward(start_end_op)) - , m_s(std::forward(s)...) - , m_ref_level(compute_max(s.ref_level()...)) - , m_level(compute_max(s.level()...)) - , m_min_level(m_level) - { - std::apply( - [this](auto&&... args) - { - (args.ref_level(m_ref_level), ...); - }, - m_s); - } - - auto& on(auto level) - { - auto ilevel = static_cast(level); - if (ilevel > m_ref_level) - { - ref_level(ilevel); - } - m_min_level = std::min(m_min_level, ilevel); - m_level = ilevel; - return *this; - } - - template - void operator()(Func&& func) - { - apply(*this, std::forward(func)); - } - - template - void apply_op(ApplyOp&&... op) - { - auto func = [&](auto& interval, auto& index) - { - (op(m_level, interval, index), ...); - }; - apply(*this, func); - } - - inline void to_stream(std::ostream& os) - { - apply_op( - [&](const auto level, const auto& i, const auto& index) - { - os << "level: " << level << ", i: " << i << ", index: " << index << std::endl; - }); - } - - template - auto get_local_set(auto level, auto& index, Func_goback_beg&& goback_fct_beg, Func_goback_end&& goback_fct_end) - { - int shift = static_cast(this->ref_level()) - static_cast(this->level()); - m_start_end_op(m_level, m_min_level, m_ref_level); - - return std::apply( - [this, &index, shift, level, &goback_fct_beg, &goback_fct_end](auto&&... args) - { - return SetTraverser(shift, - get_operator(m_operator), - args.template get_local_set( - level, - index, - m_start_end_op.template goback(std::forward(goback_fct_beg)), - m_start_end_op.template goback(std::forward(goback_fct_end)))...); - }, - m_s); - } - - template - auto get_local_set(auto level, auto& index) - { - return get_local_set(level, index, default_function_(), default_function_()); - } - - template - auto get_start_and_stop_function(Func_start&& start_fct, Func_end&& end_fct) - { - m_start_end_op(m_level, m_min_level, m_ref_level); - - return std::apply( - [this, &start_fct, &end_fct](auto&& arg, auto&&... args) - { - if constexpr (std::is_same_v) - { - return std::make_tuple(std::move(arg.template get_start_and_stop_function( - m_start_end_op.template start(std::forward(start_fct)), - m_start_end_op.template end(std::forward(end_fct)))), - std::move(args.template get_start_and_stop_function( - m_start_end_op.template start(std::forward(start_fct)), - m_start_end_op.template end(std::forward(end_fct))))...); - } - else - { - return std::make_tuple(std::move(arg.template get_start_and_stop_function( - m_start_end_op.template start(std::forward(start_fct)), - m_start_end_op.template end(std::forward(end_fct)))), - std::move(args.template get_start_and_stop_function( - m_start_end_op.template start(std::forward(start_fct)), - m_start_end_op.template end(std::forward(end_fct))))...); - } - }, - m_s); - } - - template - auto get_start_and_stop_function() - { - return get_start_and_stop_function(default_function(), default_function()); - } - - auto level() const - { - return m_level; - } - - auto ref_level() const - { - return m_ref_level; - } - - void ref_level(auto level) - { - m_ref_level = level; - std::apply( - [this](auto&&... args) - { - (args.ref_level(m_ref_level), ...); - }, - m_s); - } - - bool empty() - { - return empty_check(*this); - } - - bool exist() const - { - return std::apply( - [this](auto&&... args) - { - return m_operator.exist(args...); - }, - m_s); - } - - protected: - - Op m_operator; - StartEndOp m_start_end_op; - set_type m_s; - std::size_t m_ref_level; - std::size_t m_level; - std::size_t m_min_level; - }; - - template - requires IsLCA - struct Self - { - static constexpr std::size_t dim = lca_t::dim; - using interval_t = typename lca_t::interval_t; - using value_t = typename interval_t::value_t; - - explicit Self(const lca_t& lca) - : m_lca(lca) - , m_level(lca.level()) - , m_ref_level(m_level) - , m_min_level(m_level) - { - } - - template - void operator()(Func&& func) - { - apply(*this, std::forward(func)); - } - - template - void apply_op(ApplyOp&&... op) - { - auto func = [&](auto& interval, auto& index) - { - (op(m_level, interval, index), ...); - }; - apply(*this, func); - } - - template - auto get_local_set(auto level, auto& index, Func_goback_beg&& goback_fct_beg, Func_goback_end&& goback_fct_end) - { - if (m_lca[d - 1].empty()) - { - return IntervalListVisitor(IntervalListRange(m_lca[d - 1], 0, 0)); - } - - if constexpr (dim == d) - { - m_offsets[d - 1].clear(); - m_offsets[d - 1].push_back({0, m_lca[d - 1].size()}); - - return IntervalListVisitor(m_lca.level(), - m_level, - m_ref_level, - IntervalListRange(m_lca[d - 1], 0, static_cast(m_lca[d - 1].size()))); - } - else - { - if (m_offsets[d].empty() || m_lca[d].empty()) - { - return IntervalListVisitor(IntervalListRange(m_lca[d - 1], 0, 0)); - } - - auto new_goback_fct_beg = m_func.template goback(std::forward(goback_fct_beg)); - - if (level <= m_level && level >= m_lca.level()) - { - m_offsets[d - 1].clear(); - - auto current_index = start_shift(new_goback_fct_beg(level, index[d - 1]).second, - static_cast(m_lca.level()) - static_cast(m_level)); - auto j = find_on_dim(m_lca, d, m_offsets[d][0][0], m_offsets[d][0][1], current_index); - - if (j == std::numeric_limits::max()) - { - return IntervalListVisitor(IntervalListRange(m_lca[d - 1], 0, 0)); - } - - auto io = static_cast(m_lca[d][j].index + current_index); - auto& offsets = m_lca.offsets(d); - m_offsets[d - 1].push_back({offsets[io], offsets[io + 1]}); - - return IntervalListVisitor(m_lca.level(), - m_level, - m_ref_level, - IntervalListRange(m_lca[d - 1], - static_cast(offsets[io]), - static_cast(offsets[io + 1]))); - } - else - { - auto new_goback_fct_end = m_func.template goback(std::forward(goback_fct_end)); - - auto min_index = start_shift(new_goback_fct_beg(level, index[d - 1]).second, - static_cast(m_lca.level()) - static_cast(m_level)); - - auto max_index = end_shift(new_goback_fct_end(level, index[d - 1] + 1).second, - static_cast(m_lca.level()) - static_cast(m_level)); - - m_work[d - 1].clear(); - m_offsets[d - 1].clear(); - - if constexpr (d == dim - 1) - { - auto j_min = lower_bound_interval(m_lca[d].begin() + static_cast(m_offsets[d][0][0]), - m_lca[d].begin() + static_cast(m_offsets[d][0][1]), - min_index); - auto j_max = upper_bound_interval(j_min, m_lca[d].begin() + static_cast(m_offsets[d][0][1]), max_index) - - 1; - - if (j_min != m_lca[d].end() && j_min <= j_max) - { - auto start_offset = static_cast(j_min->index + j_min->start); - if (j_min->contains(min_index)) - { - start_offset = static_cast(j_min->index + min_index); - } - - auto end_offset = static_cast(j_max->index + j_max->end); - if (j_max->contains(max_index)) - { - end_offset = static_cast(j_max->index + max_index); - } - - if (start_offset == end_offset) - { - return IntervalListVisitor(IntervalListRange(m_lca[d - 1], 0, 0)); - } - - m_offsets[d - 1].push_back({start_offset, end_offset}); - - ListOfIntervals list_of_intervals; - for (std::size_t o = m_lca.offsets(d)[start_offset]; o < m_lca.offsets(d)[end_offset]; ++o) - { - auto start = m_lca[d - 1][o].start; - auto end = m_lca[d - 1][o].end; - list_of_intervals.add_interval({start, end}); - } - - for (auto& i : list_of_intervals) - { - m_work[d - 1].push_back(i); - } - } - } - else - { - ListOfIntervals list_of_intervals; - - for (auto& offset : m_offsets[d]) - { - for (std::size_t io = offset[0]; io < offset[1]; ++io) - { - auto j_min = lower_bound_interval( - m_lca[d].begin() + static_cast(m_lca.offsets(d + 1)[io]), - m_lca[d].begin() + static_cast(m_lca.offsets(d + 1)[io + 1]), - min_index); - auto j_max = upper_bound_interval( - j_min, - m_lca[d].begin() + static_cast(m_lca.offsets(d + 1)[io + 1]), - max_index) - - 1; - - if (j_min != m_lca[d].begin() + static_cast(m_lca.offsets(d + 1)[io + 1]) && j_min <= j_max) - { - auto start_offset = static_cast(j_min->index + j_min->start); - if (j_min->contains(min_index)) - { - start_offset = static_cast(j_min->index + min_index); - } - - auto end_offset = static_cast(j_max->index + j_max->end); - if (j_max->contains(max_index)) - { - end_offset = static_cast(j_max->index + max_index); - } - - if (start_offset == end_offset) - { - continue; - } - - m_offsets[d - 1].push_back({start_offset, end_offset}); - - for (std::size_t o = m_lca.offsets(d)[start_offset]; o < m_lca.offsets(d)[end_offset]; ++o) - { - auto start = m_lca[d - 1][o].start; - auto end = m_lca[d - 1][o].end; - list_of_intervals.add_interval({start, end}); - } - } - } - - for (auto& i : list_of_intervals) - { - m_work[d - 1].push_back(i); - } - } - } - if (m_work[d - 1].empty()) - { - return IntervalListVisitor(IntervalListRange(m_lca[d - 1], 0, 0)); - } - return IntervalListVisitor(m_lca.level(), m_level, m_ref_level, IntervalListRange(m_lca[d - 1], m_work[d - 1])); - } - } - } - - template - auto get_local_set(auto level, auto& index) - - { - return get_local_set(level, index, default_function_(), default_function_()); - } - - template - auto get_start_and_stop_function(Func_start&& start_fct, Func_end&& end_fct) - { - m_func(m_level, m_min_level, m_ref_level); - auto new_start_fct = m_func.template start(std::forward(start_fct)); - auto new_end_fct = m_func.template end(std::forward(end_fct)); - return std::make_tuple(std::move(new_start_fct), std::move(new_end_fct)); - } - - template - auto get_start_and_stop_function() - - { - return get_start_and_stop_function(default_function(), default_function()); - } - - auto ref_level() const - { - return m_ref_level; - } - - void ref_level(auto level) - { - m_ref_level = level; - } - - auto level() const - { - return m_level; - } - - auto& on(auto level) - { - m_ref_level = std::max(m_ref_level, static_cast(level)); - m_min_level = std::min(m_min_level, static_cast(level)); - m_level = static_cast(level); - return *this; - } - - bool exist() const - { - return !m_lca.empty(); - } - - bool empty() const - { - return !m_lca.empty(); - } - - const lca_t& m_lca; - std::size_t m_level; - std::size_t m_ref_level; - std::size_t m_min_level; - start_end_function m_func; - std::array, dim - 1> m_work; - std::array>, dim> m_offsets; - }; - - namespace detail - { - template - auto transform(const LevelCellArray& lca) - { - return Self>(lca); - } - - template - auto transform(LevelCellArray& lca) - { - return Self>(lca); - } - - template - auto transform(E&& e) - { - return std::forward(e); - } - } - - template - inline std::ostream& operator<<(std::ostream& out, Subset& subset) - { - subset.to_stream(out); - return out; - } - - template - auto intersection(sets_t&&... sets) - { - static constexpr std::size_t dim = get_set_dim_v; - return std::apply( - [](auto&&... args) - { - return Subset(IntersectionOp(), start_end_function(), std::forward(args)...); - }, - std::make_tuple(detail::transform(std::forward(sets))...)); - } - - template - auto union_(sets_t&&... sets) - { - static constexpr std::size_t dim = get_set_dim_v; - return std::apply( - [](auto&&... args) - { - return Subset(UnionOp(), start_end_function(), std::forward(args)...); - }, - std::make_tuple(detail::transform(std::forward(sets))...)); - } - - template - auto difference(sets_t&&... sets) - { - static constexpr std::size_t dim = get_set_dim_v; - return std::apply( - [](auto&&... args) - { - return Subset(DifferenceOp(), start_end_function(), std::forward(args)...); - }, - std::make_tuple(detail::transform(std::forward(sets))...)); - } - - template - auto translate(set_t&& set, const stencil_t& stencil) - { - constexpr std::size_t dim = std::decay_t::dim; - return Subset(SelfOp(), - start_end_translate_function(xt::xtensor_fixed>(stencil)), - detail::transform(std::forward(set))); - } - - template - auto contraction(set_t&& set, int c) - { - constexpr std::size_t dim = std::decay_t::dim; - return Subset(SelfOp(), start_end_contraction_function(c), detail::transform(std::forward(set))); - } - - template - auto self(lca_t&& lca) - { - return Self>(std::forward(lca)); - } - - //----------------------------------------------------------------// - // Contract // - //----------------------------------------------------------------// - - template - auto contract_rec(const SubsetOrLCA& set, int width, const std::array& contract_directions) - { - static constexpr std::size_t dim = SubsetOrLCA::dim; - if constexpr (direction_index < dim) - { - using direction_t = xt::xtensor_fixed>; - - auto contracted_in_other_dirs = contract_rec(set, width, contract_directions); - direction_t dir; - dir.fill(0); - dir[direction_index] = contract_directions[direction_index] ? width : 0; - - return intersection(contracted_in_other_dirs, translate(set, dir), translate(set, -dir)); - } - else - { - if constexpr (IsLCA) - { - return self(set); - } - else - { - return set; - } - } - } - - /** - * @brief Contract a set in the specified directions. - * - * @tparam SubsetOrLCA The type of the set to contract. - * @param set The set or LevelCellArray to contract. - * @param width The contraction width. - * @param contract_directions An array indicating which directions to contract (true for contraction, false for no contraction). - * @return A new set that is contracted in the specified directions. - */ - template - auto contract(const SubsetOrLCA& set, std::size_t width, const std::array& contract_directions) - { - return contract_rec<0>(set, static_cast(width), contract_directions); - } - - /** - * @brief Contract a set in all directions. - * - * This function is a convenience wrapper that contracts the set in all dimensions. - * - * @tparam SubsetOrLCA The type of the set to contract. - * @param set The set or LevelCellArray to contract. - * @param width The contraction width. - * @return A new set that is contracted in all directions. - */ - template - auto contract(const SubsetOrLCA& set, std::size_t width) - { - std::array contract_directions; - std::fill(contract_directions.begin(), contract_directions.end(), true); - return contract(set, width, contract_directions); - } - - //--------------------------------------------------------------// - // Expand // - //--------------------------------------------------------------// - - template - auto expand_rec(const SubsetOrLCA& set, const std::array& expand_directions) - { - static constexpr std::size_t dim = SubsetOrLCA::dim; - using direction_t = xt::xtensor_fixed>; - - if constexpr (width == 0 && direction_index == dim) - { - if constexpr (IsLCA) - { - return self(set); - } - else - { - return set; - } - } - else if constexpr (width > 0) - { - auto expanded_layer = expand_rec(set, expand_directions); - direction_t dir; - dir.fill(0); - dir[direction_index] = expand_directions[direction_index] ? width : 0; - - return union_(expanded_layer, translate(set, dir), translate(set, -dir)); - } - else if constexpr (direction_index < dim) - { - auto expanded_in_other_dirs = expand_rec(set, expand_directions); - direction_t dir; - dir.fill(0); - dir[direction_index] = expand_directions[direction_index] ? width : 0; - - return union_(expanded_in_other_dirs, translate(set, dir), translate(set, -dir)); - } - } - - /** - * @brief Expand a set in the specified directions. - * - * @tparam SubsetOrLCA The type of the set to expand. - * @param set The set or LevelCellArray to expand. - * @param width The expansion width. - * @param expand_directions An array indicating which directions to expand (true for expansion, false for no expansion). - * @return A new set that is expanded in the specified directions. - */ - template - auto expand(const SubsetOrLCA& set, const std::array& expand_directions) - { - return expand_rec<0, width>(set, expand_directions); - } - - /** - * @brief Expand a set in all Cartesian directions (no diagonals!). - * - * This function is a convenience wrapper that expands the set in all dimensions. - * - * @tparam SubsetOrLCA The type of the set to expand. - * @param set The set or LevelCellArray to expand. - * @param width The expansion width. - * @return A new set that is expanded in all directions. - */ - template - auto expand(const SubsetOrLCA& set) - { - std::array expand_directions; - std::fill(expand_directions.begin(), expand_directions.end(), true); - return expand(set, expand_directions); - } -} +#include "box_view.hpp" +#include "contraction.hpp" +#include "expansion.hpp" +#include "lca_view.hpp" +#include "nary_set_operator.hpp" +#include "projection.hpp" +#include "set_base.hpp" +#include "translation.hpp" +#include "utils.hpp" diff --git a/include/samurai/subset/projection.hpp b/include/samurai/subset/projection.hpp new file mode 100644 index 000000000..26d32a108 --- /dev/null +++ b/include/samurai/subset/projection.hpp @@ -0,0 +1,193 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "../samurai_config.hpp" +#include "../static_algorithm.hpp" +#include "set_base.hpp" +#include "traversers/projection_traverser.hpp" +#include "utils.hpp" + +namespace samurai +{ + + namespace detail + { + template + struct ProjectionWork; + + template + struct ProjectionWork> + { + template + using child_traverser_t = typename Set::template traverser_t; + + template + using array_of_child_traverser_t = std::vector>; + + using Type = std::tuple...>; + }; + } // namespace detail + + template + class Projection; + + template + struct SetTraits> + { + static_assert(IsSet::value); + + template + using traverser_t = ProjectionTraverser>; + + struct Workspace + { + typename detail::ProjectionWork>::Type projection_workspace; + typename Set::Workspace child_workspace; + }; + + static constexpr std::size_t dim() + { + return Set::dim; + } + }; + + template + class Projection : public SetBase> + { + using Self = Projection; + using ChildTraverserArray = typename detail::ProjectionWork>::Type; + + public: + + SAMURAI_SET_TYPEDEFS + + Projection(const Set& set, const std::size_t level) + : m_set(set) + , m_level(level) + { + if (m_level < m_set.level()) + { + m_projectionType = ProjectionType::COARSEN; + m_shift = m_set.level() - m_level; + } + else + { + m_projectionType = ProjectionType::REFINE; + m_shift = m_level - m_set.level(); + } + } + + inline std::size_t level_impl() const + { + return m_level; + } + + inline bool exist_impl() const + { + return m_set.exist(); + } + + inline bool empty_impl() const + { + return m_set.empty(); + } + + template + inline void + init_workspace_impl(const std::size_t n_traversers, std::integral_constant d_ic, Workspace& workspace) const + { + const std::size_t my_work_size_per_traverser = (m_projectionType == ProjectionType::COARSEN and d != Base::dim - 1) + ? (1 << m_shift) + : 1; + const std::size_t my_work_size = n_traversers * my_work_size_per_traverser; + + auto& childTraversers = std::get(workspace.projection_workspace); + childTraversers.clear(); + childTraversers.reserve(my_work_size); + + m_set.init_workspace(my_work_size, d_ic, workspace.child_workspace); + } + + template + inline traverser_t + get_traverser_impl(const yz_index_t& index, std::integral_constant d_ic, Workspace& workspace) const + { + auto& childTraversers = std::get(workspace.projection_workspace); + if (m_projectionType == ProjectionType::COARSEN) + { + if constexpr (d != Base::dim - 1) + { + const auto childTraversers_begin = childTraversers.end(); + fill_traverser_array(index, d_ic, workspace); + return traverser_t(childTraversers_begin, childTraversers.end(), m_shift); + } + else + { + childTraversers.push_back(m_set.get_traverser(utils::pow2(index, m_shift), d_ic, workspace.child_workspace)); + return traverser_t(std::prev(childTraversers.end()), m_projectionType, m_shift); + } + } + else + { + childTraversers.push_back(m_set.get_traverser(utils::powMinus2(index, m_shift), d_ic, workspace.child_workspace)); + return traverser_t(std::prev(childTraversers.end()), m_projectionType, m_shift); + } + } + + template + inline traverser_t + get_traverser_unordered_impl(const yz_index_t& index, std::integral_constant d_ic, Workspace& workspace) const + { + auto& childTraversers = std::get(workspace.projection_workspace); + if (m_projectionType == ProjectionType::COARSEN) + { + if constexpr (d != Base::dim - 1) + { + const auto childTraversers_begin = childTraversers.end(); + fill_traverser_array(index, d_ic, workspace); + return traverser_t(childTraversers_begin, childTraversers.end(), m_shift); + } + else + { + childTraversers.push_back(m_set.get_traverser_unordered(utils::pow2(index, m_shift), d_ic, workspace.child_workspace)); + return traverser_t(std::prev(childTraversers.end()), m_projectionType, m_shift); + } + } + else + { + childTraversers.push_back(m_set.get_traverser_unordered(utils::powMinus2(index, m_shift), d_ic, workspace.child_workspace)); + return traverser_t(std::prev(childTraversers.end()), m_projectionType, m_shift); + } + } + + private: + + template + void fill_traverser_array(const yz_index_t& index, std::integral_constant d_ic, Workspace& workspace) const + { + auto& childTraversers = std::get(workspace.projection_workspace); + + const value_t ymin = index[d] << m_shift; + const value_t ybound = (index[d] + 1) << m_shift; + + yz_index_t projected_index(utils::pow2(index, m_shift)); + + for (projected_index[d] = ymin; projected_index[d] != ybound; ++projected_index[d]) + { + childTraversers.push_back(m_set.get_traverser_unordered(projected_index, d_ic, workspace.child_workspace)); + if (childTraversers.back().is_empty()) + { + childTraversers.pop_back(); + } + } + } + + Set m_set; + std::size_t m_level; + ProjectionType m_projectionType; + std::size_t m_shift; + }; + +} // namespace samurai diff --git a/include/samurai/subset/projection_loi.hpp b/include/samurai/subset/projection_loi.hpp new file mode 100644 index 000000000..0c4fedc54 --- /dev/null +++ b/include/samurai/subset/projection_loi.hpp @@ -0,0 +1,216 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "../samurai_config.hpp" +#include "../static_algorithm.hpp" +#include "set_base.hpp" +#include "traversers/last_dim_projection_loi_traverser.hpp" +#include "traversers/projection_loi_traverser.hpp" +#include "utils.hpp" + +#include + +namespace samurai +{ + + namespace detail + { + template + struct ProjectionLOIWork; + + template + struct ProjectionLOIWork> + { + template + using child_traverser_t = typename Set::template traverser_t; + + template + using work_t = ListOfIntervals::value_t>; + + using Type = std::tuple...>; + }; + } // namespace detail + + template + class ProjectionLOI; + + template + struct SetTraits> + { + static_assert(IsSet::value); + + template + using child_traverser_t = typename Set::template traverser_t; + + template + using traverser_t = std::conditional_t>, + ProjectionLOITraverser>>; + + struct Workspace + { + typename detail::ProjectionLOIWork>::Type projection_workspace; + typename Set::Workspace child_workspace; + }; + + static constexpr std::size_t dim() + { + return Set::dim; + } + }; + + template + class ProjectionLOI : public SetBase> + { + using Self = ProjectionLOI; + using ListOfIntervals = typename detail::ProjectionLOIWork>::Type; + + public: + + SAMURAI_SET_TYPEDEFS + + using ChildWorkspace = typename Set::Workspace; + + ProjectionLOI(const Set& set, const std::size_t level) + : m_set(set) + , m_level(level) + { + if (m_level < m_set.level()) + { + m_projectionType = ProjectionType::COARSEN; + m_shift = m_set.level() - m_level; + } + else + { + m_projectionType = ProjectionType::REFINE; + m_shift = m_level - m_set.level(); + } + } + + inline std::size_t level_impl() const + { + return m_level; + } + + inline bool exist_impl() const + { + return m_set.exist(); + } + + inline bool empty_impl() const + { + return m_set.empty(); + } + + template + inline void + init_workspace_impl(const std::size_t n_traversers, std::integral_constant d_ic, Workspace& workspace) const + { + assert(n_traversers == 1); + + m_set.init_workspace(n_traversers, d_ic, workspace.child_workspace); + } + + template + inline traverser_t + get_traverser_impl(const yz_index_t& index, std::integral_constant d_ic, Workspace& workspace) const + { + auto& listOfIntervals = std::get(workspace.projection_workspace); + listOfIntervals.clear(); + + if (m_projectionType == ProjectionType::COARSEN) + { + if constexpr (d != Base::dim - 1) + { + yz_index_t index_min(utils::pow2(index, m_shift)); + yz_index_t index_max(utils::sumAndPow2(index, 1, m_shift)); + + yz_index_t index_rec; + ChildWorkspace child_workspace; + fill_list_of_interval_rec(index_min, + index_max, + index_rec, + d_ic, + std::integral_constant{}, + child_workspace, + workspace); + + return traverser_t(m_set.get_traverser(utils::pow2(index, m_shift), d_ic, workspace.child_workspace), + listOfIntervals.cbegin(), + listOfIntervals.cend()); + } + else + { + return traverser_t(m_set.get_traverser(utils::pow2(index, m_shift), d_ic, workspace.child_workspace), + m_projectionType, + m_shift); + } + } + else + { + return traverser_t(m_set.get_traverser(utils::powMinus2(index, m_shift), d_ic, workspace.child_workspace), + m_projectionType, + m_shift); + } + } + + template + inline traverser_t + get_traverser_unordered_impl(const yz_index_t& index, std::integral_constant d_ic, Workspace& workspace) const + { + return get_traverser_impl(index, d_ic, workspace); + } + + private: + + template + inline void fill_list_of_interval_rec(const yz_index_t& index_min, + const yz_index_t& index_max, + yz_index_t& index, + std::integral_constant d_ic, + std::integral_constant dCur_ic, + ChildWorkspace& child_workspace, + Workspace& workspace_to_fill) const + { + using child_traverser_t = typename Set::template traverser_t; + using child_current_interval_t = typename child_traverser_t::current_interval_t; + + m_set.init_workspace(1, dCur_ic, child_workspace); + + for (child_traverser_t traverser = m_set.get_traverser(index, dCur_ic, child_workspace); !traverser.is_empty(); + traverser.next_interval()) + { + child_current_interval_t interval = traverser.current_interval(); + + if constexpr (dCur == d) + { + std::get(workspace_to_fill.projection_workspace).add_interval(interval >> m_shift); + } + else + { + const auto index_start = std::max(interval.start, index_min[dCur - 1]); + const auto index_end = std::min(interval.end, index_max[dCur - 1]); + + for (index[dCur - 1] = index_start; index[dCur - 1] < index_end; ++index[dCur - 1]) + { + fill_list_of_interval_rec(index_min, + index_max, + index, + d_ic, + std::integral_constant{}, + child_workspace, + workspace_to_fill); + } + } + } + } + + Set m_set; + std::size_t m_level; + ProjectionType m_projectionType; + std::size_t m_shift; + }; + +} // namespace samurai diff --git a/include/samurai/subset/projection_type.hpp b/include/samurai/subset/projection_type.hpp new file mode 100644 index 000000000..5489a0c0b --- /dev/null +++ b/include/samurai/subset/projection_type.hpp @@ -0,0 +1,13 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +namespace samurai +{ + enum class ProjectionType + { + COARSEN, + REFINE + }; +} // namespace samurai diff --git a/include/samurai/subset/set_base.hpp b/include/samurai/subset/set_base.hpp new file mode 100644 index 000000000..5c4d2e0e1 --- /dev/null +++ b/include/samurai/subset/set_base.hpp @@ -0,0 +1,222 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include + +#include + +#include "../samurai_config.hpp" +#include "traversers/set_traverser_base.hpp" + +namespace samurai +{ + //////////////////////////////////////////////////////////////////////// + //// Forward Declarations + //////////////////////////////////////////////////////////////////////// + + /** + * Traits used by `SetBase` + * it must define: + * 1. a template type traverser_t + * 2. a Workspace class + * 3. a constexpr dim() method + */ + template + struct SetTraits; + + template + class SetBase; + + template + class Projection; + + template + class ProjectionLOI; + + template + void apply(const SetBase& set, Func&& func); + + //////////////////////////////////////////////////////////////////////// + //// Class definition + //////////////////////////////////////////////////////////////////////// + + template + class SetBase + { + using DerivedTraits = SetTraits; + + public: + + template + using traverser_t = typename DerivedTraits::template traverser_t; + using Workspace = typename DerivedTraits::Workspace; + + using interval_t = typename traverser_t<0>::interval_t; + using value_t = typename interval_t::value_t; + + using yz_index_t = xt::xtensor_fixed>; + + using to_lca_t = LevelCellArray; + using to_lca_coord_t = typename to_lca_t::coords_t; + + using ProjectionMethod = std::conditional_t, Projection>; + + static constexpr std::size_t dim = DerivedTraits::dim(); + + const Derived& derived_cast() const + { + return static_cast(*this); + } + + Derived& derived_cast() + { + return static_cast(*this); + } + + inline std::size_t level() const + { + return derived_cast().level_impl(); + } + + inline bool exist() const + { + return derived_cast().exist_impl(); + } + + inline bool empty() const + { + return derived_cast().empty_impl(); + } + + template + inline void init_workspace(const std::size_t n_traversers, std::integral_constant d_ic, Workspace& workspace) const + { + derived_cast().init_workspace_impl(n_traversers, d_ic, workspace); + } + + //// Only works for increasing index + template + inline traverser_t get_traverser(const yz_index_t& index, std::integral_constant d_ic, Workspace& workspace) const + { + return derived_cast().get_traverser_impl(index, d_ic, workspace); + } + + //// Works for random indexes but might be slower + template + inline traverser_t + get_traverser_unordered(const yz_index_t& index, std::integral_constant d_ic, Workspace& workspace) const + { + return derived_cast().get_traverser_unordered_impl(index, d_ic, workspace); + } + + inline ProjectionMethod on(const std::size_t level); + + template + void operator()(Func&& func) const + { + apply(derived_cast(), std::forward(func)); + } + + template + void apply_op(ApplyOp&&... op) const + { + const std::size_t l = level(); + + auto func = [l, &op...](auto& interval, auto& index) + { + (op(l, interval, index), ...); + }; + apply(derived_cast(), func); + } + + to_lca_t to_lca() const + { + return to_lca_t(*this); + } + + to_lca_t to_lca(const to_lca_coord_t& origin_point, const double scaling_factor) const + { + return to_lca_t(*this, origin_point, scaling_factor); + } + + protected: + + inline bool empty_default_impl() const + { + yz_index_t index; + Workspace workspace; + return empty_default_impl_rec(index, std::integral_constant{}, workspace); + } + + template + bool empty_default_impl_rec(yz_index_t& index, std::integral_constant d_ic, Workspace& workspace) const + { + using current_interval_t = typename traverser_t::current_interval_t; + + init_workspace(1, d_ic, workspace); + + for (traverser_t traverser = get_traverser(index, d_ic, workspace); !traverser.is_empty(); traverser.next_interval()) + { + current_interval_t interval = traverser.current_interval(); + + if constexpr (d == 0) + { + return false; + } + else + { + for (index[d - 1] = interval.start; index[d - 1] != interval.end; ++index[d - 1]) + { + if (not empty_default_impl_rec(index, std::integral_constant{}, workspace)) + { + return false; + } + } + } + } + + return true; + } + }; + +#define SAMURAI_SET_TYPEDEFS \ + using Base = SetBase; \ + \ + template \ + using traverser_t = typename Base::template traverser_t; \ + \ + using Workspace = typename Base::Workspace; \ + \ + using interval_t = typename Base::interval_t; \ + using value_t = typename Base::value_t; \ + \ + using yz_index_t = typename Base::yz_index_t; + + template + struct IsSet : std::bool_constant, T>::value> + { + }; + + template + const Set& self(const SetBase& set) + { + return set.derived_cast(); + } + +} // namespace samurai + +#include "projection.hpp" +#include "projection_loi.hpp" + +namespace samurai +{ + + template + auto SetBase::on(const std::size_t level) -> ProjectionMethod + { + return ProjectionMethod(derived_cast(), level); + } + +} diff --git a/include/samurai/subset/start_end_fct.hpp b/include/samurai/subset/start_end_fct.hpp deleted file mode 100644 index 6c304ee28..000000000 --- a/include/samurai/subset/start_end_fct.hpp +++ /dev/null @@ -1,275 +0,0 @@ -// Copyright 2018-2025 the samurai's authors -// SPDX-License-Identifier: BSD-3-Clause - -#pragma once - -#include "utils.hpp" - -namespace samurai -{ - - inline auto default_function() - { - return [](auto, auto i, auto) - { - return i; - }; - } - - inline auto default_function_() - { - return [](auto level, auto i) - { - return std::make_pair(level, i); - }; - } - - template - struct start_end_function - { - auto& operator()(std::size_t level, std::size_t min_level, std::size_t max_level) - { - m_level = level; - m_min_level = min_level; - m_shift = static_cast(max_level) - static_cast(min_level); - return *this; - } - - template - inline auto start(const Func& f) const - { - auto new_f = [&, f](auto level, auto i, auto dec) - { - if constexpr (from_diff_op) - { - dec = (static_cast(level) > m_level) ? 1 : 0; - } - int value = (((i - dec) >> m_shift) << m_shift) + dec; - return f(m_level, value, dec); - }; - return new_f; - } - - template - inline auto end(const Func& f) const - { - auto new_f = [&, f](auto level, auto i, auto dec) - { - if constexpr (from_diff_op) - { - dec = (static_cast(level) > m_level) ? 0 : 1; - } - int value = (((i - dec) >> m_shift) + dec) << m_shift; - return f(m_level, value, dec); - }; - return new_f; - } - - template - inline auto goback(const Func& f) const - { - auto new_f = [&, f](auto level, auto i) - { - auto [prev_lev, v] = f(level, i); - - int min_shift = static_cast(m_min_level) - static_cast(prev_lev); - int max_shift = static_cast(m_level) - static_cast(m_min_level); - - if constexpr (end) - { - i = end_shift(end_shift(v, min_shift), max_shift); - } - else - { - i = start_shift(start_shift(v, min_shift), max_shift); - } - - return std::make_pair(m_level, i); - }; - return new_f; - } - - std::size_t m_level; - int m_shift; - std::size_t m_min_level; - }; - - template - struct start_end_translate_function - { - using container_t = xt::xtensor_fixed>; - - explicit start_end_translate_function(const container_t& t) - : m_level(0) - , m_min_level(0) - , m_max_level(0) - , m_t(t) - { - } - - auto& operator()(auto level, auto min_level, auto max_level) - { - m_level = level; - m_min_level = min_level; - m_max_level = max_level; - return *this; - } - - template - inline auto start(const Func& f) const - { - auto new_f = [&, f](auto level, auto i, auto dec) - { - int max2curr = static_cast(m_max_level) - static_cast(level); - int curr2min = static_cast(level) - static_cast(m_min_level); - int min2max = static_cast(m_max_level) - static_cast(m_min_level); - - if constexpr (from_diff_op) - { - dec = (static_cast(level) > m_level) ? 1 : 0; - } - int value = (((((i - dec) >> max2curr) + m_t[d - 1]) >> curr2min) + dec) << min2max; - - return f(m_level, value, dec); - }; - return new_f; - } - - template - inline auto end(const Func& f) const - { - auto new_f = [&, f](auto level, auto i, auto dec) - { - int max2curr = static_cast(m_max_level) - static_cast(level); - int curr2min = static_cast(level) - static_cast(m_min_level); - int min2max = static_cast(m_max_level) - static_cast(m_min_level); - - if constexpr (from_diff_op) - { - dec = (static_cast(level) > m_level) ? 0 : 1; - } - int value = (((((i - dec) >> max2curr) + m_t[d - 1]) >> curr2min) + dec) << min2max; - - return f(m_level, value, dec); - }; - return new_f; - } - - template - inline auto goback(const Func& f) const - { - auto new_f = [&, f](auto level, auto i) - { - auto [prev_lev, v] = f(level, i); - - auto min_shift = static_cast(m_min_level) - static_cast(prev_lev); - auto max_shift = static_cast(m_level) - static_cast(m_min_level); - - if constexpr (end) - { - i = end_shift(end_shift(v, min_shift), max_shift) - m_t[d - 1]; - } - else - { - i = start_shift(start_shift(v, min_shift), max_shift) - m_t[d - 1]; - } - - return std::make_pair(m_level, i); - }; - return new_f; - } - - std::size_t m_level; - std::size_t m_min_level; - std::size_t m_max_level; - xt::xtensor_fixed> m_t; - }; - - template - struct start_end_contraction_function - { - explicit start_end_contraction_function(int c) - : m_level(0) - , m_min_level(0) - , m_max_level(0) - , m_c(c) - { - } - - auto& operator()(auto level, auto min_level, auto max_level) - { - m_level = level; - m_min_level = min_level; - m_max_level = max_level; - return *this; - } - - template - inline auto start(const Func& f) const - { - auto new_f = [&, f](auto level, auto i, auto dec) - { - int max2curr = static_cast(m_max_level) - static_cast(level); - int curr2min = static_cast(level) - static_cast(m_min_level); - int min2max = static_cast(m_max_level) - static_cast(m_min_level); - - if constexpr (from_diff_op) - { - dec = (static_cast(level) > m_level) ? 1 : 0; - dec = 1; - } - int value = (((((i - dec) >> max2curr) + m_c) >> curr2min) + dec) << min2max; - return f(m_level, value, dec); - }; - return new_f; - } - - template - inline auto end(const Func& f) const - { - auto new_f = [&, f](auto level, auto i, auto dec) - { - int max2curr = static_cast(m_max_level) - static_cast(level); - int curr2min = static_cast(level) - static_cast(m_min_level); - int min2max = static_cast(m_max_level) - static_cast(m_min_level); - - if constexpr (from_diff_op) - { - dec = (static_cast(level) > m_level) ? 0 : 1; - } - int value = (((((i - dec) >> max2curr) - m_c) >> curr2min) + dec) << min2max; - return f(m_level, value, dec); - }; - return new_f; - } - - template - inline auto goback(const Func& f) const - { - auto new_f = [&, f](auto level, auto i) - { - auto [prev_lev, v] = f(level, i); - - auto min_shift = static_cast(m_min_level) - static_cast(prev_lev); - auto max_shift = static_cast(m_level) - static_cast(m_min_level); - - if constexpr (end) - { - i = end_shift(end_shift(v, min_shift), max_shift) + m_c; - } - else - { - i = start_shift(start_shift(v, min_shift), max_shift) + m_c; - } - - return std::make_pair(m_level, i); - }; - return new_f; - } - - std::size_t m_level; - std::size_t m_min_level; - std::size_t m_max_level; - int m_c; - }; -} diff --git a/include/samurai/subset/translation.hpp b/include/samurai/subset/translation.hpp new file mode 100644 index 000000000..4cc5fd3a3 --- /dev/null +++ b/include/samurai/subset/translation.hpp @@ -0,0 +1,117 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "set_base.hpp" +#include "traversers/translation_traverser.hpp" + +#include +using namespace xt::placeholders; // this makes `_` available + +namespace samurai +{ + + template + class Translation; + + template + struct SetTraits> + { + static_assert(IsSet::value); + + template + using traverser_t = TranslationTraverser>; + + struct Workspace + { + typename Set::Workspace child_workspace; + }; + + static constexpr std::size_t dim() + { + return Set::dim; + } + }; + + template + class Translation : public SetBase> + { + using Self = Translation; + + public: + + SAMURAI_SET_TYPEDEFS + + using translation_t = xt::xtensor_fixed>; + + template + Translation(const Set& set, const translation_expr_t& translation_expr) + : m_set(set) + , m_translation(translation_expr) + { + } + + inline std::size_t level_impl() const + { + return m_set.level(); + } + + inline bool exist_impl() const + { + return m_set.exist(); + } + + inline bool empty_impl() const + { + return m_set.empty(); + } + + template + inline void + init_workspace_impl(const std::size_t n_traversers, std::integral_constant d_ic, Workspace& workspace) const + { + m_set.init_workspace(n_traversers, d_ic, workspace.child_workspace); + } + + template + inline traverser_t + get_traverser_impl(const yz_index_t& index, std::integral_constant d_ic, Workspace& workspace) const + { + yz_index_t child_index; + + for (std::size_t i = 0; i != Base::dim - 1; ++i) + { + child_index[i] = index[i] - m_translation[i + 1]; + } + + return traverser_t(m_set.get_traverser_impl(child_index, d_ic, workspace.child_workspace), m_translation[d]); + } + + template + inline traverser_t + get_traverser_unordered_impl(const yz_index_t& index, std::integral_constant d_ic, Workspace& workspace) const + { + yz_index_t child_index; + + for (std::size_t i = 0; i != Base::dim - 1; ++i) + { + child_index[i] = index[i] - m_translation[i + 1]; + } + + return traverser_t(m_set.get_traverser_unordered_impl(child_index, d_ic, workspace.child_workspace), m_translation[d]); + } + + private: + + Set m_set; + translation_t m_translation; + }; + + template + auto translate(const Set& set, const typename Translation>::translation_t& translation) + { + return Translation(self(set), translation); + } + +} // namespace samurai diff --git a/include/samurai/subset/traversers/box_traverser.hpp b/include/samurai/subset/traversers/box_traverser.hpp new file mode 100644 index 000000000..f74e2c4aa --- /dev/null +++ b/include/samurai/subset/traversers/box_traverser.hpp @@ -0,0 +1,63 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "../../interval.hpp" +#include "set_traverser_base.hpp" +#include + +namespace samurai +{ + template + class Box; + + template + class BoxTraverser; + + template + struct SetTraverserTraits> + { + static_assert(std::same_as, B>); + + using interval_t = Interval; + using current_interval_t = const interval_t&; + }; + + template + class BoxTraverser : public SetTraverserBase> + { + using Self = BoxTraverser; + + public: + + SAMURAI_SET_TRAVERSER_TYPEDEFS + + BoxTraverser(const value_t& start, const value_t& end) + : m_current_interval{start, end} + , m_empty(start < end) + { + } + + inline bool is_empty_impl() const + { + return m_empty; + } + + inline void next_interval_impl() + { + m_empty = true; + } + + inline current_interval_t current_interval_impl() const + { + return m_current_interval; + } + + private: + + interval_t m_current_interval; + bool m_empty; + }; + +} // namespace samurai diff --git a/include/samurai/subset/traversers/contraction_traverser.hpp b/include/samurai/subset/traversers/contraction_traverser.hpp new file mode 100644 index 000000000..8d9876f8d --- /dev/null +++ b/include/samurai/subset/traversers/contraction_traverser.hpp @@ -0,0 +1,71 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "set_traverser_base.hpp" + +namespace samurai +{ + + template + class ContractionTraverser; + + template + struct SetTraverserTraits> + { + static_assert(IsSetTraverser::value); + + using interval_t = typename SetTraverser::interval_t; + using current_interval_t = interval_t; + }; + + template + class ContractionTraverser : public SetTraverserBase> + { + using Self = ContractionTraverser; + + public: + + SAMURAI_SET_TRAVERSER_TYPEDEFS + + ContractionTraverser(const SetTraverser& set_traverser, const value_t contraction) + : m_set_traverser(set_traverser) + , m_contraction(contraction) + { + assert(m_contraction >= 0); + advance_to_next_valid_interval(); + } + + inline bool is_empty_impl() const + { + return m_set_traverser.is_empty(); + } + + inline void next_interval_impl() + { + m_set_traverser.next_interval(); + advance_to_next_valid_interval(); + } + + inline current_interval_t current_interval_impl() const + { + return current_interval_t(m_set_traverser.current_interval().start + m_contraction, + m_set_traverser.current_interval().end - m_contraction); + } + + private: + + inline void advance_to_next_valid_interval() + { + while (!m_set_traverser.is_empty() && m_set_traverser.current_interval().size() <= size_t(2 * m_contraction)) + { + m_set_traverser.next_interval(); + } + } + + SetTraverser m_set_traverser; + value_t m_contraction; + }; + +} // namespace samurai diff --git a/include/samurai/subset/traversers/difference_id_traverser.hpp b/include/samurai/subset/traversers/difference_id_traverser.hpp new file mode 100644 index 000000000..af2f9b294 --- /dev/null +++ b/include/samurai/subset/traversers/difference_id_traverser.hpp @@ -0,0 +1,66 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "set_traverser_base.hpp" + +namespace samurai +{ + + template + class DifferenceIdTraverser; + + template + struct SetTraverserTraits> + { + static_assert(IsSetTraverser::value); + static_assert((IsSetTraverser::value and ...)); + + using interval_t = typename FirstSetTraverser::interval_t; + using current_interval_t = interval_t; + }; + + template + class DifferenceIdTraverser : public SetTraverserBase> + { + using Self = DifferenceIdTraverser; + + public: + + SAMURAI_SET_TRAVERSER_TYPEDEFS + + static constexpr std::size_t nIntervals = 1 + sizeof...(OtherSetTraversers); + + DifferenceIdTraverser(const std::array& shifts, + const FirstSetTraverser& set_traverser, + const OtherSetTraversers&...) + : m_set_traverser(set_traverser) + , m_shift(shifts[0]) + { + } + + DifferenceIdTraverser() = delete; + + inline bool is_empty_impl() const + { + return m_set_traverser.is_empty(); + } + + inline void next_interval_impl() + { + m_set_traverser.next_interval(); + } + + inline current_interval_t current_interval_impl() const + { + return current_interval_t{m_set_traverser.current_interval().start << m_shift, m_set_traverser.current_interval().end << m_shift}; + } + + private: + + FirstSetTraverser m_set_traverser; + std::size_t m_shift; + }; + +} // namespace samurai diff --git a/include/samurai/subset/traversers/difference_traverser.hpp b/include/samurai/subset/traversers/difference_traverser.hpp new file mode 100644 index 000000000..50d279d30 --- /dev/null +++ b/include/samurai/subset/traversers/difference_traverser.hpp @@ -0,0 +1,128 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "../../static_algorithm.hpp" +#include "set_traverser_base.hpp" + +namespace samurai +{ + + template + class DifferenceTraverser; + + template + struct SetTraverserTraits> + { + static_assert((IsSetTraverser::value and ...)); + + using FirstSetTraverser = std::tuple_element_t<0, std::tuple>; + using interval_t = typename FirstSetTraverser::interval_t; + using current_interval_t = const interval_t&; + }; + + template + class DifferenceTraverser : public SetTraverserBase> + { + using Self = DifferenceTraverser; + + public: + + SAMURAI_SET_TRAVERSER_TYPEDEFS + using Childrens = std::tuple; + + template + using IthChild = typename std::tuple_element::type; + + static constexpr std::size_t nIntervals = std::tuple_size_v; + + DifferenceTraverser(const std::array& shifts, const SetTraversers&... set_traversers) + : m_min_start(std::numeric_limits::min()) + , m_set_traversers(set_traversers...) + , m_shifts(shifts) + { + compute_current_interval(); + } + + inline bool is_empty_impl() const + { + return std::get<0>(m_set_traversers).is_empty(); + } + + inline void next_interval_impl() + { + advance_ref_interval(); + compute_current_interval(); + } + + inline current_interval_t current_interval_impl() const + { + return m_current_interval; + } + + private: + + inline void advance_ref_interval() + { + if (m_current_interval.end != std::get<0>(m_set_traversers).current_interval().end << m_shifts[0]) + { + // we have removed the beginning of the current interval. + // so ve remove [m_current_interval.start, m_current_interval.end) from std::get<0>(m_set_traversers).current_interval() + m_min_start = m_current_interval.end; + } + else + { + // all of the current interval has been removed. + // move to the next one. + std::get<0>(m_set_traversers).next_interval(); + } + } + + inline void compute_current_interval() + { + while (!std::get<0>(m_set_traversers).is_empty() && !try_to_compute_current_interval()) + { + advance_ref_interval(); + } + } + + inline bool try_to_compute_current_interval() + { + assert(!std::get<0>(m_set_traversers).is_empty()); + + m_current_interval.start = std::max(m_min_start, std::get<0>(m_set_traversers).current_interval().start << m_shifts[0]); + m_current_interval.end = std::get<0>(m_set_traversers).current_interval().end << m_shifts[0]; + + static_for<1, nIntervals>::apply( + [this](const auto i) + { + IthChild& set_traverser = std::get(m_set_traversers); + + while (!set_traverser.is_empty() && (set_traverser.current_interval().end << m_shifts[i]) < m_current_interval.start) + { + set_traverser.next_interval(); + } + + if (!set_traverser.is_empty() && (set_traverser.current_interval().start << m_shifts[i]) <= m_current_interval.start) + { + m_current_interval.start = set_traverser.current_interval().end << m_shifts[i]; + m_min_start = m_current_interval.start; + set_traverser.next_interval(); + } + if (!set_traverser.is_empty() && (set_traverser.current_interval().start << m_shifts[i]) <= m_current_interval.end) + { + m_current_interval.end = set_traverser.current_interval().start << m_shifts[i]; + } + }); + + return m_current_interval.is_valid(); + } + + interval_t m_current_interval; + value_t m_min_start; + Childrens m_set_traversers; + const std::array& m_shifts; + }; + +} // namespace samurai diff --git a/include/samurai/subset/traversers/expansion_traverser.hpp b/include/samurai/subset/traversers/expansion_traverser.hpp new file mode 100644 index 000000000..a88769d50 --- /dev/null +++ b/include/samurai/subset/traversers/expansion_traverser.hpp @@ -0,0 +1,99 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "set_traverser_base.hpp" + +namespace samurai +{ + + template + class ExpansionTraverser; + + template + struct SetTraverserTraits> + { + static_assert(IsSetTraverser::value); + + using interval_t = typename SetTraverserTraits::interval_t; + using current_interval_t = const interval_t&; + }; + + template + class ExpansionTraverser : public SetTraverserBase> + { + using Self = ExpansionTraverser; + + public: + + SAMURAI_SET_TRAVERSER_TYPEDEFS + + using SetTraverserIterator = typename std::vector::iterator; + + ExpansionTraverser(SetTraverserIterator begin_set_traverser, SetTraverserIterator end_set_traverser, const value_t expansion) + : m_set_traversers(begin_set_traverser, end_set_traverser) + , m_expansion(expansion) + { + assert(m_expansion > 0); + next_interval_impl(); + } + + inline bool is_empty_impl() const + { + return m_current_interval.start == std::numeric_limits::max(); + } + + inline void next_interval_impl() + { + m_current_interval.start = std::numeric_limits::max(); + + // We find the start of the interval, i.e. the smallest set_traverser.current_interval().start + for (const SetTraverser& set_traverser : m_set_traversers) + { + if (!set_traverser.is_empty() && (set_traverser.current_interval().start - m_expansion < m_current_interval.start)) + { + m_current_interval.start = set_traverser.current_interval().start - m_expansion; + m_current_interval.end = set_traverser.current_interval().end + m_expansion; + } + } + // Now we find the end of the interval, i.e. the largest set_traverser.current_interval().end + // such that set_traverser.current_interval().start - expansion < m_current_interval.end + bool is_done = false; + while (!is_done) + { + is_done = true; + // advance set traverses that are behind current interval + for (SetTraverser& set_traverser : m_set_traversers) + { + while (!set_traverser.is_empty() && set_traverser.current_interval().end + m_expansion <= m_current_interval.end) + { + set_traverser.next_interval(); + } + } + // try to find a new end + for (const SetTraverser& set_traverser : m_set_traversers) + { + // there is an overlap + if (!set_traverser.is_empty() && set_traverser.current_interval().start - m_expansion <= m_current_interval.end) + { + is_done = false; + m_current_interval.end = set_traverser.current_interval().end + m_expansion; + } + } + } + } + + inline current_interval_t current_interval_impl() const + { + return m_current_interval; + } + + private: + + std::span m_set_traversers; + value_t m_expansion; + interval_t m_current_interval; + }; + +} // namespace samurai diff --git a/include/samurai/subset/traversers/intersection_traverser.hpp b/include/samurai/subset/traversers/intersection_traverser.hpp new file mode 100644 index 000000000..86050bc23 --- /dev/null +++ b/include/samurai/subset/traversers/intersection_traverser.hpp @@ -0,0 +1,107 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "set_traverser_base.hpp" + +namespace samurai +{ + + template + class IntersectionTraverser; + + template + struct SetTraverserTraits> + { + static_assert((IsSetTraverser::value and ...)); + + using FirstSetTraverser = std::tuple_element_t<0, std::tuple>; + using interval_t = typename FirstSetTraverser::interval_t; + using current_interval_t = const interval_t&; + }; + + template + class IntersectionTraverser : public SetTraverserBase> + { + using Self = IntersectionTraverser; + + public: + + SAMURAI_SET_TRAVERSER_TYPEDEFS + using Childrens = std::tuple; + + template + using IthChild = typename std::tuple_element::type; + + static constexpr std::size_t nIntervals = std::tuple_size_v; + + IntersectionTraverser(const std::array& shifts, const SetTraversers&... set_traverser) + : m_set_traversers(set_traverser...) + , m_shifts(shifts) + { + next_interval_impl(); + } + + inline bool is_empty_impl() const + { + return !m_current_interval.is_valid(); + } + + inline void next_interval_impl() + { + m_current_interval.start = 0; + m_current_interval.end = 0; + + while (not_is_any_child_empty() && m_current_interval.start >= m_current_interval.end) + { + m_current_interval.start = std::numeric_limits::min(); + m_current_interval.end = std::numeric_limits::max(); + + enumerate_const_items( + m_set_traversers, + [this](const auto i, const auto& set_traverser) + { + if (!set_traverser.is_empty()) + { + m_current_interval.start = std::max(m_current_interval.start, + set_traverser.current_interval().start << m_shifts[i]); + m_current_interval.end = std::min(m_current_interval.end, set_traverser.current_interval().end << m_shifts[i]); + } + }); + + enumerate_items( + m_set_traversers, + [this](const auto i, auto& set_traverser) + { + if (!set_traverser.is_empty() && ((set_traverser.current_interval().end << m_shifts[i]) == m_current_interval.end)) + { + set_traverser.next_interval(); + } + }); + } + } + + inline current_interval_t current_interval_impl() const + { + return m_current_interval; + } + + private: + + inline bool not_is_any_child_empty() const + { + return std::apply( + [](const auto&... set_traversers) + { + return (!set_traversers.is_empty() and ...); + }, + m_set_traversers); + } + + interval_t m_current_interval; + Childrens m_set_traversers; + const std::array& m_shifts; + }; + +} // namespace samurai diff --git a/include/samurai/subset/traversers/last_dim_expansion_traverser.hpp b/include/samurai/subset/traversers/last_dim_expansion_traverser.hpp new file mode 100644 index 000000000..ab47787a7 --- /dev/null +++ b/include/samurai/subset/traversers/last_dim_expansion_traverser.hpp @@ -0,0 +1,76 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "set_traverser_base.hpp" + +namespace samurai +{ + + template + class LastDimExpansionTraverser; + + template + struct SetTraverserTraits> + { + static_assert(IsSetTraverser::value); + + using interval_t = typename SetTraverserTraits::interval_t; + using current_interval_t = const interval_t&; + }; + + template + class LastDimExpansionTraverser : public SetTraverserBase> + { + using Self = LastDimExpansionTraverser; + + public: + + SAMURAI_SET_TRAVERSER_TYPEDEFS + + LastDimExpansionTraverser(const SetTraverser& set_traverser, const value_t expansion) + : m_set_traverser(set_traverser) + , m_expansion(expansion) + { + assert(m_expansion > 0); + next_interval_impl(); + } + + inline bool is_empty_impl() const + { + return m_isEmpty; + } + + inline void next_interval_impl() + { + m_isEmpty = m_set_traverser.is_empty(); + + if (!m_isEmpty) + { + m_current_interval.start = m_set_traverser.current_interval().start - m_expansion; + m_current_interval.end = m_set_traverser.current_interval().end + m_expansion; + + m_set_traverser.next_interval(); + while (!m_set_traverser.is_empty() and m_set_traverser.current_interval().start - m_expansion <= m_current_interval.end) + { + m_current_interval.end = m_set_traverser.current_interval().end + m_expansion; + m_set_traverser.next_interval(); + } + } + } + + inline current_interval_t current_interval_impl() const + { + return m_current_interval; + } + + private: + + SetTraverser m_set_traverser; + value_t m_expansion; + interval_t m_current_interval; + bool m_isEmpty; + }; + +} // namespace samurai diff --git a/include/samurai/subset/traversers/last_dim_projection_loi_traverser.hpp b/include/samurai/subset/traversers/last_dim_projection_loi_traverser.hpp new file mode 100644 index 000000000..07847cc03 --- /dev/null +++ b/include/samurai/subset/traversers/last_dim_projection_loi_traverser.hpp @@ -0,0 +1,118 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "set_traverser_base.hpp" + +namespace samurai +{ + template + class LastDimProjectionLOITraverser; + + template + struct SetTraverserTraits> + { + static_assert(IsSetTraverser::value); + + using interval_t = typename SetTraverser::interval_t; + using current_interval_t = const interval_t&; + }; + + template + class LastDimProjectionLOITraverser : public SetTraverserBase> + { + using Self = LastDimProjectionLOITraverser; + + public: + + SAMURAI_SET_TRAVERSER_TYPEDEFS + + LastDimProjectionLOITraverser(const SetTraverser& set_traverser, const ProjectionType projectionType, const std::size_t shift) + : m_set_traverser(set_traverser) + , m_projectionType(projectionType) + , m_shift(shift) + , m_isEmpty(set_traverser.is_empty()) + { + if (m_projectionType == ProjectionType::COARSEN) + { + next_interval_coarsen(); + } + else if (!m_isEmpty) + { + m_current_interval.start = m_set_traverser.current_interval().start << shift; + m_current_interval.end = m_set_traverser.current_interval().end << shift; + } + } + + inline bool is_empty_impl() const + { + return m_isEmpty; + } + + inline void next_interval_impl() + { + if (m_projectionType == ProjectionType::COARSEN) + { + next_interval_coarsen(); + } + else + { + m_set_traverser.next_interval(); + m_isEmpty = m_set_traverser.is_empty(); + if (!m_isEmpty) + { + m_current_interval.start = m_set_traverser.current_interval().start << m_shift; + m_current_interval.end = m_set_traverser.current_interval().end << m_shift; + } + } + } + + inline current_interval_t current_interval_impl() const + { + return m_current_interval; + } + + private: + + inline void next_interval_coarsen() + { + if (!m_set_traverser.is_empty()) + { + m_current_interval.start = coarsen_start(m_set_traverser.current_interval()); + m_current_interval.end = coarsen_end(m_set_traverser.current_interval()); + + m_set_traverser.next_interval(); + + // when coarsening, two disjoint intervals may be merged. + // we need to check if the next_interval overlaps + for (; !m_set_traverser.is_empty() && coarsen_start(m_set_traverser.current_interval()) <= m_current_interval.end; + m_set_traverser.next_interval()) + { + m_current_interval.end = coarsen_end(m_set_traverser.current_interval()); + } + } + else + { + m_isEmpty = true; + } + } + + inline value_t coarsen_start(const interval_t& interval) const + { + return interval.start >> m_shift; + } + + inline value_t coarsen_end(const interval_t& interval) const + { + return ((interval.end - 1) >> m_shift) + 1; + } + + SetTraverser m_set_traverser; + ProjectionType m_projectionType; + std::size_t m_shift; + interval_t m_current_interval; + bool m_isEmpty; + }; + +} // namespace samurai diff --git a/include/samurai/subset/traversers/lca_traverser.hpp b/include/samurai/subset/traversers/lca_traverser.hpp new file mode 100644 index 000000000..00a1f3620 --- /dev/null +++ b/include/samurai/subset/traversers/lca_traverser.hpp @@ -0,0 +1,62 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "set_traverser_base.hpp" + +namespace samurai +{ + template + class LevelCellArray; + + template + class LCATraverser; + + template + struct SetTraverserTraits> + { + static_assert(std::same_as, LCA>); + + using interval_t = typename LCA::interval_t; + using current_interval_t = const interval_t&; + }; + + template + class LCATraverser : public SetTraverserBase> + { + using Self = LCATraverser; + + public: + + SAMURAI_SET_TRAVERSER_TYPEDEFS + using interval_iterator = typename std::vector::const_iterator; + + LCATraverser(const interval_iterator first, const interval_iterator end) + : m_first_interval(first) + , m_end_interval(end) + { + } + + inline bool is_empty_impl() const + { + return m_first_interval == m_end_interval; + } + + inline void next_interval_impl() + { + ++m_first_interval; + } + + inline current_interval_t current_interval_impl() const + { + return *m_first_interval; + } + + private: + + interval_iterator m_first_interval; + interval_iterator m_end_interval; + }; + +} // namespace samurai diff --git a/include/samurai/subset/traversers/projection_loi_traverser.hpp b/include/samurai/subset/traversers/projection_loi_traverser.hpp new file mode 100644 index 000000000..7244cad31 --- /dev/null +++ b/include/samurai/subset/traversers/projection_loi_traverser.hpp @@ -0,0 +1,88 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "../../list_of_intervals.hpp" +#include "set_traverser_base.hpp" + +namespace samurai +{ + template + class ProjectionLOITraverser; + + template + struct SetTraverserTraits> + { + static_assert(IsSetTraverser::value); + + using interval_t = typename SetTraverser::interval_t; + using current_interval_t = interval_t; + }; + + template + class ProjectionLOITraverser : public SetTraverserBase> + { + using Self = ProjectionLOITraverser; + using SetTraverserIterator = typename ListOfIntervals::const_iterator; + + public: + + SAMURAI_SET_TRAVERSER_TYPEDEFS + + ProjectionLOITraverser(SetTraverser set_traverser, [[maybe_unused]] const ProjectionType projectionType, const std::size_t shift) + : m_set_traverser(set_traverser) + , m_shift(shift) + , m_projectionType(ProjectionType::REFINE) + { + assert(projectionType == m_projectionType); + } + + ProjectionLOITraverser(SetTraverser set_traverser, SetTraverserIterator first_interval, SetTraverserIterator bound_interval) + : m_set_traverser(set_traverser) + , m_shift(0) + , m_first_interval(first_interval) + , m_bound_interval(bound_interval) + , m_projectionType(ProjectionType::COARSEN) + { + } + + inline bool is_empty_impl() const + { + if (m_projectionType == ProjectionType::COARSEN) + { + return m_first_interval == m_bound_interval; + } + else + { + return m_set_traverser.is_empty(); + } + } + + inline void next_interval_impl() + { + if (m_projectionType == ProjectionType::COARSEN) + { + ++m_first_interval; + } + else + { + m_set_traverser.next_interval(); + } + } + + inline current_interval_t current_interval_impl() const + { + return (m_projectionType == ProjectionType::COARSEN) ? *m_first_interval : m_set_traverser.current_interval() << m_shift; + } + + private: + + SetTraverser m_set_traverser; // only used when refining + std::size_t m_shift; // only used when refining + SetTraverserIterator m_first_interval; // only use when coarsening + SetTraverserIterator m_bound_interval; // only use when coarsening + ProjectionType m_projectionType; + }; + +} // namespace samurai diff --git a/include/samurai/subset/traversers/projection_traverser.hpp b/include/samurai/subset/traversers/projection_traverser.hpp new file mode 100644 index 000000000..7cdcb353a --- /dev/null +++ b/include/samurai/subset/traversers/projection_traverser.hpp @@ -0,0 +1,165 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include + +#include "../projection_type.hpp" +#include "set_traverser_base.hpp" + +namespace samurai +{ + template + class ProjectionTraverser; + + template + struct SetTraverserTraits> + { + static_assert(IsSetTraverser::value); + + using interval_t = typename SetTraverser::interval_t; + using current_interval_t = const interval_t&; + }; + + template + class ProjectionTraverser : public SetTraverserBase> + { + using Self = ProjectionTraverser; + using SetTraverserIterator = typename std::vector::iterator; + + public: + + SAMURAI_SET_TRAVERSER_TYPEDEFS + + ProjectionTraverser(SetTraverserIterator set_traverser, const ProjectionType projectionType, const std::size_t shift) + : m_set_traversers(set_traverser, set_traverser + 1) + , m_projectionType(projectionType) + , m_shift(shift) + , m_isEmpty(set_traverser->is_empty()) + { + if (!m_isEmpty) + { + if (m_projectionType == ProjectionType::COARSEN) + { + m_current_interval.start = coarsen_start(m_set_traversers.front().current_interval()); + m_current_interval.end = coarsen_end(m_set_traversers.front().current_interval()); + + m_set_traversers.front().next_interval(); + + // when coarsening, two disjoint intervals may be merged. + // we need to check if the next_interval overlaps + for (; !m_set_traversers.front().is_empty() + && coarsen_start(m_set_traversers.front().current_interval()) <= m_current_interval.end; + m_set_traversers.front().next_interval()) + { + m_current_interval.end = coarsen_end(m_set_traversers.front().current_interval()); + } + } + else + { + m_current_interval.start = m_set_traversers.front().current_interval().start << shift; + m_current_interval.end = m_set_traversers.front().current_interval().end << shift; + } + } + } + + /* + * This constructor only works for coarsening + */ + ProjectionTraverser(SetTraverserIterator begin_set_traversers, SetTraverserIterator end_set_traversers, const std::size_t shift) + : m_set_traversers(begin_set_traversers, end_set_traversers) + , m_projectionType(ProjectionType::COARSEN) + , m_shift(shift) + { + next_interval_coarsen(); + } + + inline bool is_empty_impl() const + { + return m_isEmpty; + } + + inline void next_interval_impl() + { + if (m_projectionType == ProjectionType::COARSEN) + { + next_interval_coarsen(); + } + else + { + m_set_traversers.front().next_interval(); + m_isEmpty = m_set_traversers.front().is_empty(); + if (!m_isEmpty) + { + m_current_interval.start = m_set_traversers.front().current_interval().start << m_shift; + m_current_interval.end = m_set_traversers.front().current_interval().end << m_shift; + } + } + } + + inline current_interval_t current_interval_impl() const + { + return m_current_interval; + } + + private: + + inline void next_interval_coarsen() + { + m_current_interval.start = std::numeric_limits::max(); + // We find the start of the interval, i.e. the smallest set_traverser.current_interval().start >> m_shift + for (const SetTraverser& set_traverser : m_set_traversers) + { + if (!set_traverser.is_empty() && (coarsen_start(set_traverser.current_interval()) < m_current_interval.start)) + { + m_current_interval.start = coarsen_start(set_traverser.current_interval()); + m_current_interval.end = coarsen_end(set_traverser.current_interval()); + } + } + // Now we find the end of the interval, i.e. the largest set_traverser.current_interval().end >> m_shift + // such that (set_traverser.current_interval().start >> m_shift) < m_current_interval.end + bool is_done = false; + while (!is_done) + { + is_done = true; + // advance set traverses that are behind current interval + for (SetTraverser& set_traverser : m_set_traversers) + { + while (!set_traverser.is_empty() && (coarsen_end(set_traverser.current_interval()) <= m_current_interval.end)) + { + set_traverser.next_interval(); + } + } + // try to find a new end + for (const SetTraverser& set_traverser : m_set_traversers) + { + // there is an overlap + if (!set_traverser.is_empty() && (coarsen_start(set_traverser.current_interval()) <= m_current_interval.end)) + { + is_done = false; + m_current_interval.end = coarsen_end(set_traverser.current_interval()); + } + } + } + m_isEmpty = (m_current_interval.start == std::numeric_limits::max()); + } + + inline value_t coarsen_start(const interval_t& interval) const + { + return interval.start >> m_shift; + } + + inline value_t coarsen_end(const interval_t& interval) const + { + return ((interval.end - 1) >> m_shift) + 1; + } + + std::span m_set_traversers; + ProjectionType m_projectionType; + std::size_t m_shift; + interval_t m_current_interval; + bool m_isEmpty; + }; + +} // namespace samurai diff --git a/include/samurai/subset/traversers/set_traverser_base.hpp b/include/samurai/subset/traversers/set_traverser_base.hpp new file mode 100644 index 000000000..6bcd63e96 --- /dev/null +++ b/include/samurai/subset/traversers/set_traverser_base.hpp @@ -0,0 +1,67 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include +#include + +namespace samurai +{ + /** + * Traits used by `SetTraverserBase` + * it must define: + * 1. an interval_t class + * 2. a current_interval_t class + */ + template + struct SetTraverserTraits; + + template + class SetTraverserBase + { + public: + + using interval_t = typename SetTraverserTraits::interval_t; + using current_interval_t = typename SetTraverserTraits::current_interval_t; + using value_t = typename interval_t::value_t; + + const Derived& derived_cast() const + { + return static_cast(*this); + } + + Derived& derived_cast() + { + return static_cast(*this); + } + + inline bool is_empty() const + { + return derived_cast().is_empty_impl(); + } + + inline void next_interval() + { + assert(!is_empty()); + derived_cast().next_interval_impl(); + } + + inline current_interval_t current_interval() const + { + return derived_cast().current_interval_impl(); + } + }; + + template + struct IsSetTraverser : std::bool_constant, T>::value> + { + }; + +#define SAMURAI_SET_TRAVERSER_TYPEDEFS \ + using Base = SetTraverserBase; \ + using interval_t = typename Base::interval_t; \ + using current_interval_t = typename Base::current_interval_t; \ + using value_t = typename Base::value_t; + +} // namespace samurai diff --git a/include/samurai/subset/traversers/translation_traverser.hpp b/include/samurai/subset/traversers/translation_traverser.hpp new file mode 100644 index 000000000..3f7b6cbb4 --- /dev/null +++ b/include/samurai/subset/traversers/translation_traverser.hpp @@ -0,0 +1,60 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "set_traverser_base.hpp" + +namespace samurai +{ + + template + class TranslationTraverser; + + template + struct SetTraverserTraits> + { + static_assert(IsSetTraverser::value); + + using interval_t = typename SetTraverserTraits::interval_t; + using current_interval_t = interval_t; + }; + + template + class TranslationTraverser : public SetTraverserBase> + { + using Self = TranslationTraverser; + + public: + + SAMURAI_SET_TRAVERSER_TYPEDEFS + + TranslationTraverser(const SetTraverser& set_traverser, const value_t& translation) + : m_set_traverser(set_traverser) + , m_translation(translation) + { + } + + inline bool is_empty_impl() const + { + return m_set_traverser.is_empty(); + } + + inline void next_interval_impl() + { + m_set_traverser.next_interval(); + } + + inline current_interval_t current_interval_impl() const + { + return current_interval_t{m_set_traverser.current_interval().start + m_translation, + m_set_traverser.current_interval().end + m_translation}; + } + + private: + + SetTraverser m_set_traverser; + value_t m_translation; + }; + +} // namespace samurai diff --git a/include/samurai/subset/traversers/union_traverser.hpp b/include/samurai/subset/traversers/union_traverser.hpp new file mode 100644 index 000000000..2b8d4fb98 --- /dev/null +++ b/include/samurai/subset/traversers/union_traverser.hpp @@ -0,0 +1,108 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "set_traverser_base.hpp" + +namespace samurai +{ + + template + class UnionTraverser; + + template + struct SetTraverserTraits> + { + static_assert((IsSetTraverser::value and ...)); + + using FirstSetTraverser = std::tuple_element_t<0, std::tuple>; + using interval_t = typename FirstSetTraverser::interval_t; + using current_interval_t = const interval_t&; + }; + + template + class UnionTraverser : public SetTraverserBase> + { + using Self = UnionTraverser; + + public: + + SAMURAI_SET_TRAVERSER_TYPEDEFS + using Childrens = std::tuple; + + template + using IthChild = typename std::tuple_element::type; + + static constexpr std::size_t nIntervals = std::tuple_size::value; + + UnionTraverser(const std::array& shifts, const SetTraversers&... set_traversers) + : m_set_traversers(set_traversers...) + , m_shifts(shifts) + { + next_interval_impl(); + } + + inline bool is_empty_impl() const + { + return m_current_interval.start == std::numeric_limits::max(); + } + + inline void next_interval_impl() + { + m_current_interval.start = std::numeric_limits::max(); + // We find the start of the interval, i.e. the smallest set_traverser.current_interval().start << m_shifts[i] + enumerate_const_items( + m_set_traversers, + [this](const auto i, const auto& set_traverser) + { + if (!set_traverser.is_empty() && ((set_traverser.current_interval().start << m_shifts[i]) < m_current_interval.start)) + { + m_current_interval.start = set_traverser.current_interval().start << m_shifts[i]; + m_current_interval.end = set_traverser.current_interval().end << m_shifts[i]; + } + }); + // Now we find the end of the interval, i.e. the largest set_traverser.current_interval().end << m_shifts[i] + // such that (set_traverser.current_interval().start << m_shifts[i]) < m_current_interval.end + bool is_done = false; + while (!is_done) + { + is_done = true; + // advance set traverses that are behind current interval + enumerate_items( + m_set_traversers, + [this](const auto i, auto& set_traverser) + { + while (!set_traverser.is_empty() && (set_traverser.current_interval().end << m_shifts[i]) <= m_current_interval.end) + { + set_traverser.next_interval(); + } + }); + // try to find a new end + enumerate_const_items( + m_set_traversers, + [&is_done, this](const auto i, const auto& set_traverser) + { + // there is an overlap + if (!set_traverser.is_empty() && (set_traverser.current_interval().start << m_shifts[i]) <= m_current_interval.end) + { + is_done = false; + m_current_interval.end = set_traverser.current_interval().end << m_shifts[i]; + } + }); + } + } + + inline current_interval_t current_interval_impl() const + { + return m_current_interval; + } + + private: + + interval_t m_current_interval; + Childrens m_set_traversers; + const std::array& m_shifts; + }; + +} // namespace samurai diff --git a/include/samurai/subset/utils.hpp b/include/samurai/subset/utils.hpp index 4c774c8ef..dd621c8ed 100644 --- a/include/samurai/subset/utils.hpp +++ b/include/samurai/subset/utils.hpp @@ -3,73 +3,133 @@ #pragma once -#include -#include -#include -#include +#include +#include -// namespace samurai::experimental namespace samurai { - template - static constexpr T sentinel = std::numeric_limits::max(); - - template - inline T end_shift(T value, T shift) + //////////////////////////////////////////////////////////////////////// + //// misc + //////////////////////////////////////////////////////////////////////// + template + const T& vmin(const T& a) { - return shift >= 0 ? value << shift : ((value - 1) >> -shift) + 1; + return a; } - template - inline T start_shift(T value, T shift) + template + const T& vmax(const T& a) { - return shift >= 0 ? value << shift : value >> -shift; + return a; } - constexpr auto compute_min(auto const& value, auto const&... args) + template + T vmin(const T& a, const T& b, const Ts&... others) { - if constexpr (sizeof...(args) == 0u) // Single argument case! + if constexpr (sizeof...(Ts) == 0u) { - return value; + return std::min(a, b); } - else // For the Ts... + else { - const auto min = compute_min(args...); - return value < min ? value : min; + return a < b ? vmin(a, others...) : vmin(b, others...); } } - constexpr auto compute_max(auto const& value, auto const&... args) + template + T vmax(const T& a, const T& b, const Ts&... others) { - if constexpr (sizeof...(args) == 0u) // Single argument case! + if constexpr (sizeof...(Ts) == 0u) { - return value; + return std::max(a, b); } - else // For the Ts... + else { - const auto max = compute_max(args...); - return value > max ? value : max; + return a > b ? vmax(a, others...) : vmax(b, others...); } } - template - struct get_interval_type + namespace utils { - using type = typename S1::interval_t; - }; - template - using get_interval_t = typename get_interval_type::type; + template + std::array pow2(const std::array& input, const std::size_t shift = 1) + { + std::array output; + for (std::size_t i = 0; i != N; ++i) + { + output[i] = input[i] << shift; + } - template - struct get_set_dim - { - static constexpr std::size_t value = S1::dim; - }; + return output; + } + + template + std::array sumAndPow2(const std::array& input, const T& value, const std::size_t shift = 1) + { + std::array output; + for (std::size_t i = 0; i != N; ++i) + { + output[i] = (input[i] + value) << shift; + } + + return output; + } + + template + std::array powMinus2(const std::array& input, const std::size_t shift = 1) + { + std::array output; + for (std::size_t i = 0; i != N; ++i) + { + output[i] = input[i] >> shift; + } + + return output; + } + + template + xt::xtensor_fixed> pow2(const xt::xtensor_fixed>& input, const std::size_t shift = 1) + { + xt::xtensor_fixed> output; + for (std::size_t i = 0; i != N; ++i) + { + output[i] = input[i] << shift; + } + + return output; + } + + template + xt::xtensor_fixed> + sumAndPow2(const xt::xtensor_fixed>& input, const T& value, const std::size_t shift = 1) + { + xt::xtensor_fixed> output; + for (std::size_t i = 0; i != N; ++i) + { + output[i] = (input[i] + value) << shift; + } - template - constexpr std::size_t get_set_dim_v = get_set_dim...>::value; + return output; + } + template + xt::xtensor_fixed> powMinus2(const xt::xtensor_fixed>& input, const std::size_t shift = 1) + { + xt::xtensor_fixed> output; + for (std::size_t i = 0; i != N; ++i) + { + output[i] = input[i] >> shift; + } + + return output; + } + + } + + //////////////////////////////////////////////////////////////////////// + //// intervals args + //////////////////////////////////////////////////////////////////////// template ::value_type::value_type> ForwardIt lower_bound_interval(ForwardIt begin, ForwardIt end, const T& value) { @@ -94,20 +154,40 @@ namespace samurai }); } - template - constexpr void zip_apply_impl(F&& f, Tuple1&& t1, Tuple2&& t2, std::index_sequence) + //////////////////////////////////////////////////////////////////////// + //// tuple iteration + //////////////////////////////////////////////////////////////////////// + namespace detail { - (f(std::get(std::forward(t1)), std::get(std::forward(t2))), ...); + template + Func enumerate_items(Tuple& tuple, Func func, std::index_sequence) + { + (func(Is, std::get(tuple)), ...); + return func; + } + + template + Func enumerate_const_items(const Tuple& tuple, Func func, std::index_sequence) + { + (func(Is, std::get(tuple)), ...); + return func; + } } - template - constexpr void zip_apply(F&& f, Tuple1&& t1, Tuple2&& t2) + template + Func enumerate_items(Tuple& tuple, Func func) { - constexpr std::size_t size1 = std::tuple_size_v>; - constexpr std::size_t size2 = std::tuple_size_v>; + constexpr std::size_t N = std::tuple_size_v>; - static_assert(size1 == size2, "Tuples must have the same size"); + return detail::enumerate_items(tuple, func, std::make_index_sequence{}); + } + + template + Func enumerate_const_items(const Tuple& tuple, Func func) + { + constexpr std::size_t N = std::tuple_size_v>; - zip_apply_impl(std::forward(f), std::forward(t1), std::forward(t2), std::make_index_sequence{}); + return detail::enumerate_const_items(tuple, func, std::make_index_sequence{}); } -} + +} // namespace samurai diff --git a/include/samurai/subset/visitor.hpp b/include/samurai/subset/visitor.hpp deleted file mode 100644 index 507c59f3f..000000000 --- a/include/samurai/subset/visitor.hpp +++ /dev/null @@ -1,379 +0,0 @@ -// Copyright 2018-2025 the samurai's authors -// SPDX-License-Identifier: BSD-3-Clause - -#pragma once - -#include -#include - -#include "utils.hpp" - -namespace samurai -{ - template - class IntervalListRange - { - public: - - using container_t = container_; - using value_t = typename container_t::value_type; - using iterator_t = typename container_t::const_iterator; - - IntervalListRange(const container_t& data, std::ptrdiff_t start, std::ptrdiff_t end) - : m_data(data) - , m_work(data) - , m_begin(data.cbegin() + start) - , m_end(data.cbegin() + end) - { - } - - IntervalListRange(const container_t& data, const container_t& w) - : m_data(data) - , m_work(w) - , m_begin(w.cbegin()) - , m_end(w.cend()) - { - } - - inline auto begin() const - { - return m_begin; - } - - inline auto end() const - { - return m_end; - } - - private: - - const container_t& m_data; - const container_t& m_work; - iterator_t m_begin; - iterator_t m_end; - }; - - template - class IntervalListVisitor - { - public: - - using container_t = container_; - using base_t = IntervalListRange; - using iterator_t = typename base_t::iterator_t; - using interval_t = typename base_t::value_t; - using value_t = typename interval_t::value_t; - - IntervalListVisitor(auto lca_level, auto level, auto max_level, const IntervalListRange& intervals) - : m_lca_level(static_cast(lca_level)) - , m_shift2dest(static_cast(max_level) - static_cast(level)) - , m_shift2ref(static_cast(max_level) - static_cast(lca_level)) - , m_intervals(intervals) - , m_first(intervals.begin()) - , m_last(intervals.end()) - , m_current(std::numeric_limits::min()) - , m_is_start(true) - { - } - - explicit IntervalListVisitor(IntervalListRange&& intervals) - : m_lca_level(std::numeric_limits::infinity()) - , m_shift2dest(std::numeric_limits::infinity()) - , m_shift2ref(std::numeric_limits::infinity()) - , m_intervals(std::move(intervals)) - , m_first(m_intervals.begin()) - , m_last(m_intervals.end()) - , m_current(sentinel) - , m_is_start(true) - { - } - - template - inline auto start(const auto& it, Func& start_fct) const - { - auto i = it->start << m_shift2ref; - return start_fct(m_lca_level, i, 0); - } - - template - inline auto end(const auto& it, Func& end_fct) const - { - auto i = it->end << m_shift2ref; - return end_fct(m_lca_level, i, 1); - } - - inline bool is_in(auto scan) const - { - // Recall that we check if scan is inside an interval defined as [start, - // end[. The end of the interval is not included. - // - // if the m_current value is the start of the interval which means m_is_start = - // true then if scan is lower than m_current, scan is not in the - // interval. - // - // if the m_current value is the end of the interval which means m_is_start = false - // then if scan is lower than m_current, scan is in the interval. - return m_current != sentinel && !((scan < m_current) ^ (!m_is_start)); - } - - inline bool is_empty() const - { - return m_current == sentinel; - } - - inline auto min() const - { - return m_current; - } - - inline auto shift() const - { - return m_shift2dest; - } - - template - inline void next_interval(StartEnd& start_and_stop) - { - auto& [start_fct, end_fct] = start_and_stop; // cppcheck-suppress variableScope - - auto i_start = start(m_first, start_fct); - auto i_end = end(m_first, end_fct); - while (m_first + 1 != m_last && i_end >= start(m_first + 1, start_fct)) - { - ++m_first; - i_end = end(m_first, end_fct); - } - m_current_interval = {i_start, i_end}; - - if (m_current_interval.is_valid()) - { - m_current = m_current_interval.start; - } - else - { - m_current = sentinel; - } - } - - template - inline void next(auto scan, StartEnd& start_and_stop) - { - if (m_current == std::numeric_limits::min()) - { - next_interval(start_and_stop); - return; - } - - if (m_current == scan) - { - if (m_is_start) - { - m_current = m_current_interval.end; - } - else - { - ++m_first; - - if (m_first == m_last) - { - m_current = sentinel; - return; - } - next_interval(start_and_stop); - } - m_is_start = !m_is_start; - } - } - - private: - - int m_lca_level; - int m_shift2dest; - int m_shift2ref; - IntervalListRange m_intervals; - iterator_t m_first; - iterator_t m_last; - value_t m_current; - interval_t m_current_interval; - bool m_is_start; - }; - - template - class SetTraverser - { - public: - - static constexpr std::size_t dim = get_set_dim_v; - using set_type = std::tuple; - using interval_t = get_interval_t; - - SetTraverser(int shift, const Operator& op, S&&... s) - : m_shift(shift) - , m_operator(op) - , m_s(std::forward(s)...) - { - } - - inline auto shift() const - { - return m_shift; - } - - inline bool is_in(auto scan) const - { - return std::apply( - [this, scan](auto&&... args) - { - return m_operator.is_in(scan, args...); - }, - m_s); - } - - inline bool is_empty() const - { - return std::apply( - [this](auto&&... args) - { - return m_operator.is_empty(args...); - }, - m_s); - } - - inline auto min() const - { - return std::apply( - [](auto&&... args) - { - return compute_min(args.min()...); - }, - m_s); - } - - template - void next(auto scan, StartEnd&& start_and_stop) - { - zip_apply( - [scan](auto& arg, auto& start_end_fct) - { - arg.next(scan, start_end_fct); - }, - m_s, - std::forward(start_and_stop)); - } - - private: - - int m_shift; - Operator m_operator; - set_type m_s; - }; - - struct IntersectionOp - { - bool is_in(auto scan, const auto&... args) const - { - return (args.is_in(scan) && ...); - } - - bool is_empty(const auto&... args) const - { - return (args.is_empty() || ...); - } - - bool exist(const auto&... args) const - { - return (args.exist() && ...); - } - }; - - struct UnionOp - { - bool is_in(auto scan, const auto&... args) const - { - return (args.is_in(scan) || ...); - } - - bool is_empty(const auto&... args) const - { - return (args.is_empty() && ...); - } - - bool exist(const auto&... args) const - { - return (args.exist() || ...); - } - }; - - struct DifferenceOp - { - bool is_in(auto scan, const auto& arg, const auto&... args) const - { - return arg.is_in(scan) && !(args.is_in(scan) || ...); - } - - bool is_empty(const auto& arg, const auto&...) const - { - return arg.is_empty(); - } - - bool exist(const auto& arg, const auto&...) const - { - return arg.exist(); - } - }; - - struct Difference2Op - { - bool is_in(auto scan, const auto& arg, const auto&...) const - { - return arg.is_in(scan); - } - - bool is_empty(const auto& arg, const auto&...) const - { - return arg.is_empty(); - } - - bool exist(const auto& arg, const auto&...) const - { - return arg.exist(); - } - }; - - template - auto get_operator(const operator_t& op) - { - return op; - } - - template - auto get_operator(const DifferenceOp& op) - { - if constexpr (d == 1) - { - return op; - } - else - { - return Difference2Op(); - } - } - - struct SelfOp - { - bool is_in(auto scan, const auto& arg) const - { - return arg.is_in(scan); - } - - bool is_empty(const auto& arg) const - { - return arg.is_empty(); - } - - bool exist(const auto& arg) const - { - return arg.exist(); - } - }; -} diff --git a/tests/test_subset.cpp b/tests/test_subset.cpp index 1528938d4..2239eb87b 100644 --- a/tests/test_subset.cpp +++ b/tests/test_subset.cpp @@ -14,6 +14,9 @@ #include #include +#include +#include + namespace samurai { TEST(subset, lower_bound) @@ -144,9 +147,9 @@ namespace samurai TEST(subset, compute_min) { - EXPECT_EQ(1, compute_min(3, 4, 1, 4)); - EXPECT_EQ(0, compute_min(0, 0, 0, 0)); - EXPECT_EQ(-1, compute_min(-1, -1, -1, -1)); + EXPECT_EQ(1, vmin(3, 4, 1, 4)); + EXPECT_EQ(0, vmin(0, 0, 0, 0)); + EXPECT_EQ(-1, vmin(-1, -1, -1, -1)); } TEST(subset, check_dim) @@ -224,8 +227,9 @@ namespace samurai EXPECT_EQ(interval_t(0, 20), i); }); - EXPECT_EQ(set.on(5).level(), 5); - apply(set, + auto set2 = set.on(5); + EXPECT_EQ(set2.level(), 5); + apply(set2, [](auto& i, auto) { EXPECT_EQ(interval_t(0, 40), i); @@ -511,6 +515,134 @@ namespace samurai } } + TEST(subset, expand) + { + using interval_t = typename LevelCellArray<2>::interval_t; + using expected_t = std::vector>; + + LevelCellArray<2> ca; + + ca.add_interval_back({0, 1}, {0}); + + { + const auto translated_ca = translate(ca, {3 + 1, 0}); + const auto joined_cas = union_(ca, translated_ca); + + const auto set = expand(joined_cas, 3); + + expected_t expected{ + {-3, {-3, 8}}, + {-2, {-3, 8}}, + {-1, {-3, 8}}, + {0, {-3, 8}}, + {1, {-3, 8}}, + {2, {-3, 8}}, + {3, {-3, 8}} + }; + + bool is_set_empty = true; + std::size_t ie = 0; + set( + [&expected, &is_set_empty, &ie](const auto& x_interval, const auto& yz) + { + is_set_empty = false; + EXPECT_EQ(expected[ie++], std::make_pair(yz[0], x_interval)); + }); + EXPECT_EQ(ie, expected.size()); + EXPECT_FALSE(is_set_empty); + } + + { + const auto translated_ca = translate(ca, {0, 3 + 1}); + const auto joined_cas = union_(ca, translated_ca); + + const auto set = expand(joined_cas, 3); + + expected_t expected{ + {-3, {-3, 4}}, + {-2, {-3, 4}}, + {-1, {-3, 4}}, + {0, {-3, 4}}, + {1, {-3, 4}}, + {2, {-3, 4}}, + {3, {-3, 4}}, + {4, {-3, 4}}, + {5, {-3, 4}}, + {6, {-3, 4}}, + {7, {-3, 4}} + }; + + bool is_set_empty = true; + std::size_t ie = 0; + set( + [&expected, &is_set_empty, &ie](const auto& x_interval, const auto& yz) + { + is_set_empty = false; + EXPECT_EQ(expected[ie++], std::make_pair(yz[0], x_interval)); + }); + EXPECT_EQ(ie, expected.size()); + EXPECT_FALSE(is_set_empty); + } + { + const auto translated_ca = translate(ca, {3 + 1, 3 + 1}); + const auto joined_cas = union_(ca, translated_ca); + + const auto set = expand(joined_cas, 3); + + expected_t expected{ + {-3, {-3, 4}}, + {-2, {-3, 4}}, + {-1, {-3, 4}}, + {0, {-3, 4}}, + {1, {-3, 8}}, + {2, {-3, 8}}, + {3, {-3, 8}}, + {4, {1, 8} }, + {5, {1, 8} }, + {6, {1, 8} }, + {7, {1, 8} } + }; + + bool is_set_empty = true; + std::size_t ie = 0; + set( + [&expected, &is_set_empty, &ie](const auto& x_interval, const auto& yz) + { + is_set_empty = false; + EXPECT_EQ(expected[ie++], std::make_pair(yz[0], x_interval)); + }); + EXPECT_EQ(ie, expected.size()); + EXPECT_FALSE(is_set_empty); + + const auto lca_joined_cas = joined_cas.to_lca(); + const auto lca_set = set.to_lca(); + } + } + + TEST(subset, contract) + { + LevelCellArray<2> ca; + + ca.add_interval_back({0, 1}, {0}); + + { + const auto translated_ca = translate(ca, {3 + 1, 0}); + const auto joined_cas = union_(ca, translated_ca); + + const auto set = contract(joined_cas, 1); + + bool is_set_empty = true; + set( + [&is_set_empty](const auto& x_interval, const auto& yz) + { + fmt::print("x_interval = {} -- yz = {}", x_interval, yz[0]); + is_set_empty = false; + }); + EXPECT_TRUE(is_set_empty); + //~ EXPECT_TRUE(set.empty()); + } + } + TEST(subset, translate) { CellList<1> cl; @@ -1172,6 +1304,13 @@ namespace samurai ca = {cl, true}; + //~ std::cout << "===========================================" << std::endl; + //~ std::cout << ca[4] << std::endl; + //~ std::cout << "===========================================" << std::endl; + //~ std::cout << self(ca[4]).on(3).to_lca() << std::endl; + //~ std::cout << "===========================================" << std::endl; + //~ std::exit(0); + // Test self-similarity at different scales bool found = false; apply(intersection(ca[3], self(ca[4]).on(3)),