diff --git a/CMakeLists.txt b/CMakeLists.txt
index 2c506705d..3e785bab9 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -3,7 +3,7 @@
 # SPDX-License-Identifier: MIT
 
 cmake_minimum_required(VERSION 3.22)
-project(DDC VERSION 0.3.0 LANGUAGES CXX)
+project(DDC VERSION 0.4.0 LANGUAGES CXX)
 
 # List of options
 
@@ -192,7 +192,7 @@ endif()
 
 if("${DDC_BUILD_KERNELS_SPLINES}")
     # Ginkgo
-    find_package(Ginkgo 1.8.0 EXACT REQUIRED)
+    find_package(Ginkgo 1.8...<2 REQUIRED)
 
     # Lapacke
     find_package(LAPACKE REQUIRED)
diff --git a/README.md b/README.md
index 0919a2b2f..76425fa39 100644
--- a/README.md
+++ b/README.md
@@ -47,7 +47,7 @@ To use DDC components, one needs the following dependencies:
 * (optional, IO interface) DDC::pdi
   * PDI 1.6...<2
 * (optional, spline interpolation) DDC::splines
-  * Ginkgo 1.8.0
+  * Ginkgo 1.8...<2
   * Kokkos Kernels 4.5.1...<5
 
 ## Getting the code and basic configuration
diff --git a/cmake/DDCConfig.cmake.in b/cmake/DDCConfig.cmake.in
index 811e34007..96237305c 100644
--- a/cmake/DDCConfig.cmake.in
+++ b/cmake/DDCConfig.cmake.in
@@ -24,7 +24,7 @@ if(@DDC_BUILD_KERNELS_FFT@)
 endif()
 
 if(@DDC_BUILD_KERNELS_SPLINES@)
-   ddc_find_dependency(Ginkgo 1.8.0 EXACT)
+   ddc_find_dependency(Ginkgo 1.8...<2)
    # DDC installs a FindLAPACKE.cmake file.
    # We choose to rely on it by prepending to CMAKE_MODULE_PATH
    # only the time of calling ddc_find_dependency.
diff --git a/include/ddc/discrete_domain.hpp b/include/ddc/discrete_domain.hpp
index dafb17d8c..925e0edac 100644
--- a/include/ddc/discrete_domain.hpp
+++ b/include/ddc/discrete_domain.hpp
@@ -264,10 +264,9 @@ class DiscreteDomain<>
 
     KOKKOS_DEFAULTED_FUNCTION constexpr DiscreteDomain() = default;
 
-    // Construct a DiscreteDomain from a reordered copy of `domain`
-    template <class... ODDims>
-    KOKKOS_FUNCTION constexpr explicit DiscreteDomain(
-            [[maybe_unused]] DiscreteDomain<ODDims...> const& domain)
+    /// Construct a DiscreteDomain by copies and merge of domains
+    template <class... DDoms, class = std::enable_if_t<(is_discrete_domain_v<DDoms> && ...)>>
+    KOKKOS_FUNCTION constexpr explicit DiscreteDomain([[maybe_unused]] DDoms const&... domains)
     {
     }
 
diff --git a/install_test/CMakeLists.txt b/install_test/CMakeLists.txt
index 457683337..e1cef912c 100644
--- a/install_test/CMakeLists.txt
+++ b/install_test/CMakeLists.txt
@@ -5,7 +5,7 @@
 cmake_minimum_required(VERSION 3.22)
 project(test-installed-ddc LANGUAGES CXX)
 
-find_package(DDC 0.3 REQUIRED COMPONENTS fft pdi splines)
+find_package(DDC 0.4 REQUIRED COMPONENTS fft pdi splines)
 
 message("DDC options:")
 message("DDC_BUILD_DOUBLE_PRECISION=${DDC_BUILD_DOUBLE_PRECISION}")
diff --git a/tests/splines/batched_2d_spline_builder.cpp b/tests/splines/batched_2d_spline_builder.cpp
index 34006a5be..54c97225a 100644
--- a/tests/splines/batched_2d_spline_builder.cpp
+++ b/tests/splines/batched_2d_spline_builder.cpp
@@ -100,9 +100,6 @@ struct DDimGPS : GrevillePoints<BSplines<X>>::interpolation_discrete_dimension_t
 {
 };
 
-template <typename X, typename I1, typename I2>
-using DDim = std::conditional_t<std::is_same_v<X, I1> || std::is_same_v<X, I2>, DDimGPS<X>, X>;
-
 #if defined(BC_PERIODIC)
 template <typename DDim1, typename DDim2>
 using evaluator_type = Evaluator2D::
@@ -121,10 +118,6 @@ using DVect = ddc::DiscreteVector<DDimX...>;
 template <typename... X>
 using Coord = ddc::Coordinate<X...>;
 
-// Extract batch dimensions from DDim (remove dimension of interest). Usefull
-template <typename I1, typename I2, typename... X>
-using BatchDims = ddc::type_seq_remove_t<ddc::detail::TypeSeq<X...>, ddc::detail::TypeSeq<I1, I2>>;
-
 // Templated function giving first coordinate of the mesh in given dimension.
 template <typename X>
 KOKKOS_FUNCTION Coord<X> x0()
@@ -171,34 +164,33 @@ void InterestDimInitializer(std::size_t const ncells)
 
 // Checks that when evaluating the spline at interpolation points one
 // recovers values that were used to build the spline
-template <typename ExecSpace, typename MemorySpace, typename I1, typename I2, typename... X>
+template <
+        typename ExecSpace,
+        typename MemorySpace,
+        typename DDimI1,
+        typename DDimI2,
+        typename... DDims>
 void Batched2dSplineTest()
 {
+    using I1 = typename DDimI1::continuous_dimension_type;
+    using I2 = typename DDimI2::continuous_dimension_type;
+
     // Instantiate execution spaces and initialize spaces
     ExecSpace const exec_space;
     std::size_t const ncells = 10;
-    InterestDimInitializer<DDim<I1, I1, I2>>(ncells);
-    InterestDimInitializer<DDim<I2, I1, I2>>(ncells);
+    InterestDimInitializer<DDimI1>(ncells);
+    InterestDimInitializer<DDimI2>(ncells);
 
     // Create the values domain (mesh)
-    ddc::DiscreteDomain<DDim<I1, I1, I2>> const interpolation_domain1
-            = GrevillePoints<BSplines<I1>>::template get_domain<DDim<I1, I1, I2>>();
-    ddc::DiscreteDomain<DDim<I2, I1, I2>> const interpolation_domain2
-            = GrevillePoints<BSplines<I2>>::template get_domain<DDim<I2, I1, I2>>();
-    ddc::DiscreteDomain<DDim<I1, I1, I2>, DDim<I2, I1, I2>> const
-            interpolation_domain(interpolation_domain1, interpolation_domain2);
-    // If we remove auto using the constructor syntax, nvcc does not compile
-    auto const dom_vals_tmp = ddc::DiscreteDomain<DDim<X, void, void>...>(
-            ddc::DiscreteDomain<DDim<
-                    X,
-                    void,
-                    void>>(DElem<DDim<X, void, void>>(0), DVect<DDim<X, void, void>>(ncells))...);
-    ddc::DiscreteDomain<DDim<X, I1, I2>...> const dom_vals
-            = ddc::replace_dim_of<DDim<I1, void, void>, DDim<I1, I1, I2>>(
-                    ddc::replace_dim_of<
-                            DDim<I2, void, void>,
-                            DDim<I2, I1, I2>>(dom_vals_tmp, interpolation_domain),
-                    interpolation_domain);
+    ddc::DiscreteDomain<DDimI1> const interpolation_domain1
+            = GrevillePoints<BSplines<I1>>::template get_domain<DDimI1>();
+    ddc::DiscreteDomain<DDimI2> const interpolation_domain2
+            = GrevillePoints<BSplines<I2>>::template get_domain<DDimI2>();
+    // The following line creates a discrete domain over all dimensions (DDims...) except DDimI1 and DDimI2.
+    auto const dom_vals_tmp = ddc::remove_dims_of_t<ddc::DiscreteDomain<DDims...>, DDimI1, DDimI2>(
+            ddc::DiscreteDomain<DDims>(DElem<DDims>(0), DVect<DDims>(ncells))...);
+    ddc::DiscreteDomain<DDims...> const
+            dom_vals(dom_vals_tmp, interpolation_domain1, interpolation_domain2);
 
 #if defined(BC_HERMITE)
     // Create the derivs domain
@@ -210,11 +202,10 @@ void Batched2dSplineTest()
             derivs_domain(derivs_domain1, derivs_domain2);
 
     auto const dom_derivs_1d
-            = ddc::replace_dim_of<DDim<I1, I1, I2>, ddc::Deriv<I1>>(dom_vals, derivs_domain1);
-    auto const dom_derivs2
-            = ddc::replace_dim_of<DDim<I2, I1, I2>, ddc::Deriv<I2>>(dom_vals, derivs_domain2);
+            = ddc::replace_dim_of<DDimI1, ddc::Deriv<I1>>(dom_vals, derivs_domain1);
+    auto const dom_derivs2 = ddc::replace_dim_of<DDimI2, ddc::Deriv<I2>>(dom_vals, derivs_domain2);
     auto const dom_derivs
-            = ddc::replace_dim_of<DDim<I2, I1, I2>, ddc::Deriv<I2>>(dom_derivs_1d, derivs_domain2);
+            = ddc::replace_dim_of<DDimI2, ddc::Deriv<I2>>(dom_derivs_1d, derivs_domain2);
 #endif
 
     // Create a SplineBuilder over BSplines<I> and batched along other dimensions using some boundary conditions
@@ -223,24 +214,24 @@ void Batched2dSplineTest()
             MemorySpace,
             BSplines<I1>,
             BSplines<I2>,
-            DDim<I1, I1, I2>,
-            DDim<I2, I1, I2>,
+            DDimI1,
+            DDimI2,
             s_bcl,
             s_bcr,
             s_bcl,
             s_bcr,
             ddc::SplineSolver::GINKGO,
-            DDim<X, I1, I2>...> const spline_builder(dom_vals);
+            DDims...> const spline_builder(dom_vals);
 
     // Compute usefull domains (dom_interpolation, dom_batch, dom_bsplines and dom_spline)
-    ddc::DiscreteDomain<DDim<I1, I1, I2>, DDim<I2, I1, I2>> const dom_interpolation
+    ddc::DiscreteDomain<DDimI1, DDimI2> const dom_interpolation
             = spline_builder.interpolation_domain();
     auto const dom_spline = spline_builder.batched_spline_domain();
 
     // Allocate and fill a chunk containing values to be passed as input to spline_builder. Those are values of cosine along interest dimension duplicated along batch dimensions
     ddc::Chunk vals_1d_host_alloc(dom_interpolation, ddc::HostAllocator<double>());
     ddc::ChunkSpan const vals_1d_host = vals_1d_host_alloc.span_view();
-    evaluator_type<DDim<I1, I1, I2>, DDim<I2, I1, I2>> const evaluator(dom_interpolation);
+    evaluator_type<DDimI1, DDimI2> const evaluator(dom_interpolation);
     evaluator(vals_1d_host);
     auto vals_1d_alloc = ddc::create_mirror_view_and_copy(exec_space, vals_1d_host);
     ddc::ChunkSpan const vals_1d = vals_1d_alloc.span_view();
@@ -250,8 +241,8 @@ void Batched2dSplineTest()
     ddc::parallel_for_each(
             exec_space,
             vals.domain(),
-            KOKKOS_LAMBDA(DElem<DDim<X, I1, I2>...> const e) {
-                vals(e) = vals_1d(DElem<DDim<I1, I1, I2>, DDim<I2, I1, I2>>(e));
+            KOKKOS_LAMBDA(DElem<DDims...> const e) {
+                vals(e) = vals_1d(DElem<DDimI1, DDimI2>(e));
             });
 
 #if defined(BC_HERMITE)
@@ -261,16 +252,14 @@ void Batched2dSplineTest()
     ddc::ChunkSpan const derivs_1d_lhs = derivs_1d_lhs_alloc.span_view();
     if (s_bcl == ddc::BoundCond::HERMITE) {
         ddc::Chunk derivs_1d_lhs1_host_alloc(
-                ddc::DiscreteDomain<
-                        ddc::Deriv<I1>,
-                        DDim<I2, I1, I2>>(derivs_domain1, interpolation_domain2),
+                ddc::DiscreteDomain<ddc::Deriv<I1>, DDimI2>(derivs_domain1, interpolation_domain2),
                 ddc::HostAllocator<double>());
         ddc::ChunkSpan const derivs_1d_lhs1_host = derivs_1d_lhs1_host_alloc.span_view();
         ddc::for_each(
                 derivs_1d_lhs1_host.domain(),
-                KOKKOS_LAMBDA(ddc::DiscreteElement<ddc::Deriv<I1>, DDim<I2, I1, I2>> const e) {
+                KOKKOS_LAMBDA(ddc::DiscreteElement<ddc::Deriv<I1>, DDimI2> const e) {
                     auto deriv_idx = ddc::DiscreteElement<ddc::Deriv<I1>>(e).uid();
-                    auto x2 = ddc::coordinate(ddc::DiscreteElement<DDim<I2, I1, I2>>(e));
+                    auto x2 = ddc::coordinate(ddc::DiscreteElement<DDimI2>(e));
                     derivs_1d_lhs1_host(e)
                             = evaluator.deriv(x0<I1>(), x2, deriv_idx + shift - 1, 0);
                 });
@@ -283,7 +272,7 @@ void Batched2dSplineTest()
                 derivs_1d_lhs.domain(),
                 KOKKOS_LAMBDA(
                         typename decltype(derivs_1d_lhs.domain())::discrete_element_type const e) {
-                    derivs_1d_lhs(e) = derivs_1d_lhs1(DElem<ddc::Deriv<I1>, DDim<I2, I1, I2>>(e));
+                    derivs_1d_lhs(e) = derivs_1d_lhs1(DElem<ddc::Deriv<I1>, DDimI2>(e));
                 });
     }
 
@@ -291,16 +280,14 @@ void Batched2dSplineTest()
     ddc::ChunkSpan const derivs_1d_rhs = derivs_1d_rhs_alloc.span_view();
     if (s_bcl == ddc::BoundCond::HERMITE) {
         ddc::Chunk derivs_1d_rhs1_host_alloc(
-                ddc::DiscreteDomain<
-                        ddc::Deriv<I1>,
-                        DDim<I2, I1, I2>>(derivs_domain1, interpolation_domain2),
+                ddc::DiscreteDomain<ddc::Deriv<I1>, DDimI2>(derivs_domain1, interpolation_domain2),
                 ddc::HostAllocator<double>());
         ddc::ChunkSpan const derivs_1d_rhs1_host = derivs_1d_rhs1_host_alloc.span_view();
         ddc::for_each(
                 derivs_1d_rhs1_host.domain(),
-                KOKKOS_LAMBDA(ddc::DiscreteElement<ddc::Deriv<I1>, DDim<I2, I1, I2>> const e) {
+                KOKKOS_LAMBDA(ddc::DiscreteElement<ddc::Deriv<I1>, DDimI2> const e) {
                     auto deriv_idx = ddc::DiscreteElement<ddc::Deriv<I1>>(e).uid();
-                    auto x2 = ddc::coordinate(ddc::DiscreteElement<DDim<I2, I1, I2>>(e));
+                    auto x2 = ddc::coordinate(ddc::DiscreteElement<DDimI2>(e));
                     derivs_1d_rhs1_host(e)
                             = evaluator.deriv(xN<I1>(), x2, deriv_idx + shift - 1, 0);
                 });
@@ -313,7 +300,7 @@ void Batched2dSplineTest()
                 derivs_1d_rhs.domain(),
                 KOKKOS_LAMBDA(
                         typename decltype(derivs_1d_rhs.domain())::discrete_element_type const e) {
-                    derivs_1d_rhs(e) = derivs_1d_rhs1(DElem<ddc::Deriv<I1>, DDim<I2, I1, I2>>(e));
+                    derivs_1d_rhs(e) = derivs_1d_rhs1(DElem<ddc::Deriv<I1>, DDimI2>(e));
                 });
     }
 
@@ -321,15 +308,13 @@ void Batched2dSplineTest()
     ddc::ChunkSpan const derivs2_lhs = derivs2_lhs_alloc.span_view();
     if (s_bcl == ddc::BoundCond::HERMITE) {
         ddc::Chunk derivs2_lhs1_host_alloc(
-                ddc::DiscreteDomain<
-                        DDim<I1, I1, I2>,
-                        ddc::Deriv<I2>>(interpolation_domain1, derivs_domain2),
+                ddc::DiscreteDomain<DDimI1, ddc::Deriv<I2>>(interpolation_domain1, derivs_domain2),
                 ddc::HostAllocator<double>());
         ddc::ChunkSpan const derivs2_lhs1_host = derivs2_lhs1_host_alloc.span_view();
         ddc::for_each(
                 derivs2_lhs1_host.domain(),
-                KOKKOS_LAMBDA(ddc::DiscreteElement<DDim<I1, I1, I2>, ddc::Deriv<I2>> const e) {
-                    auto x1 = ddc::coordinate(ddc::DiscreteElement<DDim<I1, I1, I2>>(e));
+                KOKKOS_LAMBDA(ddc::DiscreteElement<DDimI1, ddc::Deriv<I2>> const e) {
+                    auto x1 = ddc::coordinate(ddc::DiscreteElement<DDimI1>(e));
                     auto deriv_idx = ddc::DiscreteElement<ddc::Deriv<I2>>(e).uid();
                     derivs2_lhs1_host(e) = evaluator.deriv(x1, x0<I2>(), 0, deriv_idx + shift - 1);
                 });
@@ -342,7 +327,7 @@ void Batched2dSplineTest()
                 derivs2_lhs.domain(),
                 KOKKOS_LAMBDA(
                         typename decltype(derivs2_lhs.domain())::discrete_element_type const e) {
-                    derivs2_lhs(e) = derivs2_lhs1(DElem<DDim<I1, I1, I2>, ddc::Deriv<I2>>(e));
+                    derivs2_lhs(e) = derivs2_lhs1(DElem<DDimI1, ddc::Deriv<I2>>(e));
                 });
     }
 
@@ -350,15 +335,13 @@ void Batched2dSplineTest()
     ddc::ChunkSpan const derivs2_rhs = derivs2_rhs_alloc.span_view();
     if (s_bcl == ddc::BoundCond::HERMITE) {
         ddc::Chunk derivs2_rhs1_host_alloc(
-                ddc::DiscreteDomain<
-                        DDim<I1, I1, I2>,
-                        ddc::Deriv<I2>>(interpolation_domain1, derivs_domain2),
+                ddc::DiscreteDomain<DDimI1, ddc::Deriv<I2>>(interpolation_domain1, derivs_domain2),
                 ddc::HostAllocator<double>());
         ddc::ChunkSpan const derivs2_rhs1_host = derivs2_rhs1_host_alloc.span_view();
         ddc::for_each(
                 derivs2_rhs1_host.domain(),
-                KOKKOS_LAMBDA(ddc::DiscreteElement<DDim<I1, I1, I2>, ddc::Deriv<I2>> const e) {
-                    auto x1 = ddc::coordinate(ddc::DiscreteElement<DDim<I1, I1, I2>>(e));
+                KOKKOS_LAMBDA(ddc::DiscreteElement<DDimI1, ddc::Deriv<I2>> const e) {
+                    auto x1 = ddc::coordinate(ddc::DiscreteElement<DDimI1>(e));
                     auto deriv_idx = ddc::DiscreteElement<ddc::Deriv<I2>>(e).uid();
                     derivs2_rhs1_host(e) = evaluator.deriv(x1, xN<I2>(), 0, deriv_idx + shift - 1);
                 });
@@ -371,7 +354,7 @@ void Batched2dSplineTest()
                 derivs2_rhs.domain(),
                 KOKKOS_LAMBDA(
                         typename decltype(derivs2_rhs.domain())::discrete_element_type const e) {
-                    derivs2_rhs(e) = derivs2_rhs1(DElem<DDim<I1, I1, I2>, ddc::Deriv<I2>>(e));
+                    derivs2_rhs(e) = derivs2_rhs1(DElem<DDimI1, ddc::Deriv<I2>>(e));
                 });
     }
 
@@ -484,13 +467,13 @@ void Batched2dSplineTest()
             MemorySpace,
             BSplines<I1>,
             BSplines<I2>,
-            DDim<I1, I1, I2>,
-            DDim<I2, I1, I2>,
+            DDimI1,
+            DDimI2,
             extrapolation_rule_1_type,
             extrapolation_rule_1_type,
             extrapolation_rule_2_type,
             extrapolation_rule_2_type,
-            DDim<X, I1, I2>...> const
+            DDims...> const
             spline_evaluator(
                     extrapolation_rule_1,
                     extrapolation_rule_1,
@@ -503,8 +486,8 @@ void Batched2dSplineTest()
     ddc::parallel_for_each(
             exec_space,
             coords_eval.domain(),
-            KOKKOS_LAMBDA(DElem<DDim<X, I1, I2>...> const e) {
-                coords_eval(e) = ddc::coordinate(DElem<DDim<I1, I1, I2>, DDim<I2, I1, I2>>(e));
+            KOKKOS_LAMBDA(DElem<DDims...> const e) {
+                coords_eval(e) = ddc::coordinate(DElem<DDimI1, DDimI2>(e));
             });
 
 
@@ -534,7 +517,7 @@ void Batched2dSplineTest()
             spline_eval.domain(),
             0.,
             ddc::reducer::max<double>(),
-            KOKKOS_LAMBDA(DElem<DDim<X, I1, I2>...> const e) {
+            KOKKOS_LAMBDA(DElem<DDims...> const e) {
                 return Kokkos::abs(spline_eval(e) - vals(e));
             });
     double const max_norm_error_diff1 = ddc::parallel_transform_reduce(
@@ -542,9 +525,9 @@ void Batched2dSplineTest()
             spline_eval_deriv1.domain(),
             0.,
             ddc::reducer::max<double>(),
-            KOKKOS_LAMBDA(DElem<DDim<X, I1, I2>...> const e) {
-                Coord<I1> const x = ddc::coordinate(DElem<DDim<I1, I1, I2>>(e));
-                Coord<I2> const y = ddc::coordinate(DElem<DDim<I2, I1, I2>>(e));
+            KOKKOS_LAMBDA(DElem<DDims...> const e) {
+                Coord<I1> const x = ddc::coordinate(DElem<DDimI1>(e));
+                Coord<I2> const y = ddc::coordinate(DElem<DDimI2>(e));
                 return Kokkos::abs(spline_eval_deriv1(e) - evaluator.deriv(x, y, 1, 0));
             });
     double const max_norm_error_diff2 = ddc::parallel_transform_reduce(
@@ -552,9 +535,9 @@ void Batched2dSplineTest()
             spline_eval_deriv2.domain(),
             0.,
             ddc::reducer::max<double>(),
-            KOKKOS_LAMBDA(DElem<DDim<X, I1, I2>...> const e) {
-                Coord<I1> const x = ddc::coordinate(DElem<DDim<I1, I1, I2>>(e));
-                Coord<I2> const y = ddc::coordinate(DElem<DDim<I2, I1, I2>>(e));
+            KOKKOS_LAMBDA(DElem<DDims...> const e) {
+                Coord<I1> const x = ddc::coordinate(DElem<DDimI1>(e));
+                Coord<I2> const y = ddc::coordinate(DElem<DDimI2>(e));
                 return Kokkos::abs(spline_eval_deriv2(e) - evaluator.deriv(x, y, 0, 1));
             });
     double const max_norm_error_diff12 = ddc::parallel_transform_reduce(
@@ -562,9 +545,9 @@ void Batched2dSplineTest()
             spline_eval_deriv1.domain(),
             0.,
             ddc::reducer::max<double>(),
-            KOKKOS_LAMBDA(DElem<DDim<X, I1, I2>...> const e) {
-                Coord<I1> const x = ddc::coordinate(DElem<DDim<I1, I1, I2>>(e));
-                Coord<I2> const y = ddc::coordinate(DElem<DDim<I2, I1, I2>>(e));
+            KOKKOS_LAMBDA(DElem<DDims...> const e) {
+                Coord<I1> const x = ddc::coordinate(DElem<DDimI1>(e));
+                Coord<I2> const y = ddc::coordinate(DElem<DDimI2>(e));
                 return Kokkos::abs(spline_eval_deriv12(e) - evaluator.deriv(x, y, 1, 1));
             });
 
@@ -574,8 +557,7 @@ void Batched2dSplineTest()
     double const max_norm_diff2 = evaluator.max_norm(0, 1);
     double const max_norm_diff12 = evaluator.max_norm(1, 1);
 
-    SplineErrorBounds<evaluator_type<DDim<I1, I1, I2>, DDim<I2, I1, I2>>> const error_bounds(
-            evaluator);
+    SplineErrorBounds<evaluator_type<DDimI1, DDimI2>> const error_bounds(evaluator);
     EXPECT_LE(
             max_norm_error,
             std::
@@ -632,10 +614,10 @@ TEST(SUFFIX(Batched2dSplineHost), 2DXY)
     Batched2dSplineTest<
             Kokkos::DefaultHostExecutionSpace,
             Kokkos::DefaultHostExecutionSpace::memory_space,
-            DimX,
-            DimY,
-            DimX,
-            DimY>();
+            DDimGPS<DimX>,
+            DDimGPS<DimY>,
+            DDimGPS<DimX>,
+            DDimGPS<DimY>>();
 }
 
 TEST(SUFFIX(Batched2dSplineDevice), 2DXY)
@@ -643,10 +625,10 @@ TEST(SUFFIX(Batched2dSplineDevice), 2DXY)
     Batched2dSplineTest<
             Kokkos::DefaultExecutionSpace,
             Kokkos::DefaultExecutionSpace::memory_space,
-            DimX,
-            DimY,
-            DimX,
-            DimY>();
+            DDimGPS<DimX>,
+            DDimGPS<DimY>,
+            DDimGPS<DimX>,
+            DDimGPS<DimY>>();
 }
 
 TEST(SUFFIX(Batched2dSplineHost), 3DXY)
@@ -654,10 +636,10 @@ TEST(SUFFIX(Batched2dSplineHost), 3DXY)
     Batched2dSplineTest<
             Kokkos::DefaultHostExecutionSpace,
             Kokkos::DefaultHostExecutionSpace::memory_space,
-            DimX,
-            DimY,
-            DimX,
-            DimY,
+            DDimGPS<DimX>,
+            DDimGPS<DimY>,
+            DDimGPS<DimX>,
+            DDimGPS<DimY>,
             DDimBatch>();
 }
 
@@ -666,11 +648,11 @@ TEST(SUFFIX(Batched2dSplineHost), 3DXZ)
     Batched2dSplineTest<
             Kokkos::DefaultHostExecutionSpace,
             Kokkos::DefaultHostExecutionSpace::memory_space,
-            DimX,
-            DimZ,
-            DimX,
+            DDimGPS<DimX>,
+            DDimGPS<DimZ>,
+            DDimGPS<DimX>,
             DDimBatch,
-            DimZ>();
+            DDimGPS<DimZ>>();
 }
 
 TEST(SUFFIX(Batched2dSplineHost), 3DYZ)
@@ -678,11 +660,11 @@ TEST(SUFFIX(Batched2dSplineHost), 3DYZ)
     Batched2dSplineTest<
             Kokkos::DefaultHostExecutionSpace,
             Kokkos::DefaultHostExecutionSpace::memory_space,
-            DimY,
-            DimZ,
+            DDimGPS<DimY>,
+            DDimGPS<DimZ>,
             DDimBatch,
-            DimY,
-            DimZ>();
+            DDimGPS<DimY>,
+            DDimGPS<DimZ>>();
 }
 
 TEST(SUFFIX(Batched2dSplineDevice), 3DXY)
@@ -690,10 +672,10 @@ TEST(SUFFIX(Batched2dSplineDevice), 3DXY)
     Batched2dSplineTest<
             Kokkos::DefaultExecutionSpace,
             Kokkos::DefaultExecutionSpace::memory_space,
-            DimX,
-            DimY,
-            DimX,
-            DimY,
+            DDimGPS<DimX>,
+            DDimGPS<DimY>,
+            DDimGPS<DimX>,
+            DDimGPS<DimY>,
             DDimBatch>();
 }
 
@@ -702,11 +684,11 @@ TEST(SUFFIX(Batched2dSplineDevice), 3DXZ)
     Batched2dSplineTest<
             Kokkos::DefaultExecutionSpace,
             Kokkos::DefaultExecutionSpace::memory_space,
-            DimX,
-            DimZ,
-            DimX,
+            DDimGPS<DimX>,
+            DDimGPS<DimZ>,
+            DDimGPS<DimX>,
             DDimBatch,
-            DimZ>();
+            DDimGPS<DimZ>>();
 }
 
 TEST(SUFFIX(Batched2dSplineDevice), 3DYZ)
@@ -714,9 +696,9 @@ TEST(SUFFIX(Batched2dSplineDevice), 3DYZ)
     Batched2dSplineTest<
             Kokkos::DefaultExecutionSpace,
             Kokkos::DefaultExecutionSpace::memory_space,
-            DimY,
-            DimZ,
+            DDimGPS<DimY>,
+            DDimGPS<DimZ>,
             DDimBatch,
-            DimY,
-            DimZ>();
+            DDimGPS<DimY>,
+            DDimGPS<DimZ>>();
 }
diff --git a/tests/splines/batched_spline_builder.cpp b/tests/splines/batched_spline_builder.cpp
index cad56aa05..f0d3b5a2a 100644
--- a/tests/splines/batched_spline_builder.cpp
+++ b/tests/splines/batched_spline_builder.cpp
@@ -113,9 +113,6 @@ struct DDimGPS : GrevillePoints<BSplines<X>>::interpolation_discrete_dimension_t
 {
 };
 
-template <typename X, typename Y>
-using DDim = std::conditional_t<std::is_same_v<X, Y>, DDimGPS<X>, X>;
-
 template <typename DDimX>
 using evaluator_type = CosineEvaluator::Evaluator<DDimX>;
 
@@ -126,10 +123,6 @@ using DVect = ddc::DiscreteVector<DDimX...>;
 template <typename... X>
 using Coord = ddc::Coordinate<X...>;
 
-// Extract batch dimensions from DDim (remove dimension of interest). Usefull
-template <typename X, typename... Y>
-using BatchDims = ddc::type_seq_remove_t<ddc::detail::TypeSeq<Y...>, ddc::detail::TypeSeq<X>>;
-
 // Templated function giving first coordinate of the mesh in given dimension.
 template <typename X>
 KOKKOS_FUNCTION Coord<X> x0()
@@ -176,30 +169,30 @@ void InterestDimInitializer(std::size_t const ncells)
 
 // Checks that when evaluating the spline at interpolation points one
 // recovers values that were used to build the spline
-template <typename ExecSpace, typename MemorySpace, typename I, typename... X>
+template <typename ExecSpace, typename MemorySpace, typename DDimI, typename... DDims>
 void BatchedSplineTest()
 {
+    using I = typename DDimI::continuous_dimension_type;
+
     // Instantiate execution spaces and initialize spaces
     ExecSpace const exec_space;
 
     std::size_t const ncells = 10;
-    InterestDimInitializer<DDim<I, I>>(ncells);
+    InterestDimInitializer<DDimI>(ncells);
 
     // Create the values domain (mesh)
-    ddc::DiscreteDomain<DDim<I, I>> const interpolation_domain
-            = GrevillePoints<BSplines<I>>::template get_domain<DDim<I, I>>();
-    // If we remove auto using the constructor syntax, nvcc does not compile
-    auto const dom_vals_tmp = ddc::DiscreteDomain<DDim<X, void>...>(
-            ddc::DiscreteDomain<
-                    DDim<X, void>>(DElem<DDim<X, void>>(0), DVect<DDim<X, void>>(ncells))...);
-    ddc::DiscreteDomain<DDim<X, I>...> const dom_vals
-            = ddc::replace_dim_of<DDim<I, void>, DDim<I, I>>(dom_vals_tmp, interpolation_domain);
+    ddc::DiscreteDomain<DDimI> const interpolation_domain
+            = GrevillePoints<BSplines<I>>::template get_domain<DDimI>();
+    // The following line creates a discrete domain over all dimensions (DDims...) except DDimI.
+    auto const dom_vals_tmp = ddc::remove_dims_of_t<ddc::DiscreteDomain<DDims...>, DDimI>(
+            ddc::DiscreteDomain<DDims>(DElem<DDims>(0), DVect<DDims>(ncells))...);
+    ddc::DiscreteDomain<DDims...> const dom_vals(dom_vals_tmp, interpolation_domain);
 
 #if defined(BC_HERMITE)
     // Create the derivs domain
     ddc::DiscreteDomain<ddc::Deriv<I>> const
             derivs_domain(DElem<ddc::Deriv<I>>(1), DVect<ddc::Deriv<I>>(s_degree_x / 2));
-    auto const dom_derivs = ddc::replace_dim_of<DDim<I, I>, ddc::Deriv<I>>(dom_vals, derivs_domain);
+    auto const dom_derivs = ddc::replace_dim_of<DDimI, ddc::Deriv<I>>(dom_vals, derivs_domain);
 #endif
 
     // Create a SplineBuilder over BSplines<I> and batched along other dimensions using some boundary conditions
@@ -207,7 +200,7 @@ void BatchedSplineTest()
             ExecSpace,
             MemorySpace,
             BSplines<I>,
-            DDim<I, I>,
+            DDimI,
             s_bcl,
             s_bcr,
 #if defined(SOLVER_LAPACK)
@@ -215,17 +208,17 @@ void BatchedSplineTest()
 #elif defined(SOLVER_GINKGO)
             ddc::SplineSolver::GINKGO,
 #endif
-            DDim<X, I>...> const spline_builder(dom_vals);
+            DDims...> const spline_builder(dom_vals);
 
     // Compute usefull domains (dom_interpolation, dom_batch, dom_bsplines and dom_spline)
-    ddc::DiscreteDomain<DDim<I, I>> const dom_interpolation = spline_builder.interpolation_domain();
+    ddc::DiscreteDomain<DDimI> const dom_interpolation = spline_builder.interpolation_domain();
     auto const dom_batch = spline_builder.batch_domain();
     auto const dom_spline = spline_builder.batched_spline_domain();
 
     // Allocate and fill a chunk containing values to be passed as input to spline_builder. Those are values of cosine along interest dimension duplicated along batch dimensions
     ddc::Chunk vals_1d_host_alloc(dom_interpolation, ddc::HostAllocator<double>());
     ddc::ChunkSpan const vals_1d_host = vals_1d_host_alloc.span_view();
-    evaluator_type<DDim<I, I>> const evaluator(dom_interpolation);
+    evaluator_type<DDimI> const evaluator(dom_interpolation);
     evaluator(vals_1d_host);
     auto vals_1d_alloc = ddc::create_mirror_view_and_copy(exec_space, vals_1d_host);
     ddc::ChunkSpan const vals_1d = vals_1d_alloc.span_view();
@@ -235,9 +228,7 @@ void BatchedSplineTest()
     ddc::parallel_for_each(
             exec_space,
             vals.domain(),
-            KOKKOS_LAMBDA(DElem<DDim<X, I>...> const e) {
-                vals(e) = vals_1d(DElem<DDim<I, I>>(e));
-            });
+            KOKKOS_LAMBDA(DElem<DDims...> const e) { vals(e) = vals_1d(DElem<DDimI>(e)); });
 
 #if defined(BC_HERMITE)
     // Allocate and fill a chunk containing derivs to be passed as input to spline_builder.
@@ -315,10 +306,10 @@ void BatchedSplineTest()
             ExecSpace,
             MemorySpace,
             BSplines<I>,
-            DDim<I, I>,
+            DDimI,
             extrapolation_rule_type,
             extrapolation_rule_type,
-            DDim<X, I>...> const spline_evaluator_batched(extrapolation_rule, extrapolation_rule);
+            DDims...> const spline_evaluator_batched(extrapolation_rule, extrapolation_rule);
 
     // Instantiate chunk of coordinates of dom_interpolation
     ddc::Chunk coords_eval_alloc(dom_vals, ddc::KokkosAllocator<Coord<I>, MemorySpace>());
@@ -326,8 +317,8 @@ void BatchedSplineTest()
     ddc::parallel_for_each(
             exec_space,
             coords_eval.domain(),
-            KOKKOS_LAMBDA(DElem<DDim<X, I>...> const e) {
-                coords_eval(e) = ddc::coordinate(DElem<DDim<I, I>>(e));
+            KOKKOS_LAMBDA(DElem<DDims...> const e) {
+                coords_eval(e) = ddc::coordinate(DElem<DDimI>(e));
             });
 
 
@@ -350,7 +341,7 @@ void BatchedSplineTest()
             spline_eval.domain(),
             0.,
             ddc::reducer::max<double>(),
-            KOKKOS_LAMBDA(DElem<DDim<X, I>...> const e) {
+            KOKKOS_LAMBDA(DElem<DDims...> const e) {
                 return Kokkos::abs(spline_eval(e) - vals(e));
             });
 
@@ -359,8 +350,8 @@ void BatchedSplineTest()
             spline_eval_deriv.domain(),
             0.,
             ddc::reducer::max<double>(),
-            KOKKOS_LAMBDA(DElem<DDim<X, I>...> const e) {
-                Coord<I> const x = ddc::coordinate(DElem<DDim<I, I>>(e));
+            KOKKOS_LAMBDA(DElem<DDims...> const e) {
+                Coord<I> const x = ddc::coordinate(DElem<DDimI>(e));
                 return Kokkos::abs(spline_eval_deriv(e) - evaluator.deriv(x, 1));
             });
     double const max_norm_error_integ = ddc::parallel_transform_reduce(
@@ -379,7 +370,7 @@ void BatchedSplineTest()
     double const max_norm_diff = evaluator.max_norm(1);
     double const max_norm_int = evaluator.max_norm(-1);
 
-    SplineErrorBounds<evaluator_type<DDim<I, I>>> const error_bounds(evaluator);
+    SplineErrorBounds<evaluator_type<DDimI>> const error_bounds(evaluator);
     EXPECT_LE(
             max_norm_error,
             std::max(error_bounds.error_bound(dx<I>(ncells), s_degree_x), 1.0e-14 * max_norm));
@@ -428,8 +419,8 @@ TEST(SUFFIX(BatchedSplineHost), 1DX)
     BatchedSplineTest<
             Kokkos::DefaultHostExecutionSpace,
             Kokkos::DefaultHostExecutionSpace::memory_space,
-            DimX,
-            DimX>();
+            DDimGPS<DimX>,
+            DDimGPS<DimX>>();
 }
 
 TEST(SUFFIX(BatchedSplineDevice), 1DX)
@@ -437,8 +428,8 @@ TEST(SUFFIX(BatchedSplineDevice), 1DX)
     BatchedSplineTest<
             Kokkos::DefaultExecutionSpace,
             Kokkos::DefaultExecutionSpace::memory_space,
-            DimX,
-            DimX>();
+            DDimGPS<DimX>,
+            DDimGPS<DimX>>();
 }
 
 TEST(SUFFIX(BatchedSplineHost), 2DX)
@@ -446,8 +437,8 @@ TEST(SUFFIX(BatchedSplineHost), 2DX)
     BatchedSplineTest<
             Kokkos::DefaultHostExecutionSpace,
             Kokkos::DefaultHostExecutionSpace::memory_space,
-            DimX,
-            DimX,
+            DDimGPS<DimX>,
+            DDimGPS<DimX>,
             DDimBatch1>();
 }
 
@@ -456,9 +447,9 @@ TEST(SUFFIX(BatchedSplineHost), 2DY)
     BatchedSplineTest<
             Kokkos::DefaultHostExecutionSpace,
             Kokkos::DefaultHostExecutionSpace::memory_space,
-            DimY,
+            DDimGPS<DimY>,
             DDimBatch1,
-            DimY>();
+            DDimGPS<DimY>>();
 }
 
 TEST(SUFFIX(BatchedSplineDevice), 2DX)
@@ -466,8 +457,8 @@ TEST(SUFFIX(BatchedSplineDevice), 2DX)
     BatchedSplineTest<
             Kokkos::DefaultExecutionSpace,
             Kokkos::DefaultExecutionSpace::memory_space,
-            DimX,
-            DimX,
+            DDimGPS<DimX>,
+            DDimGPS<DimX>,
             DDimBatch1>();
 }
 
@@ -476,9 +467,9 @@ TEST(SUFFIX(BatchedSplineDevice), 2DY)
     BatchedSplineTest<
             Kokkos::DefaultExecutionSpace,
             Kokkos::DefaultExecutionSpace::memory_space,
-            DimY,
+            DDimGPS<DimY>,
             DDimBatch1,
-            DimY>();
+            DDimGPS<DimY>>();
 }
 
 TEST(SUFFIX(BatchedSplineHost), 3DX)
@@ -486,8 +477,8 @@ TEST(SUFFIX(BatchedSplineHost), 3DX)
     BatchedSplineTest<
             Kokkos::DefaultHostExecutionSpace,
             Kokkos::DefaultHostExecutionSpace::memory_space,
-            DimX,
-            DimX,
+            DDimGPS<DimX>,
+            DDimGPS<DimX>,
             DDimBatch1,
             DDimBatch2>();
 }
@@ -497,9 +488,9 @@ TEST(SUFFIX(BatchedSplineHost), 3DY)
     BatchedSplineTest<
             Kokkos::DefaultHostExecutionSpace,
             Kokkos::DefaultHostExecutionSpace::memory_space,
-            DimY,
+            DDimGPS<DimY>,
             DDimBatch1,
-            DimY,
+            DDimGPS<DimY>,
             DDimBatch2>();
 }
 
@@ -508,10 +499,10 @@ TEST(SUFFIX(BatchedSplineHost), 3DZ)
     BatchedSplineTest<
             Kokkos::DefaultHostExecutionSpace,
             Kokkos::DefaultHostExecutionSpace::memory_space,
-            DimZ,
+            DDimGPS<DimZ>,
             DDimBatch1,
             DDimBatch2,
-            DimZ>();
+            DDimGPS<DimZ>>();
 }
 
 TEST(SUFFIX(BatchedSplineDevice), 3DX)
@@ -519,8 +510,8 @@ TEST(SUFFIX(BatchedSplineDevice), 3DX)
     BatchedSplineTest<
             Kokkos::DefaultExecutionSpace,
             Kokkos::DefaultExecutionSpace::memory_space,
-            DimX,
-            DimX,
+            DDimGPS<DimX>,
+            DDimGPS<DimX>,
             DDimBatch1,
             DDimBatch2>();
 }
@@ -530,9 +521,9 @@ TEST(SUFFIX(BatchedSplineDevice), 3DY)
     BatchedSplineTest<
             Kokkos::DefaultExecutionSpace,
             Kokkos::DefaultExecutionSpace::memory_space,
-            DimY,
+            DDimGPS<DimY>,
             DDimBatch1,
-            DimY,
+            DDimGPS<DimY>,
             DDimBatch2>();
 }
 
@@ -541,10 +532,10 @@ TEST(SUFFIX(BatchedSplineDevice), 3DZ)
     BatchedSplineTest<
             Kokkos::DefaultExecutionSpace,
             Kokkos::DefaultExecutionSpace::memory_space,
-            DimZ,
+            DDimGPS<DimZ>,
             DDimBatch1,
             DDimBatch2,
-            DimZ>();
+            DDimGPS<DimZ>>();
 }
 
 
@@ -553,8 +544,8 @@ TEST(SUFFIX(BatchedSplineHost), 4DX)
     BatchedSplineTest<
             Kokkos::DefaultHostExecutionSpace,
             Kokkos::DefaultHostExecutionSpace::memory_space,
-            DimX,
-            DimX,
+            DDimGPS<DimX>,
+            DDimGPS<DimX>,
             DDimBatch1,
             DDimBatch2,
             DDimBatch3>();
@@ -565,9 +556,9 @@ TEST(SUFFIX(BatchedSplineHost), 4DY)
     BatchedSplineTest<
             Kokkos::DefaultHostExecutionSpace,
             Kokkos::DefaultHostExecutionSpace::memory_space,
-            DimY,
+            DDimGPS<DimY>,
             DDimBatch1,
-            DimY,
+            DDimGPS<DimY>,
             DDimBatch2,
             DDimBatch3>();
 }
@@ -577,10 +568,10 @@ TEST(SUFFIX(BatchedSplineHost), 4DZ)
     BatchedSplineTest<
             Kokkos::DefaultHostExecutionSpace,
             Kokkos::DefaultHostExecutionSpace::memory_space,
-            DimZ,
+            DDimGPS<DimZ>,
             DDimBatch1,
             DDimBatch2,
-            DimZ,
+            DDimGPS<DimZ>,
             DDimBatch3>();
 }
 
@@ -589,11 +580,11 @@ TEST(SUFFIX(BatchedSplineHost), 4DT)
     BatchedSplineTest<
             Kokkos::DefaultHostExecutionSpace,
             Kokkos::DefaultHostExecutionSpace::memory_space,
-            DimT,
+            DDimGPS<DimT>,
             DDimBatch1,
             DDimBatch2,
             DDimBatch3,
-            DimT>();
+            DDimGPS<DimT>>();
 }
 
 TEST(SUFFIX(BatchedSplineDevice), 4DX)
@@ -601,8 +592,8 @@ TEST(SUFFIX(BatchedSplineDevice), 4DX)
     BatchedSplineTest<
             Kokkos::DefaultExecutionSpace,
             Kokkos::DefaultExecutionSpace::memory_space,
-            DimX,
-            DimX,
+            DDimGPS<DimX>,
+            DDimGPS<DimX>,
             DDimBatch1,
             DDimBatch2,
             DDimBatch3>();
@@ -613,9 +604,9 @@ TEST(SUFFIX(BatchedSplineDevice), 4DY)
     BatchedSplineTest<
             Kokkos::DefaultExecutionSpace,
             Kokkos::DefaultExecutionSpace::memory_space,
-            DimY,
+            DDimGPS<DimY>,
             DDimBatch1,
-            DimY,
+            DDimGPS<DimY>,
             DDimBatch2,
             DDimBatch3>();
 }
@@ -625,10 +616,10 @@ TEST(SUFFIX(BatchedSplineDevice), 4DZ)
     BatchedSplineTest<
             Kokkos::DefaultExecutionSpace,
             Kokkos::DefaultExecutionSpace::memory_space,
-            DimZ,
+            DDimGPS<DimZ>,
             DDimBatch1,
             DDimBatch2,
-            DimZ,
+            DDimGPS<DimZ>,
             DDimBatch3>();
 }
 
@@ -637,9 +628,9 @@ TEST(SUFFIX(BatchedSplineDevice), 4DT)
     BatchedSplineTest<
             Kokkos::DefaultExecutionSpace,
             Kokkos::DefaultExecutionSpace::memory_space,
-            DimT,
+            DDimGPS<DimT>,
             DDimBatch1,
             DDimBatch2,
             DDimBatch3,
-            DimT>();
+            DDimGPS<DimT>>();
 }
diff --git a/tests/splines/extrapolation_rule.cpp b/tests/splines/extrapolation_rule.cpp
index 3620a118a..cb64624fb 100644
--- a/tests/splines/extrapolation_rule.cpp
+++ b/tests/splines/extrapolation_rule.cpp
@@ -86,12 +86,6 @@ struct DDimGPS2 : GrevillePoints2<BSplines<X>>::interpolation_discrete_dimension
 {
 };
 
-template <typename X, typename I1, typename I2>
-using DDim = std::conditional_t<
-        std::is_same_v<X, I1>,
-        DDimGPS1<X>,
-        std::conditional_t<std::is_same_v<X, I2>, DDimGPS2<X>, X>>;
-
 #if defined(BC_PERIODIC)
 template <typename DDim1, typename DDim2>
 using evaluator_type = Evaluator2D::
@@ -110,10 +104,6 @@ using DVect = ddc::DiscreteVector<DDimX...>;
 template <typename... X>
 using Coord = ddc::Coordinate<X...>;
 
-// Extract batch dimensions from DDim (remove dimension of interest). Usefull
-template <typename I1, typename I2, typename... X>
-using BatchDims = ddc::type_seq_remove_t<ddc::detail::TypeSeq<X...>, ddc::detail::TypeSeq<I1, I2>>;
-
 // Templated function giving first coordinate of the mesh in given dimension.
 template <typename X>
 KOKKOS_FUNCTION Coord<X> x0()
@@ -159,35 +149,35 @@ void InterestDimInitializer(std::size_t const ncells)
 
 // Checks that when evaluating the spline at interpolation points one
 // recovers values that were used to build the spline
-template <typename ExecSpace, typename MemorySpace, typename I1, typename I2, typename... X>
+template <
+        typename ExecSpace,
+        typename MemorySpace,
+        typename DDimI1,
+        typename DDimI2,
+        typename... DDims>
 void ExtrapolationRuleSplineTest()
 {
+    using I1 = typename DDimI1::continuous_dimension_type;
+    using I2 = typename DDimI2::continuous_dimension_type;
+
     // Instantiate execution spaces and initialize spaces
     ExecSpace const exec_space;
     std::size_t const ncells = 10;
-    InterestDimInitializer<DDim<I1, I1, I2>>(ncells);
-    ddc::init_discrete_space<DDim<I1, I1, I2>>(
-            GrevillePoints1<BSplines<I1>>::template get_sampling<DDim<I1, I1, I2>>());
-    InterestDimInitializer<DDim<I2, I1, I2>>(ncells);
-    ddc::init_discrete_space<DDim<I2, I1, I2>>(
-            GrevillePoints2<BSplines<I2>>::template get_sampling<DDim<I2, I1, I2>>());
+    InterestDimInitializer<DDimI1>(ncells);
+    ddc::init_discrete_space<DDimI1>(
+            GrevillePoints1<BSplines<I1>>::template get_sampling<DDimI1>());
+    InterestDimInitializer<DDimI2>(ncells);
+    ddc::init_discrete_space<DDimI2>(
+            GrevillePoints2<BSplines<I2>>::template get_sampling<DDimI2>());
 
     // Create the values domain (mesh)
-    ddc::DiscreteDomain<DDim<I1, I1, I2>, DDim<I2, I1, I2>> const interpolation_domain(
-            GrevillePoints1<BSplines<I1>>::template get_domain<DDim<I1, I1, I2>>(),
-            GrevillePoints2<BSplines<I2>>::template get_domain<DDim<I2, I1, I2>>());
-    // If we remove auto using the constructor syntax, nvcc does not compile
-    auto const dom_vals_tmp = ddc::DiscreteDomain<DDim<X, void, void>...>(
-            ddc::DiscreteDomain<DDim<
-                    X,
-                    void,
-                    void>>(DElem<DDim<X, void, void>>(0), DVect<DDim<X, void, void>>(ncells))...);
-    ddc::DiscreteDomain<DDim<X, I1, I2>...> const dom_vals
-            = ddc::replace_dim_of<DDim<I1, void, void>, DDim<I1, I1, I2>>(
-                    ddc::replace_dim_of<
-                            DDim<I2, void, void>,
-                            DDim<I2, I1, I2>>(dom_vals_tmp, interpolation_domain),
-                    interpolation_domain);
+    ddc::DiscreteDomain<DDimI1, DDimI2> const interpolation_domain(
+            GrevillePoints1<BSplines<I1>>::template get_domain<DDimI1>(),
+            GrevillePoints2<BSplines<I2>>::template get_domain<DDimI2>());
+    // The following line creates a discrete domain over all dimensions (DDims...) except DDimI1 and DDimI2.
+    auto const dom_vals_tmp = ddc::remove_dims_of_t<ddc::DiscreteDomain<DDims...>, DDimI1, DDimI2>(
+            ddc::DiscreteDomain<DDims>(DElem<DDims>(0), DVect<DDims>(ncells))...);
+    ddc::DiscreteDomain<DDims...> const dom_vals(dom_vals_tmp, interpolation_domain);
 
     // Create a SplineBuilder over BSplines<I> and batched along other dimensions using some boundary conditions
     ddc::SplineBuilder2D<
@@ -195,24 +185,24 @@ void ExtrapolationRuleSplineTest()
             MemorySpace,
             BSplines<I1>,
             BSplines<I2>,
-            DDim<I1, I1, I2>,
-            DDim<I2, I1, I2>,
+            DDimI1,
+            DDimI2,
             s_bcl1,
             s_bcr1,
             s_bcl2,
             s_bcr2,
             ddc::SplineSolver::GINKGO,
-            DDim<X, I1, I2>...> const spline_builder(dom_vals);
+            DDims...> const spline_builder(dom_vals);
 
     // Compute usefull domains (dom_interpolation, dom_batch, dom_bsplines and dom_spline)
-    ddc::DiscreteDomain<DDim<I1, I1, I2>, DDim<I2, I1, I2>> const dom_interpolation
+    ddc::DiscreteDomain<DDimI1, DDimI2> const dom_interpolation
             = spline_builder.interpolation_domain();
     auto const dom_spline = spline_builder.batched_spline_domain();
 
     // Allocate and fill a chunk containing values to be passed as input to spline_builder. Those are values of cosine along interest dimension duplicated along batch dimensions
     ddc::Chunk vals_1d_host_alloc(dom_interpolation, ddc::HostAllocator<double>());
     ddc::ChunkSpan const vals_1d_host = vals_1d_host_alloc.span_view();
-    evaluator_type<DDim<I1, I1, I2>, DDim<I2, I1, I2>> const evaluator(dom_interpolation);
+    evaluator_type<DDimI1, DDimI2> const evaluator(dom_interpolation);
     evaluator(vals_1d_host);
     auto vals_1d_alloc = ddc::create_mirror_view_and_copy(exec_space, vals_1d_host);
     ddc::ChunkSpan const vals_1d = vals_1d_alloc.span_view();
@@ -222,8 +212,8 @@ void ExtrapolationRuleSplineTest()
     ddc::parallel_for_each(
             exec_space,
             vals.domain(),
-            KOKKOS_LAMBDA(DElem<DDim<X, I1, I2>...> const e) {
-                vals(e) = vals_1d(DElem<DDim<I1, I1, I2>, DDim<I2, I1, I2>>(e));
+            KOKKOS_LAMBDA(DElem<DDims...> const e) {
+                vals(e) = vals_1d(DElem<DDimI1, DDimI2>(e));
             });
 
     // Instantiate chunk of spline coefs to receive output of spline_builder
@@ -272,13 +262,13 @@ void ExtrapolationRuleSplineTest()
             MemorySpace,
             BSplines<I1>,
             BSplines<I2>,
-            DDim<I1, I1, I2>,
-            DDim<I2, I1, I2>,
+            DDimI1,
+            DDimI2,
             extrapolation_rule_dim_1_type,
             extrapolation_rule_dim_1_type,
             extrapolation_rule_dim_2_type,
             extrapolation_rule_dim_2_type,
-            DDim<X, I1, I2>...> const
+            DDims...> const
             spline_evaluator_batched(
                     extrapolation_rule_left_dim_1,
                     extrapolation_rule_right_dim_1,
@@ -295,9 +285,8 @@ void ExtrapolationRuleSplineTest()
     ddc::parallel_for_each(
             exec_space,
             coords_eval.domain(),
-            KOKKOS_LAMBDA(DElem<DDim<X, I1, I2>...> const e) {
-                coords_eval(e)
-                        = ddc::coordinate(DElem<DDim<I1, I1, I2>, DDim<I2, I1, I2>>(e)) + displ;
+            KOKKOS_LAMBDA(DElem<DDims...> const e) {
+                coords_eval(e) = ddc::coordinate(DElem<DDimI1, DDimI2>(e)) + displ;
             });
 
 
@@ -314,35 +303,28 @@ void ExtrapolationRuleSplineTest()
             spline_eval.domain(),
             0.,
             ddc::reducer::max<double>(),
-            KOKKOS_LAMBDA(DElem<DDim<X, I1, I2>...> const e) {
+            KOKKOS_LAMBDA(DElem<DDims...> const e) {
 #if defined(ER_NULL)
                 return Kokkos::abs(spline_eval(e));
 #elif defined(ER_CONSTANT)
                 typename decltype(ddc::remove_dims_of(
                         vals.domain(),
-                        vals.template domain<DDim<I1, I1, I2>>()))::discrete_element_type const
+                        vals.template domain<DDimI1>()))::discrete_element_type const
                         e_without_interest(e);
-                double tmp;
-                if (Coord<I2>(coords_eval(e)) > xN<I2>()) {
 #if defined(BC_PERIODIC)
-                    tmp = vals(ddc::DiscreteElement<DDim<X, I1, I2>...>(
-                            vals.template domain<DDim<I1, I1, I2>>().back(),
-                            e_without_interest));
+                double const tmp = vals(vals.template domain<DDimI1>().back(), e_without_interest);
 #else
+                double tmp;
+                if (Coord<I2>(coords_eval(e)) > xN<I2>()) {
                     typename decltype(ddc::remove_dims_of(
                             vals.domain(),
-                            vals.template domain<DDim<I1, I1, I2>, DDim<I2, I1, I2>>()))::
-                            discrete_element_type const e_batch(e);
-                    tmp = vals(ddc::DiscreteElement<DDim<X, I1, I2>...>(
-                            vals.template domain<DDim<I1, I1, I2>>().back(),
-                            vals.template domain<DDim<I2, I1, I2>>().back(),
-                            e_batch));
-#endif
+                            vals.template domain<DDimI1, DDimI2>()))::discrete_element_type const
+                            e_batch(e);
+                    tmp = vals(vals.template domain<DDimI1, DDimI2>().back(), e_batch);
                 } else {
-                    tmp = vals(ddc::DiscreteElement<DDim<X, I1, I2>...>(
-                            vals.template domain<DDim<I1, I1, I2>>().back(),
-                            e_without_interest));
+                    tmp = vals(vals.template domain<DDimI1>().back(), e_without_interest);
                 }
+#endif
                 return Kokkos::abs(spline_eval(e) - tmp);
 #endif
             });
@@ -377,10 +359,10 @@ TEST(SUFFIX(ExtrapolationRuleSplineHost), 2DXY)
     ExtrapolationRuleSplineTest<
             Kokkos::DefaultHostExecutionSpace,
             Kokkos::DefaultHostExecutionSpace::memory_space,
-            DimX,
-            DimY,
-            DimX,
-            DimY>();
+            DDimGPS1<DimX>,
+            DDimGPS2<DimY>,
+            DDimGPS1<DimX>,
+            DDimGPS2<DimY>>();
 }
 
 TEST(SUFFIX(ExtrapolationRuleSplineDevice), 2DXY)
@@ -388,8 +370,8 @@ TEST(SUFFIX(ExtrapolationRuleSplineDevice), 2DXY)
     ExtrapolationRuleSplineTest<
             Kokkos::DefaultExecutionSpace,
             Kokkos::DefaultExecutionSpace::memory_space,
-            DimX,
-            DimY,
-            DimX,
-            DimY>();
+            DDimGPS1<DimX>,
+            DDimGPS2<DimY>,
+            DDimGPS1<DimX>,
+            DDimGPS2<DimY>>();
 }