From d6fcdc9043abcc4015bd0ba12f9fd755c9c0ca79 Mon Sep 17 00:00:00 2001
From: Loic Gouarin <loic.gouarin@gmail.com>
Date: Fri, 28 Jul 2023 11:41:25 +0200
Subject: [PATCH] feat: add cell_in in functionBC (#121)

- modify template parameter of a Cell
- add the cell inside the domain in the method which computes the values
in FunctionBC
---
 demos/FiniteVolume/level_set_from_scratch.cpp |   2 +-
 demos/FiniteVolume/stokes_2d.cpp              |  12 +-
 demos/highorder/main.cpp                      |   2 +-
 demos/multigrid/test_cases.hpp                |  27 ++--
 demos/p4est/simple_2d.cpp                     |   4 +-
 demos/pablo/bubble_2d.cpp                     |   2 +-
 .../AMR_1D_Burgers/step_4/update_field.hpp    |   2 +-
 .../AMR_1D_Burgers/step_4/update_mesh.hpp     |   2 +-
 demos/tutorial/graduation_case_3.cpp          |   2 +-
 include/samurai/algorithm.hpp                 |  29 ++---
 include/samurai/algorithm/update.hpp          |   4 +-
 include/samurai/bc.hpp                        | 119 +++++++++---------
 include/samurai/cell.hpp                      |  56 +++++----
 include/samurai/cell_array.hpp                |  11 +-
 include/samurai/field.hpp                     |  30 ++---
 include/samurai/field_expression.hpp          |   4 +-
 include/samurai/graduation.hpp                |   2 +-
 include/samurai/interface.hpp                 |   6 +-
 include/samurai/io/from_geometry.hpp          |   4 +-
 include/samurai/level_cell_array.hpp          |  12 +-
 include/samurai/level_cell_list.hpp           |   4 +-
 include/samurai/mesh.hpp                      |   9 +-
 include/samurai/numeric/gauss_legendre.hpp    |   8 +-
 .../samurai/petsc/fv/FV_scheme_assembly.hpp   |  68 +++++-----
 include/samurai/stencil.hpp                   |  23 ++--
 test/test_cell.cpp                            |   7 +-
 26 files changed, 225 insertions(+), 226 deletions(-)

diff --git a/demos/FiniteVolume/level_set_from_scratch.cpp b/demos/FiniteVolume/level_set_from_scratch.cpp
index 6d2d70098..f6b99c883 100644
--- a/demos/FiniteVolume/level_set_from_scratch.cpp
+++ b/demos/FiniteVolume/level_set_from_scratch.cpp
@@ -334,7 +334,7 @@ bool update_mesh(Field& f, Field_u& u, Tag& tag)
     samurai::for_each_interval(mesh[SimpleID::cells],
                                [&](std::size_t level, const auto& interval, const auto& index_yz)
                                {
-                                   auto itag = static_cast<std::size_t>(interval.start + interval.index);
+                                   auto itag = interval.start + interval.index;
                                    for (int i = interval.start; i < interval.end; ++i)
                                    {
                                        if (tag[itag] & static_cast<int>(samurai::CellFlag::refine))
diff --git a/demos/FiniteVolume/stokes_2d.cpp b/demos/FiniteVolume/stokes_2d.cpp
index ea9373427..6f21b0d05 100644
--- a/demos/FiniteVolume/stokes_2d.cpp
+++ b/demos/FiniteVolume/stokes_2d.cpp
@@ -230,7 +230,7 @@ int main(int argc, char* argv[])
 
         // Boundary conditions
         samurai::make_bc<samurai::Dirichlet>(velocity,
-                                             [](const auto& coord)
+                                             [](const auto&, const auto& coord)
                                              {
                                                  const auto& x = coord[0];
                                                  const auto& y = coord[1];
@@ -240,7 +240,7 @@ int main(int argc, char* argv[])
                                              });
 
         samurai::make_bc<samurai::Neumann>(pressure,
-                                           [](const auto& coord)
+                                           [](const auto&, const auto& coord)
                                            {
                                                const auto& x = coord[0];
                                                const auto& y = coord[1];
@@ -393,12 +393,12 @@ int main(int argc, char* argv[])
 
         // Boundary conditions
         samurai::make_bc<samurai::Dirichlet>(velocity_np1,
-                                             [&](const auto& coord)
+                                             [&](const auto&, const auto& coord)
                                              {
                                                  return exact_velocity(0, coord);
                                              });
         samurai::make_bc<samurai::Neumann>(pressure_np1,
-                                           [&](const auto& coord)
+                                           [&](const auto&, const auto& coord)
                                            {
                                                return exact_normal_grad_pressure(0, coord);
                                            });
@@ -480,13 +480,13 @@ int main(int argc, char* argv[])
             // Boundary conditions
             velocity_np1.get_bc().clear();
             samurai::make_bc<samurai::Dirichlet>(velocity_np1,
-                                                 [&](const auto& coord)
+                                                 [&](const auto&, const auto& coord)
                                                  {
                                                      return exact_velocity(t_np1, coord);
                                                  });
             pressure_np1.get_bc().clear();
             samurai::make_bc<samurai::Neumann>(pressure_np1,
-                                               [&](const auto& coord)
+                                               [&](const auto&, const auto& coord)
                                                {
                                                    return exact_normal_grad_pressure(t_np1, coord);
                                                });
diff --git a/demos/highorder/main.cpp b/demos/highorder/main.cpp
index 5ef594863..1e1c960bd 100644
--- a/demos/highorder/main.cpp
+++ b/demos/highorder/main.cpp
@@ -267,7 +267,7 @@ int main(int argc, char* argv[])
 
         auto u = samurai::make_field<double, 1>("u", mesh);
         samurai::make_bc<samurai::Dirichlet>(u,
-                                             [](const auto& coord)
+                                             [](const auto&, const auto& coord)
                                              {
                                                  const auto& x = coord[0];
                                                  const auto& y = coord[1];
diff --git a/demos/multigrid/test_cases.hpp b/demos/multigrid/test_cases.hpp
index 1c4c56545..f91afce02 100644
--- a/demos/multigrid/test_cases.hpp
+++ b/demos/multigrid/test_cases.hpp
@@ -8,11 +8,12 @@ class TestCase
   public:
 
     static constexpr std::size_t dim = Field::dim;
-    using coords                     = xt::xtensor_fixed<double, xt::xshape<dim>>;
+    using cell_t                     = typename Field::cell_t;
+    using coords_t                   = typename cell_t::coords_t;
 
-    using field_value_t    = typename samurai::detail::return_type<double, Field::size>::type;
-    using field_function_t = typename samurai::FunctionBc<dim, double, Field::size>::function_t;
-    using boundary_cond_t  = field_function_t;
+    using boundary_cond_t  = typename samurai::FunctionBc<Field>::function_t;
+    using field_value_t    = typename samurai::FunctionBc<Field>::value_t;
+    using field_function_t = std::function<field_value_t(const coords_t&)>;
 
     virtual bool solution_is_known()
     {
@@ -41,9 +42,12 @@ class TestCase
     {
         if (solution_is_known())
         {
-            return solution();
+            return [&](const auto&, const auto& coords)
+            {
+                return solution()(coords);
+            };
         }
-        return [](const coords&)
+        return [](const cell_t&, const coords_t&)
         {
             if constexpr (Field::size == 1)
             {
@@ -80,7 +84,8 @@ class PolynomialTestCase : public TestCase<Field>
     using field_value_t              = typename TestCase<Field>::field_value_t;
     using field_function_t           = typename TestCase<Field>::field_function_t;
     using boundary_cond_t            = typename TestCase<Field>::boundary_cond_t;
-    using coords                     = typename TestCase<Field>::coords;
+    using coords_t                   = typename TestCase<Field>::coords_t;
+    using cell_t                     = typename TestCase<Field>::cell_t;
 
     bool solution_is_known() override
     {
@@ -145,7 +150,7 @@ class PolynomialTestCase : public TestCase<Field>
 
     boundary_cond_t dirichlet() override
     {
-        return [](const coords&)
+        return [](const cell_t&, const coords_t&)
         {
             if constexpr (Field::size == 1)
             {
@@ -164,7 +169,7 @@ class PolynomialTestCase : public TestCase<Field>
     {
         if constexpr (dim == 1)
         {
-            return [](const coords& coord)
+            return [](const auto&, const coords_t& coord)
             {
                 const auto& x = coord[0];
                 if (x == 0 || x == 1)
@@ -180,11 +185,11 @@ class PolynomialTestCase : public TestCase<Field>
         }
         else if constexpr (dim == 2)
         {
-            return [](const auto& coord)
+            return [](const auto&, const auto& coord)
             {
                 const auto& x = coord[0];
                 const auto& y = coord[1];
-                double value;
+                double value  = 0.;
                 if (x == 0 || x == 1)
                 {
                     value = y * (y - 1);
diff --git a/demos/p4est/simple_2d.cpp b/demos/p4est/simple_2d.cpp
index 9c33610ac..9a40380b5 100644
--- a/demos/p4est/simple_2d.cpp
+++ b/demos/p4est/simple_2d.cpp
@@ -65,7 +65,7 @@ void refine_1(mesh_t& mesh, std::size_t max_level)
         samurai::for_each_interval(mesh,
                                    [&](std::size_t level, const auto& interval, const auto& index_yz)
                                    {
-                                       auto itag = static_cast<std::size_t>(interval.start + interval.index);
+                                       auto itag = interval.start + interval.index;
                                        for (coord_index_t i = interval.start; i < interval.end; ++i)
                                        {
                                            if (cell_tag[itag])
@@ -165,7 +165,7 @@ void refine_2(mesh_t& mesh, std::size_t max_level)
         samurai::for_each_interval(mesh[mesh_id_t::cells],
                                    [&](std::size_t level, const auto& interval, const auto& index_yz)
                                    {
-                                       auto itag = static_cast<std::size_t>(interval.start + interval.index);
+                                       auto itag = interval.start + interval.index;
                                        for (coord_index_t i = interval.start; i < interval.end; ++i)
                                        {
                                            if (cell_tag[itag])
diff --git a/demos/pablo/bubble_2d.cpp b/demos/pablo/bubble_2d.cpp
index 50964e2cb..3b08d06f9 100644
--- a/demos/pablo/bubble_2d.cpp
+++ b/demos/pablo/bubble_2d.cpp
@@ -85,7 +85,7 @@ void update_mesh(Mesh& mesh,
     samurai::for_each_interval(mesh,
                                [&](std::size_t level, const auto& interval, const auto& index_yz)
                                {
-                                   auto itag = static_cast<std::size_t>(interval.start + interval.index);
+                                   auto itag = interval.start + interval.index;
                                    for (int i = interval.start; i < interval.end; ++i)
                                    {
                                        if (tag[itag] & static_cast<int>(samurai::CellFlag::refine))
diff --git a/demos/tutorial/AMR_1D_Burgers/step_4/update_field.hpp b/demos/tutorial/AMR_1D_Burgers/step_4/update_field.hpp
index d1fafdfd1..696aef945 100644
--- a/demos/tutorial/AMR_1D_Burgers/step_4/update_field.hpp
+++ b/demos/tutorial/AMR_1D_Burgers/step_4/update_field.hpp
@@ -48,7 +48,7 @@ void update_field(Field& f, const Tag& tag, Mesh& new_mesh)
     samurai::for_each_interval(mesh[mesh_id_t::cells],
                                [&](std::size_t level, const auto& interval, const auto&)
                                {
-                                   auto itag = static_cast<std::size_t>(interval.start + interval.index);
+                                   auto itag = interval.start + interval.index;
                                    for (coord_index_t i = interval.start; i < interval.end; ++i)
                                    {
                                        if (tag[itag] & static_cast<int>(samurai::CellFlag::refine))
diff --git a/demos/tutorial/AMR_1D_Burgers/step_4/update_mesh.hpp b/demos/tutorial/AMR_1D_Burgers/step_4/update_mesh.hpp
index 558595e3d..782161a62 100644
--- a/demos/tutorial/AMR_1D_Burgers/step_4/update_mesh.hpp
+++ b/demos/tutorial/AMR_1D_Burgers/step_4/update_mesh.hpp
@@ -41,7 +41,7 @@ bool update_mesh(Field& f, const Tag& tag)
     samurai::for_each_interval(mesh[mesh_id_t::cells],
                                [&](std::size_t level, const auto& interval, const auto&)
                                {
-                                   auto itag = static_cast<std::size_t>(interval.start + interval.index);
+                                   auto itag = interval.start + interval.index;
                                    for (int i = interval.start; i < interval.end; ++i)
                                    {
                                        if (tag[itag] & static_cast<int>(samurai::CellFlag::refine))
diff --git a/demos/tutorial/graduation_case_3.cpp b/demos/tutorial/graduation_case_3.cpp
index 73893a106..432141d3e 100644
--- a/demos/tutorial/graduation_case_3.cpp
+++ b/demos/tutorial/graduation_case_3.cpp
@@ -139,7 +139,7 @@ int main(int argc, char* argv[])
                                    [&](std::size_t level, const auto& interval, const auto& index)
                                    {
                                        auto j    = index[0];
-                                       auto itag = static_cast<std::size_t>(interval.start + interval.index);
+                                       auto itag = interval.start + interval.index;
                                        for (int i = interval.start; i < interval.end; ++i)
                                        {
                                            if (tag[itag] && level < max_level)
diff --git a/include/samurai/algorithm.hpp b/include/samurai/algorithm.hpp
index 9147eda3d..5f73a8ee9 100644
--- a/include/samurai/algorithm.hpp
+++ b/include/samurai/algorithm.hpp
@@ -159,8 +159,9 @@ namespace samurai
     template <std::size_t dim, class TInterval, class Func>
     inline void for_each_cell(const LevelCellArray<dim, TInterval>& lca, Func&& f)
     {
-        using coord_index_t = typename TInterval::coord_index_t;
-        xt::xtensor_fixed<coord_index_t, xt::xshape<dim>> index;
+        using cell_t        = Cell<dim, TInterval>;
+        using index_value_t = typename cell_t::value_t;
+        typename cell_t::indices_t index;
 
         for (auto it = lca.cbegin(); it != lca.cend(); ++it)
         {
@@ -169,10 +170,10 @@ namespace samurai
                 index[d + 1] = it.index()[d];
             }
 
-            for (coord_index_t i = it->start; i < it->end; ++i)
+            for (index_value_t i = it->start; i < it->end; ++i)
             {
                 index[0] = i;
-                Cell<coord_index_t, dim> cell{lca.level(), index, static_cast<std::size_t>(it->index + i)};
+                cell_t cell{lca.level(), index, it->index + i};
                 f(cell);
             }
         }
@@ -181,10 +182,9 @@ namespace samurai
     template <std::size_t dim, class TInterval, class Func, class F, class... CT>
     inline void for_each_cell(const LevelCellArray<dim, TInterval>& lca, subset_operator<F, CT...> set, Func&& f)
     {
-        using set_t         = subset_operator<F, CT...>;
-        using coord_index_t = typename set_t::coord_index_t;
-
-        xt::xtensor_fixed<coord_index_t, xt::xshape<dim>> index;
+        using cell_t        = Cell<dim, TInterval>;
+        using index_value_t = typename cell_t::value_t;
+        typename cell_t::indices_t index;
 
         set(
             [&](const auto& interval, const auto& index_yz)
@@ -192,10 +192,10 @@ namespace samurai
                 index[0]                         = interval.start;
                 auto cell_index                  = lca.get_index(index);
                 xt::view(index, xt::range(1, _)) = index_yz;
-                for (coord_index_t i = interval.start; i < interval.end; ++i)
+                for (index_value_t i = interval.start; i < interval.end; ++i)
                 {
                     index[0] = i;
-                    Cell<coord_index_t, dim> cell{set.level(), index, cell_index++};
+                    cell_t cell{set.level(), index, cell_index++};
                     f(cell);
                 }
             });
@@ -223,18 +223,19 @@ namespace samurai
     template <class Mesh, class coord_type, class Func>
     inline void for_each_cell(const Mesh& mesh, std::size_t level, const typename Mesh::interval_t& i, const coord_type& index, Func&& f)
     {
-        using coord_index_t              = typename Mesh::interval_t::coord_index_t;
         static constexpr std::size_t dim = Mesh::dim;
+        using cell_t                     = Cell<dim, typename Mesh::interval_t>;
+        using index_value_t              = typename cell_t::value_t;
+        typename cell_t::indices_t coord;
 
-        xt::xtensor_fixed<coord_index_t, xt::xshape<dim>> coord;
         coord[0] = i.start;
         for (std::size_t d = 0; d < dim - 1; ++d)
         {
             coord[d + 1] = index[d];
         }
         auto cell_index = mesh.get_index(level, coord);
-        Cell<coord_index_t, dim> cell{level, coord, cell_index};
-        for (coord_index_t ii = 0; ii < static_cast<coord_index_t>(i.size()); ++ii)
+        cell_t cell{level, coord, cell_index};
+        for (index_value_t ii = 0; ii < static_cast<index_value_t>(i.size()); ++ii)
         {
             f(cell);
             cell.indices[0]++; // increment x coordinate
diff --git a/include/samurai/algorithm/update.hpp b/include/samurai/algorithm/update.hpp
index 574dfaa8a..c6d2d924a 100644
--- a/include/samurai/algorithm/update.hpp
+++ b/include/samurai/algorithm/update.hpp
@@ -489,7 +489,7 @@ namespace samurai
         for_each_interval(mesh[mesh_id_t::cells],
                           [&](std::size_t level, const auto& interval, const auto& index)
                           {
-                              auto itag = static_cast<std::size_t>(interval.start + interval.index);
+                              auto itag = interval.start + interval.index;
                               for (auto i = interval.start; i < interval.end; ++i)
                               {
                                   if (tag[itag] & static_cast<int>(CellFlag::refine))
@@ -555,7 +555,7 @@ namespace samurai
         for_each_interval(mesh[mesh_id_t::cells],
                           [&](std::size_t level, const auto& interval, const auto& index)
                           {
-                              auto itag = static_cast<std::size_t>(interval.start + interval.index);
+                              auto itag = interval.start + interval.index;
                               for (value_t i = interval.start; i < interval.end; ++i)
                               {
                                   if (tag[itag] & static_cast<int>(CellFlag::refine))
diff --git a/include/samurai/bc.hpp b/include/samurai/bc.hpp
index bcc61b6a6..ba3bc9907 100644
--- a/include/samurai/bc.hpp
+++ b/include/samurai/bc.hpp
@@ -80,11 +80,13 @@ namespace samurai
     ////////////////////////
     // BcValue definition //
     ////////////////////////
-    template <std::size_t dim, class T, std::size_t size>
+    template <class Field>
     struct BcValue
     {
-        using value_t  = detail::return_type_t<T, size>;
-        using coords_t = xt::xtensor_fixed<T, xt::xshape<dim>>;
+        static constexpr std::size_t dim = Field::dim;
+        using value_t                    = detail::return_type_t<typename Field::value_type, Field::size>;
+        using coords_t                   = xt::xtensor_fixed<double, xt::xshape<dim>>;
+        using cell_t                     = typename Field::cell_t;
 
         virtual ~BcValue()                 = default;
         BcValue(const BcValue&)            = delete;
@@ -92,30 +94,31 @@ namespace samurai
         BcValue(BcValue&&)                 = delete;
         BcValue& operator=(BcValue&&)      = delete;
 
-        virtual value_t get_value(const coords_t&) const = 0;
-        virtual std::unique_ptr<BcValue> clone() const   = 0;
-        virtual BCVType type() const                     = 0;
+        virtual value_t get_value(const cell_t&, const coords_t&) const = 0;
+        virtual std::unique_ptr<BcValue> clone() const                  = 0;
+        virtual BCVType type() const                                    = 0;
 
       protected:
 
         BcValue() = default;
     };
 
-    template <std::size_t dim, class T, std::size_t size>
-    class ConstantBc : public BcValue<dim, T, size>
+    template <class Field>
+    class ConstantBc : public BcValue<Field>
     {
       public:
 
-        using base_t   = BcValue<dim, T, size>;
+        using base_t   = BcValue<Field>;
         using value_t  = typename base_t::value_t;
         using coords_t = typename base_t::coords_t;
+        using cell_t   = typename base_t::cell_t;
 
         template <class... CT>
         ConstantBc(const CT... v);
 
         ConstantBc();
 
-        value_t get_value(const coords_t&) const override;
+        value_t get_value(const cell_t&, const coords_t&) const override;
         std::unique_ptr<base_t> clone() const override;
         BCVType type() const override;
 
@@ -124,19 +127,20 @@ namespace samurai
         value_t m_v;
     };
 
-    template <std::size_t dim, class T, std::size_t size>
-    class FunctionBc : public BcValue<dim, T, size>
+    template <class Field>
+    class FunctionBc : public BcValue<Field>
     {
       public:
 
-        using base_t     = BcValue<dim, T, size>;
+        using base_t     = BcValue<Field>;
         using value_t    = typename base_t::value_t;
         using coords_t   = typename base_t::coords_t;
-        using function_t = std::function<value_t(const coords_t&)>;
+        using cell_t     = typename base_t::cell_t;
+        using function_t = std::function<value_t(const cell_t&, const coords_t&)>;
 
         FunctionBc(const function_t& f);
 
-        value_t get_value(const coords_t& coords) const override;
+        value_t get_value(const cell_t& cell_in, const coords_t& coords) const override;
         std::unique_ptr<base_t> clone() const override;
         BCVType type() const override;
 
@@ -149,57 +153,57 @@ namespace samurai
     // BcValue implementation //
     ////////////////////////////
 
-    template <std::size_t dim, class T, std::size_t size>
+    template <class Field>
     template <class... CT>
-    ConstantBc<dim, T, size>::ConstantBc(const CT... v)
+    ConstantBc<Field>::ConstantBc(const CT... v)
         : m_v{v...}
     {
     }
 
-    template <std::size_t dim, class T, std::size_t size>
-    ConstantBc<dim, T, size>::ConstantBc()
+    template <class Field>
+    ConstantBc<Field>::ConstantBc()
     {
-        detail::fill(m_v, T{0});
+        detail::fill(m_v, typename Field::value_type{0});
     }
 
-    template <std::size_t dim, class T, std::size_t size>
-    inline auto ConstantBc<dim, T, size>::get_value(const coords_t&) const -> value_t
+    template <class Field>
+    inline auto ConstantBc<Field>::get_value(const cell_t&, const coords_t&) const -> value_t
     {
         return m_v;
     }
 
-    template <std::size_t dim, class T, std::size_t size>
-    auto ConstantBc<dim, T, size>::clone() const -> std::unique_ptr<base_t>
+    template <class Field>
+    auto ConstantBc<Field>::clone() const -> std::unique_ptr<base_t>
     {
         return std::make_unique<ConstantBc>(m_v);
     }
 
-    template <std::size_t dim, class T, std::size_t size>
-    inline BCVType ConstantBc<dim, T, size>::type() const
+    template <class Field>
+    inline BCVType ConstantBc<Field>::type() const
     {
         return BCVType::constant;
     }
 
-    template <std::size_t dim, class T, std::size_t size>
-    FunctionBc<dim, T, size>::FunctionBc(const function_t& f)
+    template <class Field>
+    FunctionBc<Field>::FunctionBc(const function_t& f)
         : m_func(f)
     {
     }
 
-    template <std::size_t dim, class T, std::size_t size>
-    inline auto FunctionBc<dim, T, size>::get_value(const coords_t& coords) const -> value_t
+    template <class Field>
+    inline auto FunctionBc<Field>::get_value(const cell_t& cell_in, const coords_t& coords) const -> value_t
     {
-        return m_func(coords);
+        return m_func(cell_in, coords);
     }
 
-    template <std::size_t dim, class T, std::size_t size>
-    auto FunctionBc<dim, T, size>::clone() const -> std::unique_ptr<base_t>
+    template <class Field>
+    auto FunctionBc<Field>::clone() const -> std::unique_ptr<base_t>
     {
         return std::make_unique<FunctionBc>(m_func);
     }
 
-    template <std::size_t dim, class T, std::size_t size>
-    inline BCVType FunctionBc<dim, T, size>::type() const
+    template <class Field>
+    inline BCVType FunctionBc<Field>::type() const
     {
         return BCVType::function;
     }
@@ -513,10 +517,11 @@ namespace samurai
         static constexpr std::size_t size = Field::size;
         using interval_t                  = typename Field::interval_t;
 
-        using bcvalue_t    = BcValue<dim, typename Field::value_type, size>;
+        using bcvalue_t    = BcValue<Field>;
         using bcvalue_impl = std::unique_ptr<bcvalue_t>;
         using value_t      = typename bcvalue_t::value_t;
         using coords_t     = typename bcvalue_t::coords_t;
+        using cell_t       = typename bcvalue_t::cell_t;
 
         using bcregion_t = BcRegion<dim, interval_t>;
         using lca_t      = typename bcregion_t::lca_t;
@@ -550,7 +555,7 @@ namespace samurai
                            xt::xtensor_fixed<typename interval_t::value_t, xt::xshape<dim - 1>> index);
 
         value_t constant_value();
-        value_t value(const coords_t& coords) const;
+        value_t value(const cell_t& cell_in, const coords_t& coords) const;
         const auto& value() const;
         BCVType get_value_type() const;
 
@@ -635,12 +640,16 @@ namespace samurai
         if (p_bcvalue->type() == BCVType::function)
         {
             coords_t coords;
+            cell_t cell_in;
+
             const double dx = 1. / (1 << level);
 
-            coords[0] = dx * i.start + 0.5 * (1 + dir[0]) * dx;
+            cell_in.level = level;
+            coords[0]     = dx * i.start + 0.5 * (1 + dir[0]) * dx;
             for (std::size_t d = 1; d < dim; ++d)
             {
-                coords[d] = dx * index[d - 1] + 0.5 * (1 + dir[d]) * dx;
+                coords[d]          = dx * index[d - 1] + 0.5 * (1 + dir[d]) * dx;
+                cell_in.indices[d] = index[d - 1] + dir[d];
             }
 
             if constexpr (size == 1)
@@ -652,10 +661,14 @@ namespace samurai
                 m_value.resize({i.size(), size});
             }
 
+            cell_in.indices[0] = i.start + dir[0];
+            cell_in.index      = i.index;
             for (std::size_t ii = 0; ii < i.size(); ++ii)
             {
-                xt::view(m_value, ii) = p_bcvalue->get_value(coords);
+                xt::view(m_value, ii) = p_bcvalue->get_value(cell_in, coords);
                 coords[0] += dx;
+                ++cell_in.indices[0];
+                ++cell_in.index;
             }
         }
     }
@@ -663,7 +676,7 @@ namespace samurai
     template <class Field>
     inline auto Bc<Field>::constant_value() -> value_t
     {
-        return p_bcvalue->get_value({});
+        return p_bcvalue->get_value({}, {});
     }
 
     template <class Field>
@@ -673,9 +686,9 @@ namespace samurai
     }
 
     template <class Field>
-    inline auto Bc<Field>::value(const coords_t& coords) const -> value_t
+    inline auto Bc<Field>::value(const cell_t& cell_in, const coords_t& coords) const -> value_t
     {
-        return p_bcvalue->get_value(coords);
+        return p_bcvalue->get_value(cell_in, coords);
     }
 
     template <class Field>
@@ -710,34 +723,22 @@ namespace samurai
     }
 
     template <template <class> class bc_type, class Field>
-    auto make_bc(Field& field, typename FunctionBc<Field::dim, typename Field::value_type, Field::size>::function_t func)
+    auto make_bc(Field& field, typename FunctionBc<Field>::function_t func)
     {
-        using value_t              = typename Field::value_type;
-        constexpr std::size_t dim  = Field::dim;
-        constexpr std::size_t size = Field::size;
-
         auto& mesh = detail::get_mesh(field.mesh());
-        return field.attach_bc(bc_type<Field>(mesh, FunctionBc<dim, value_t, size>(func)));
+        return field.attach_bc(bc_type<Field>(mesh, FunctionBc<Field>(func)));
     }
 
     template <template <class> class bc_type, class Field>
     auto make_bc(Field& field)
     {
-        using value_t              = typename Field::value_type;
-        constexpr std::size_t dim  = Field::dim;
-        constexpr std::size_t size = Field::size;
-
         auto& mesh = detail::get_mesh(field.mesh());
-        return field.attach_bc(bc_type<Field>(mesh, ConstantBc<dim, value_t, size>()));
+        return field.attach_bc(bc_type<Field>(mesh, ConstantBc<Field>()));
     }
 
     template <template <class> class bc_type, class Field, class... T>
     auto make_bc(Field& field, typename Field::value_type v1, T... v)
     {
-        using value_t              = typename Field::value_type;
-        constexpr std::size_t dim  = Field::dim;
-        constexpr std::size_t size = Field::size;
-
         static_assert(std::is_same_v<typename Field::value_type, std::common_type_t<typename Field::value_type, T...>>,
                       "The constant value type must be the same as the field value_type");
         static_assert(Field::size == sizeof...(T) + 1,
@@ -745,7 +746,7 @@ namespace samurai
                       "number of element in the field");
 
         auto& mesh = detail::get_mesh(field.mesh());
-        return field.attach_bc(bc_type<Field>(mesh, ConstantBc<dim, value_t, size>(v1, v...)));
+        return field.attach_bc(bc_type<Field>(mesh, ConstantBc<Field>(v1, v...)));
     }
 
     //////////////
diff --git a/include/samurai/cell.hpp b/include/samurai/cell.hpp
index 933f93604..794cb26d8 100644
--- a/include/samurai/cell.hpp
+++ b/include/samurai/cell.hpp
@@ -23,30 +23,33 @@ namespace samurai
      *  A cell is defined by its level, its integer coordinates,
      *  and its index in the data array.
      *
-     *  @tparam TCoord_index The type of the coordinates.
      *  @tparam dim_ The dimension of the cell.
+     *  @tparam TInterval The type of the interval.
      */
-    template <class TCoord_index, std::size_t dim_>
+    template <std::size_t dim_, class TInterval>
     struct Cell
     {
         static constexpr std::size_t dim = dim_;
-        using value_t                    = TCoord_index;
+        using interval_t                 = TInterval;
+        using value_t                    = typename interval_t::value_t;
+        using index_t                    = typename interval_t::index_t;
+        using indices_t                  = xt::xtensor_fixed<value_t, xt::xshape<dim>>;
+        using coords_t                   = xt::xtensor_fixed<double, xt::xshape<dim>>;
 
         Cell() = default;
 
-        template <class T>
-        Cell(std::size_t level, const T& indices, std::size_t index);
+        Cell(std::size_t level, const indices_t& indices, index_t index);
 
-        xt::xtensor_fixed<double, xt::xshape<dim>> corner() const;
+        coords_t corner() const;
         double corner(std::size_t i) const;
 
         /// The center of the cell.
-        xt::xtensor_fixed<double, xt::xshape<dim>> center() const;
+        coords_t center() const;
         double center(std::size_t i) const;
 
         /// The center of the face in the requested Cartesian direction.
         template <class Vector>
-        xt::xtensor_fixed<double, xt::xshape<dim>> face_center(const Vector& direction) const;
+        coords_t face_center(const Vector& direction) const;
 
         void to_stream(std::ostream& os) const;
 
@@ -54,18 +57,17 @@ namespace samurai
         std::size_t level = 0;
 
         /// The integer coordinates of the cell.
-        xt::xtensor_fixed<value_t, xt::xshape<dim>> indices;
+        indices_t indices;
 
         /// The index where the cell is in the data array.
-        std::size_t index = 0;
+        index_t index = 0;
 
         /// The length of the cell.
         double length = 0;
     };
 
-    template <class TCoord_index, std::size_t dim_>
-    template <class T>
-    inline Cell<TCoord_index, dim_>::Cell(std::size_t level_, const T& indices_, std::size_t index_)
+    template <std::size_t dim_, class TInterval>
+    inline Cell<dim_, TInterval>::Cell(std::size_t level_, const indices_t& indices_, index_t index_)
         : level(level_)
         , indices(indices_)
         , index(index_)
@@ -76,14 +78,14 @@ namespace samurai
     /**
      * The minimum corner of the cell.
      */
-    template <class TCoord_index, std::size_t dim_>
-    inline auto Cell<TCoord_index, dim_>::corner() const -> xt::xtensor_fixed<double, xt::xshape<dim>>
+    template <std::size_t dim_, class TInterval>
+    inline auto Cell<dim_, TInterval>::corner() const -> coords_t
     {
         return length * indices;
     }
 
-    template <class TCoord_index, std::size_t dim_>
-    inline double Cell<TCoord_index, dim_>::corner(std::size_t i) const
+    template <std::size_t dim_, class TInterval>
+    inline double Cell<dim_, TInterval>::corner(std::size_t i) const
     {
         return length * indices[i];
     }
@@ -91,34 +93,34 @@ namespace samurai
     /**
      * The minimum corner of the cell.
      */
-    template <class TCoord_index, std::size_t dim_>
-    inline auto Cell<TCoord_index, dim_>::center() const -> xt::xtensor_fixed<double, xt::xshape<dim>>
+    template <std::size_t dim_, class TInterval>
+    inline auto Cell<dim_, TInterval>::center() const -> coords_t
     {
         return length * (indices + 0.5);
     }
 
-    template <class TCoord_index, std::size_t dim_>
-    inline double Cell<TCoord_index, dim_>::center(std::size_t i) const
+    template <std::size_t dim_, class TInterval>
+    inline double Cell<dim_, TInterval>::center(std::size_t i) const
     {
         return length * (indices[i] + 0.5);
     }
 
-    template <class TCoord_index, std::size_t dim_>
+    template <std::size_t dim_, class TInterval>
     template <class Vector>
-    inline auto Cell<TCoord_index, dim_>::face_center(const Vector& direction) const -> xt::xtensor_fixed<double, xt::xshape<dim>>
+    inline auto Cell<dim_, TInterval>::face_center(const Vector& direction) const -> coords_t
     {
         assert(abs(xt::sum(direction)(0)) == 1); // We only want a Cartesian unit vector
         return center() + (length / 2) * direction;
     }
 
-    template <class TCoord_index, std::size_t dim_>
-    inline void Cell<TCoord_index, dim_>::to_stream(std::ostream& os) const
+    template <std::size_t dim_, class TInterval>
+    inline void Cell<dim_, TInterval>::to_stream(std::ostream& os) const
     {
         os << "Cell -> level: " << level << " indices: " << indices << " center: " << center() << " index: " << index;
     }
 
-    template <class TCoord_index, std::size_t dim>
-    inline std::ostream& operator<<(std::ostream& out, const Cell<TCoord_index, dim>& cell)
+    template <std::size_t dim, class TInterval>
+    inline std::ostream& operator<<(std::ostream& out, const Cell<dim, TInterval>& cell)
     {
         cell.to_stream(out);
         return out;
diff --git a/include/samurai/cell_array.hpp b/include/samurai/cell_array.hpp
index f1b9cafda..4fafc4885 100644
--- a/include/samurai/cell_array.hpp
+++ b/include/samurai/cell_array.hpp
@@ -73,6 +73,7 @@ namespace samurai
         using self_type  = CellArray<dim_, TInterval, max_size_>;
         using interval_t = TInterval;
         using value_t    = typename interval_t::value_t;
+        using index_t    = typename interval_t::index_t;
         using lca_type   = LevelCellArray<dim, TInterval>;
         using cl_type    = CellList<dim, TInterval, max_size>;
 
@@ -96,8 +97,8 @@ namespace samurai
         const interval_t& get_interval(std::size_t level, const xt::xtensor_fixed<value_t, xt::xshape<dim>>& coord) const;
 
         template <typename... T>
-        std::size_t get_index(std::size_t level, value_t i, T... index) const;
-        std::size_t get_index(std::size_t level, const xt::xtensor_fixed<value_t, xt::xshape<dim>>& coord) const;
+        index_t get_index(std::size_t level, value_t i, T... index) const;
+        index_t get_index(std::size_t level, const xt::xtensor_fixed<value_t, xt::xshape<dim>>& coord) const;
 
         std::size_t nb_cells() const;
         std::size_t nb_cells(std::size_t level) const;
@@ -279,14 +280,15 @@ namespace samurai
 
     template <std::size_t dim_, class TInterval, std::size_t max_size_>
     template <typename... T>
-    inline std::size_t CellArray<dim_, TInterval, max_size_>::get_index(std::size_t level, value_t i, T... index) const
+    inline auto CellArray<dim_, TInterval, max_size_>::get_index(std::size_t level, value_t i, T... index) const -> index_t
     {
         return m_cells[level].get_index(i, index...);
     }
 
     template <std::size_t dim_, class TInterval, std::size_t max_size_>
-    inline std::size_t
+    inline auto
     CellArray<dim_, TInterval, max_size_>::get_index(std::size_t level, const xt::xtensor_fixed<value_t, xt::xshape<dim>>& coord) const
+        -> index_t
     {
         return m_cells[level].get_index(coord);
     }
@@ -357,7 +359,6 @@ namespace samurai
     template <std::size_t dim_, class TInterval, std::size_t max_size_>
     inline void CellArray<dim_, TInterval, max_size_>::update_index()
     {
-        using index_t        = typename interval_t::index_t;
         std::size_t acc_size = 0;
         for_each_interval(*this,
                           [&](auto, auto& interval, auto)
diff --git a/include/samurai/field.hpp b/include/samurai/field.hpp
index 0395c917c..ea62606af 100644
--- a/include/samurai/field.hpp
+++ b/include/samurai/field.hpp
@@ -74,27 +74,28 @@ namespace samurai
         {
             static constexpr std::size_t dim = mesh_t::dim;
             using interval_t                 = typename mesh_t::interval_t;
-            using cell_t                     = Cell<typename interval_t::coord_index_t, dim>;
+            using index_t                    = typename interval_t::index_t;
+            using cell_t                     = Cell<dim, interval_t>;
             using data_type                  = xt::xtensor<value_t, 1>;
 
-            inline const value_t& operator[](std::size_t i) const
+            inline const value_t& operator[](index_t i) const
             {
-                return this->derived_cast().m_data[i];
+                return this->derived_cast().m_data[static_cast<std::size_t>(i)];
             }
 
-            inline value_t& operator[](std::size_t i)
+            inline value_t& operator[](index_t i)
             {
-                return this->derived_cast().m_data[i];
+                return this->derived_cast().m_data[static_cast<std::size_t>(i)];
             }
 
             inline const value_t& operator[](const cell_t& cell) const
             {
-                return this->derived_cast().m_data[cell.index];
+                return this->derived_cast().m_data[static_cast<std::size_t>(cell.index)];
             }
 
             inline value_t& operator[](const cell_t& cell)
             {
-                return this->derived_cast().m_data[cell.index];
+                return this->derived_cast().m_data[static_cast<std::size_t>(cell.index)];
             }
 
             inline const value_t& operator()(std::size_t i) const
@@ -139,17 +140,18 @@ namespace samurai
         {
             static constexpr std::size_t dim = mesh_t::dim;
             using interval_t                 = typename mesh_t::interval_t;
-            using cell_t                     = Cell<typename interval_t::coord_index_t, dim>;
+            using index_t                    = typename interval_t::index_t;
+            using cell_t                     = Cell<dim, interval_t>;
             using data_type                  = xt::xtensor<value_t, 2>;
 
-            inline auto operator[](std::size_t i) const
+            inline auto operator[](index_t i) const
             {
-                return xt::view(this->derived_cast().m_data, i);
+                return xt::view(this->derived_cast().m_data, static_cast<std::size_t>(i));
             }
 
-            inline auto operator[](std::size_t i)
+            inline auto operator[](index_t i)
             {
-                return xt::view(this->derived_cast().m_data, i);
+                return xt::view(this->derived_cast().m_data, static_cast<std::size_t>(i));
             }
 
             inline auto operator[](const cell_t& cell) const
@@ -235,7 +237,7 @@ namespace samurai
         {
             static constexpr std::size_t dim = mesh_t::dim;
             using interval_t                 = typename mesh_t::interval_t;
-            using cell_t                     = Cell<typename interval_t::coord_index_t, dim>;
+            using cell_t                     = Cell<dim, interval_t>;
             using data_type                  = xt::xtensor<value_t, 2>;
 
             inline auto operator[](std::size_t i) const
@@ -355,7 +357,7 @@ namespace samurai
 
         using inner_types::dim;
         using interval_t = typename mesh_t::interval_t;
-        using cell_t     = Cell<typename interval_t::coord_index_t, dim>;
+        using cell_t     = typename inner_types::cell_t;
 
         using iterator               = Field_iterator<self_type, false>;
         using const_iterator         = Field_iterator<const self_type, true>;
diff --git a/include/samurai/field_expression.hpp b/include/samurai/field_expression.hpp
index 40b42d23c..a91acee0b 100644
--- a/include/samurai/field_expression.hpp
+++ b/include/samurai/field_expression.hpp
@@ -91,8 +91,8 @@ namespace samurai
             return expr;
         }
 
-        template <class coord_index_t, std::size_t dim>
-        inline auto operator()(const Cell<coord_index_t, dim>& cell) const
+        template <std::size_t dim>
+        inline auto operator()(const Cell<dim, interval_t>& cell) const
         {
             return evaluate(std::make_index_sequence<sizeof...(CT)>(), cell);
         }
diff --git a/include/samurai/graduation.hpp b/include/samurai/graduation.hpp
index f204bb889..22dd5e220 100644
--- a/include/samurai/graduation.hpp
+++ b/include/samurai/graduation.hpp
@@ -75,7 +75,7 @@ namespace samurai
             for_each_interval(mesh,
                               [&](std::size_t level, const auto& interval, const auto& index_yz)
                               {
-                                  std::size_t itag = static_cast<std::size_t>(interval.start + interval.index);
+                                  auto itag = interval.start + interval.index;
                                   for (auto i = interval.start; i < interval.end; ++i, ++itag)
                                   {
                                       if (tag[itag])
diff --git a/include/samurai/interface.hpp b/include/samurai/interface.hpp
index 44e6dbb80..2b789cda2 100644
--- a/include/samurai/interface.hpp
+++ b/include/samurai/interface.hpp
@@ -133,8 +133,7 @@ namespace samurai
         static constexpr std::size_t dim = Mesh::dim;
         using mesh_id_t                  = typename Mesh::mesh_id_t;
         using mesh_interval_t            = typename Mesh::mesh_interval_t;
-        using coord_index_t              = typename Mesh::config::interval_t::coord_index_t;
-        using cell_t                     = typename samurai::Cell<coord_index_t, dim>;
+        using cell_t                     = Cell<dim, typename Mesh::interval_t>;
 
         if (level >= mesh.max_level())
         {
@@ -204,8 +203,7 @@ namespace samurai
         static constexpr std::size_t dim = Mesh::dim;
         using mesh_id_t                  = typename Mesh::mesh_id_t;
         using mesh_interval_t            = typename Mesh::mesh_interval_t;
-        using coord_index_t              = typename Mesh::config::interval_t::coord_index_t;
-        using cell_t                     = typename samurai::Cell<coord_index_t, dim>;
+        using cell_t                     = Cell<dim, typename Mesh::interval_t>;
 
         if (level >= mesh.max_level())
         {
diff --git a/include/samurai/io/from_geometry.hpp b/include/samurai/io/from_geometry.hpp
index dfa9b339c..5cd503a87 100644
--- a/include/samurai/io/from_geometry.hpp
+++ b/include/samurai/io/from_geometry.hpp
@@ -44,8 +44,8 @@ namespace samurai
 
         using Mesh = CellArray<3>;
         Mesh mesh;
-        using cell_t     = Cell<typename Mesh::value_t, 3>;
         using interval_t = typename Mesh::interval_t;
+        using cell_t     = Cell<3, interval_t>;
         double scale;
         std::tie(scale, mesh) = cgal::init_mesh<3>(start_level, cgal_mesh);
 
@@ -128,7 +128,7 @@ namespace samurai
                                   }
                                   else
                                   {
-                                      std::size_t itag = static_cast<std::size_t>(interval.start + interval.index);
+                                      auto itag = interval.start + interval.index;
                                       for (typename interval_t::value_t i = interval.start; i < interval.end; ++i, ++itag)
                                       {
                                           if ((tag[itag] & static_cast<int>(CellFlag::refine)) && level < max_level)
diff --git a/include/samurai/level_cell_array.hpp b/include/samurai/level_cell_array.hpp
index 82e1fb389..5c1ebce6b 100644
--- a/include/samurai/level_cell_array.hpp
+++ b/include/samurai/level_cell_array.hpp
@@ -111,8 +111,8 @@ namespace samurai
         const interval_t& get_interval(const xt::xtensor_fixed<value_t, xt::xshape<dim>>& coord) const;
 
         template <typename... T>
-        std::size_t get_index(const value_t& i, T... index) const;
-        std::size_t get_index(const xt::xtensor_fixed<value_t, xt::xshape<dim>>& coord) const;
+        index_t get_index(const value_t& i, T... index) const;
+        index_t get_index(const xt::xtensor_fixed<value_t, xt::xshape<dim>>& coord) const;
 
         void update_index();
 
@@ -467,15 +467,15 @@ namespace samurai
 
     template <std::size_t Dim, class TInterval>
     template <typename... T>
-    inline std::size_t LevelCellArray<Dim, TInterval>::get_index(const value_t& i, T... index) const
+    inline auto LevelCellArray<Dim, TInterval>::get_index(const value_t& i, T... index) const -> index_t
     {
-        return static_cast<std::size_t>(get_interval({i, i + 1}, index...).index + i);
+        return get_interval({i, i + 1}, index...).index + i;
     }
 
     template <std::size_t Dim, class TInterval>
-    inline std::size_t LevelCellArray<Dim, TInterval>::get_index(const xt::xtensor_fixed<value_t, xt::xshape<dim>>& coord) const
+    inline auto LevelCellArray<Dim, TInterval>::get_index(const xt::xtensor_fixed<value_t, xt::xshape<dim>>& coord) const -> index_t
     {
-        return static_cast<std::size_t>(get_interval(coord).index + coord[0]);
+        return get_interval(coord).index + coord[0];
     }
 
     /**
diff --git a/include/samurai/level_cell_list.hpp b/include/samurai/level_cell_list.hpp
index a5456088c..5266dace3 100644
--- a/include/samurai/level_cell_list.hpp
+++ b/include/samurai/level_cell_list.hpp
@@ -86,7 +86,7 @@ namespace samurai
 
         void to_stream(std::ostream& os) const;
 
-        void add_cell(const Cell<coord_index_t, dim>& cell);
+        void add_cell(const Cell<dim, interval_t>& cell);
 
       private:
 
@@ -152,7 +152,7 @@ namespace samurai
     }
 
     template <std::size_t Dim, class TInterval>
-    inline void LevelCellList<Dim, TInterval>::add_cell(const Cell<coord_index_t, dim>& cell)
+    inline void LevelCellList<Dim, TInterval>::add_cell(const Cell<dim, interval_t>& cell)
     {
         using namespace xt::placeholders;
 
diff --git a/include/samurai/mesh.hpp b/include/samurai/mesh.hpp
index b0e514341..d3ce0e389 100644
--- a/include/samurai/mesh.hpp
+++ b/include/samurai/mesh.hpp
@@ -48,6 +48,7 @@ namespace samurai
         using mesh_id_t  = typename config::mesh_id_t;
         using interval_t = typename config::interval_t;
         using value_t    = typename interval_t::value_t;
+        using index_t    = typename interval_t::index_t;
 
         using cl_type  = CellList<dim, interval_t, max_refinement_level>;
         using lcl_type = typename cl_type::lcl_type;
@@ -80,8 +81,8 @@ namespace samurai
         const interval_t& get_interval(std::size_t level, const xt::xtensor_fixed<value_t, xt::xshape<dim>>& coord) const;
 
         template <typename... T>
-        std::size_t get_index(std::size_t level, value_t i, T... index) const;
-        std::size_t get_index(std::size_t level, const xt::xtensor_fixed<value_t, xt::xshape<dim>>& coord) const;
+        index_t get_index(std::size_t level, value_t i, T... index) const;
+        index_t get_index(std::size_t level, const xt::xtensor_fixed<value_t, xt::xshape<dim>>& coord) const;
 
         void to_stream(std::ostream& os) const;
 
@@ -279,13 +280,13 @@ namespace samurai
 
     template <class D, class Config>
     template <typename... T>
-    inline std::size_t Mesh_base<D, Config>::get_index(std::size_t level, value_t i, T... index) const
+    inline auto Mesh_base<D, Config>::get_index(std::size_t level, value_t i, T... index) const -> index_t
     {
         return m_cells[mesh_id_t::reference].get_index(level, i, index...);
     }
 
     template <class D, class Config>
-    inline std::size_t Mesh_base<D, Config>::get_index(std::size_t level, const xt::xtensor_fixed<value_t, xt::xshape<dim>>& coord) const
+    inline auto Mesh_base<D, Config>::get_index(std::size_t level, const xt::xtensor_fixed<value_t, xt::xshape<dim>>& coord) const -> index_t
     {
         return m_cells[mesh_id_t::reference].get_index(level, coord);
     }
diff --git a/include/samurai/numeric/gauss_legendre.hpp b/include/samurai/numeric/gauss_legendre.hpp
index 00ad32c93..427223b9d 100644
--- a/include/samurai/numeric/gauss_legendre.hpp
+++ b/include/samurai/numeric/gauss_legendre.hpp
@@ -27,8 +27,8 @@ namespace samurai
             init();
         }
 
-        template <std::size_t func_result_size, class TCoord_index, std::size_t dim, class Func>
-        auto quadrature(const Cell<TCoord_index, dim>& cell, Func&& f)
+        template <std::size_t func_result_size, std::size_t dim, class TInterval, class Func>
+        auto quadrature(const Cell<dim, TInterval>& cell, Func&& f)
         {
             const double half_h = cell.length / 2;
             if constexpr (func_result_size == 1)
@@ -48,8 +48,8 @@ namespace samurai
 
       private:
 
-        template <class TCoord_index, std::size_t dim, class FuncResultType, class Func>
-        void compute_quadrature_sum(const Cell<TCoord_index, dim>& cell, FuncResultType& sum, Func&& f)
+        template <std::size_t dim, class TInterval, class FuncResultType, class Func>
+        void compute_quadrature_sum(const Cell<dim, TInterval>& cell, FuncResultType& sum, Func&& f)
         {
             static_assert(dim >= 1 && dim <= 3,
                           "The Gauss-Legendre quadrature is not implemented "
diff --git a/include/samurai/petsc/fv/FV_scheme_assembly.hpp b/include/samurai/petsc/fv/FV_scheme_assembly.hpp
index 9d6fbcb43..c27048d6b 100644
--- a/include/samurai/petsc/fv/FV_scheme_assembly.hpp
+++ b/include/samurai/petsc/fv/FV_scheme_assembly.hpp
@@ -35,6 +35,7 @@ namespace samurai
             using interval_t                                       = typename mesh_t::interval_t;
             using field_value_type                                 = typename field_t::value_type; // double
             using coord_index_t                                    = typename interval_t::coord_index_t;
+            using index_t                                          = typename interval_t::index_t;
             static constexpr std::size_t dim                       = field_t::dim;
             static constexpr std::size_t field_size                = field_t::size;
             static constexpr std::size_t output_field_size         = cfg_t::output_field_size;
@@ -43,7 +44,7 @@ namespace samurai
             static constexpr std::size_t bdry_stencil_size         = bdry_cfg_t::stencil_size;
             static constexpr std::size_t nb_bdry_ghosts            = bdry_cfg_t::nb_ghosts;
             static constexpr DirichletEnforcement dirichlet_enfcmt = bdry_cfg_t::dirichlet_enfcmt;
-            using cell_t                                           = samurai::Cell<coord_index_t, dim>;
+            using cell_t                                           = Cell<dim, interval_t>;
 
             using dirichlet_t = Dirichlet<field_t>;
             using neumann_t   = Neumann<field_t>;
@@ -61,9 +62,10 @@ namespace samurai
             std::vector<bool> m_is_row_empty;
 
             // Ghost recursion
-            using cell_coeff_pair_t     = std::pair<std::size_t, double>;
+            using cell_coeff_pair_t     = std::pair<index_t, double>;
             using CellLinearCombination = std::vector<cell_coeff_pair_t>;
-            std::map<std::size_t, CellLinearCombination> m_ghost_recursion;
+            using recursion_t           = std::map<index_t, CellLinearCombination>;
+            recursion_t m_ghost_recursion;
 
           public:
 
@@ -97,7 +99,7 @@ namespace samurai
             {
                 static constexpr PetscInt number_of_children = (1 << dim);
 
-                std::map<std::size_t, CellLinearCombination> recursion;
+                recursion_t recursion;
 
                 for_each_projection_ghost_and_children_cells<std::size_t>(
                     mesh(),
@@ -236,7 +238,7 @@ namespace samurai
                 }
                 else
                 {
-                    return m_col_shift + static_cast<PetscInt>(cell.index * field_size + field_j);
+                    return m_col_shift + static_cast<PetscInt>(cell.index * static_cast<index_t>(field_size) + field_j);
                 }
             }
 
@@ -253,7 +255,7 @@ namespace samurai
                 }
                 else
                 {
-                    return m_row_shift + static_cast<PetscInt>(cell.index * output_field_size + field_i);
+                    return m_row_shift + static_cast<PetscInt>(cell.index * static_cast<index_t>(output_field_size) + field_i);
                 }
             }
 
@@ -576,12 +578,12 @@ namespace samurai
                         double bc_value;
                         if constexpr (field_size == 1)
                         {
-                            bc_value = bc->value(boundary_point);
+                            bc_value = bc->value({}, boundary_point);
                         }
                         else
                         {
-                            bc_value = bc->value(boundary_point)(field_i); // TODO: call get_value() only once instead of
-                                                                           // once per field_i
+                            bc_value = bc->value({}, boundary_point)(field_i); // TODO: call get_value() only once instead of
+                                                                               // once per field_i
                         }
 
                         if constexpr (dirichlet_enfcmt == DirichletEnforcement::Elimination)
@@ -789,8 +791,6 @@ namespace samurai
 
             void assemble_prediction_1D(Mat& A)
             {
-                using index_t = int;
-
                 // double scalar = 1;
                 //  if constexpr (is_FluxBasedScheme_v<Scheme>)
                 //  {
@@ -824,7 +824,7 @@ namespace samurai
                                     double value           = -interpx[ci];
                                     auto coarse_cell_index = this->col_index(
                                         static_cast<PetscInt>(
-                                            this->mesh().get_index(ghost.level - 1, ig + static_cast<index_t>(ci - prediction_order))),
+                                            this->mesh().get_index(ghost.level - 1, ig + static_cast<coord_index_t>(ci - prediction_order))),
                                         field_i);
                                     MatSetValue(A, ghost_index, coarse_cell_index, scaling * value, INSERT_VALUES);
                                 }
@@ -837,8 +837,6 @@ namespace samurai
             template <class Cell>
             auto prediction_linear_combination_1D(const Cell& ghost)
             {
-                using index_t = int;
-
                 std::vector<cell_coeff_pair_t> linear_comb;
 
                 auto ii      = ghost.indices(0);
@@ -854,7 +852,8 @@ namespace samurai
                 {
                     if (ci != prediction_order)
                     {
-                        auto coarse_cell_index = this->mesh().get_index(ghost.level - 1, ig + static_cast<index_t>(ci - prediction_order));
+                        auto coarse_cell_index = this->mesh().get_index(ghost.level - 1,
+                                                                        ig + static_cast<coord_index_t>(ci - prediction_order));
                         linear_comb.emplace_back(coarse_cell_index, interpx[ci]);
                     }
                 }
@@ -863,7 +862,6 @@ namespace samurai
 
             void assemble_prediction_2D(Mat& A)
             {
-                using index_t = int;
                 samurai::for_each_prediction_ghost(
                     mesh(),
                     [&](auto& ghost)
@@ -896,11 +894,11 @@ namespace samurai
                                     if (ci != prediction_order || cj != prediction_order)
                                     {
                                         double value           = -interpx[ci] * interpy[cj];
-                                        auto coarse_cell_index = this->col_index(
-                                            static_cast<PetscInt>(this->mesh().get_index(ghost.level - 1,
-                                                                                         ig + static_cast<index_t>(ci - prediction_order),
-                                                                                         jg + static_cast<index_t>(cj - prediction_order))),
-                                            field_i);
+                                        auto coarse_cell_index = this->col_index(static_cast<PetscInt>(this->mesh().get_index(
+                                                                                     ghost.level - 1,
+                                                                                     ig + static_cast<coord_index_t>(ci - prediction_order),
+                                                                                     jg + static_cast<coord_index_t>(cj - prediction_order))),
+                                                                                 field_i);
                                         MatSetValue(A, ghost_index, coarse_cell_index, scaling * value, INSERT_VALUES);
                                     }
                                 }
@@ -913,8 +911,6 @@ namespace samurai
             template <class Cell>
             auto prediction_linear_combination_2D(const Cell& ghost)
             {
-                using index_t = int;
-
                 std::vector<cell_coeff_pair_t> linear_comb;
 
                 auto ii      = ghost.indices(0);
@@ -937,8 +933,8 @@ namespace samurai
                         if (ci != prediction_order || cj != prediction_order)
                         {
                             auto coarse_cell_index = this->mesh().get_index(ghost.level - 1,
-                                                                            ig + static_cast<index_t>(ci - prediction_order),
-                                                                            jg + static_cast<index_t>(cj - prediction_order));
+                                                                            ig + static_cast<coord_index_t>(ci - prediction_order),
+                                                                            jg + static_cast<coord_index_t>(cj - prediction_order));
                             linear_comb.emplace_back(coarse_cell_index, interpx[ci] * interpy[cj]);
                         }
                     }
@@ -948,7 +944,6 @@ namespace samurai
 
             void assemble_prediction_3D(Mat& A)
             {
-                using index_t = int;
                 samurai::for_each_prediction_ghost(
                     mesh(),
                     [&](auto& ghost)
@@ -987,12 +982,13 @@ namespace samurai
                                         if (ci != prediction_order || cj != prediction_order || ck != prediction_order)
                                         {
                                             double value           = -interpx[ci] * interpy[cj] * interpz[ck];
-                                            auto coarse_cell_index = this->col_index(static_cast<PetscInt>(this->mesh().get_index(
-                                                                                         ghost.level - 1,
-                                                                                         ig + static_cast<index_t>(ci - prediction_order),
-                                                                                         jg + static_cast<index_t>(cj - prediction_order),
-                                                                                         kg + static_cast<index_t>(ck - prediction_order))),
-                                                                                     field_i);
+                                            auto coarse_cell_index = this->col_index(
+                                                static_cast<PetscInt>(
+                                                    this->mesh().get_index(ghost.level - 1,
+                                                                           ig + static_cast<coord_index_t>(ci - prediction_order),
+                                                                           jg + static_cast<coord_index_t>(cj - prediction_order),
+                                                                           kg + static_cast<coord_index_t>(ck - prediction_order))),
+                                                field_i);
                                             MatSetValue(A, ghost_index, coarse_cell_index, scaling * value, INSERT_VALUES);
                                         }
                                     }
@@ -1006,8 +1002,6 @@ namespace samurai
             template <class Cell>
             auto prediction_linear_combination_3D(const Cell& ghost)
             {
-                using index_t = int;
-
                 std::vector<cell_coeff_pair_t> linear_comb;
 
                 auto ii      = ghost.indices(0);
@@ -1036,9 +1030,9 @@ namespace samurai
                             if (ci != prediction_order || cj != prediction_order || ck != prediction_order)
                             {
                                 auto coarse_cell_index = this->mesh().get_index(ghost.level - 1,
-                                                                                ig + static_cast<index_t>(ci - prediction_order),
-                                                                                jg + static_cast<index_t>(cj - prediction_order),
-                                                                                kg + static_cast<index_t>(ck - prediction_order));
+                                                                                ig + static_cast<coord_index_t>(ci - prediction_order),
+                                                                                jg + static_cast<coord_index_t>(cj - prediction_order),
+                                                                                kg + static_cast<coord_index_t>(ck - prediction_order));
                                 linear_comb.emplace_back(coarse_cell_index, interpx[ci] * interpy[cj] * interpz[ck]);
                             }
                         }
diff --git a/include/samurai/stencil.hpp b/include/samurai/stencil.hpp
index 20189b354..a0efc7bb4 100644
--- a/include/samurai/stencil.hpp
+++ b/include/samurai/stencil.hpp
@@ -70,13 +70,13 @@ namespace samurai
         static constexpr std::size_t dim = Mesh::dim;
         using mesh_interval_t            = typename Mesh::mesh_interval_t;
         using coord_index_t              = typename Mesh::config::interval_t::coord_index_t;
-        using Cell                       = typename samurai::Cell<coord_index_t, dim>;
+        using cell_t                     = Cell<dim, typename Mesh::interval_t>;
 
       private:
 
         const Mesh& m_mesh;                         // NOLINT(cppcoreguidelines-avoid-const-or-ref-data-members)
         const Stencil<stencil_size, dim> m_stencil; // NOLINT(cppcoreguidelines-avoid-const-or-ref-data-members)
-        std::array<Cell, stencil_size> m_cells;
+        std::array<cell_t, stencil_size> m_cells;
         unsigned int m_origin_cell;
 
       public:
@@ -93,14 +93,14 @@ namespace samurai
         void init(const mesh_interval_t& mesh_interval)
         {
             double length = cell_length(mesh_interval.level);
-            for (Cell& cell : m_cells)
+            for (cell_t& cell : m_cells)
             {
                 cell.level  = mesh_interval.level;
                 cell.length = length;
             }
 
             // origin of the stencil
-            Cell& origin_cell      = m_cells[m_origin_cell];
+            cell_t& origin_cell    = m_cells[m_origin_cell];
             origin_cell.indices[0] = mesh_interval.i.start;
             for (unsigned int d = 0; d < dim - 1; ++d)
             {
@@ -118,7 +118,7 @@ namespace samurai
                 auto d = xt::view(m_stencil, id);
 
                 // Translate the coordinates according the direction d
-                Cell& cell = m_cells[id];
+                cell_t& cell = m_cells[id];
                 for (unsigned int k = 0; k < dim; ++k)
                 {
                     cell.indices[k] = origin_cell.indices[k] + d[k];
@@ -137,14 +137,7 @@ namespace samurai
                 if (same_row) // same row as the stencil origin
                 {
                     // translation on the row
-                    if (d[0] >= 0)
-                    {
-                        cell.index = origin_cell.index + static_cast<std::size_t>(d[0]);
-                    }
-                    else
-                    {
-                        cell.index = origin_cell.index - static_cast<std::size_t>(-d[0]);
-                    }
+                    cell.index = origin_cell.index + d[0];
                 }
                 else
                 {
@@ -155,14 +148,14 @@ namespace samurai
 
         void move_next()
         {
-            for (Cell& cell : m_cells)
+            for (cell_t& cell : m_cells)
             {
                 cell.index++;      // increment cell index
                 cell.indices[0]++; // increment x-coordinate
             }
         }
 
-        std::array<Cell, stencil_size>& cells()
+        std::array<cell_t, stencil_size>& cells()
         {
             return m_cells;
         }
diff --git a/test/test_cell.cpp b/test/test_cell.cpp
index 7cd3d82b8..49e3373b7 100644
--- a/test/test_cell.cpp
+++ b/test/test_cell.cpp
@@ -3,20 +3,21 @@
 #include <xtensor/xfixed.hpp>
 
 #include <samurai/cell.hpp>
+#include <samurai/interval.hpp>
 
 namespace samurai
 {
     TEST(cell, length)
     {
         auto indices = xt::xtensor_fixed<int, xt::xshape<2>>({1, 1});
-        Cell<int, 2> c{1, indices, 0};
+        Cell<2, Interval<int>> c{1, indices, 0};
         EXPECT_EQ(c.length, 0.5);
     }
 
     TEST(cell, center)
     {
         auto indices = xt::xtensor_fixed<int, xt::xshape<2>>({1, 1});
-        Cell<int, 2> c{1, indices, 0};
+        Cell<2, Interval<int>> c{1, indices, 0};
         xt::xarray<double> expected{.75, .75};
         EXPECT_EQ(c.center(), expected);
     }
@@ -24,7 +25,7 @@ namespace samurai
     TEST(cell, first_corner)
     {
         auto indices = xt::xtensor_fixed<int, xt::xshape<2>>({1, 1});
-        Cell<int, 2> c{1, indices, 0};
+        Cell<2, Interval<int>> c{1, indices, 0};
         xt::xarray<double> expected{.5, .5};
         EXPECT_EQ(c.corner(), expected);
     }