diff --git a/examples/free_scalar_field/2d_free_scalar_field.cpp b/examples/free_scalar_field/2d_free_scalar_field.cpp
index e3c1828..3552145 100644
--- a/examples/free_scalar_field/2d_free_scalar_field.cpp
+++ b/examples/free_scalar_field/2d_free_scalar_field.cpp
@@ -273,6 +273,11 @@ using BetaLow = sil::tensor::Covariant<Beta>;
 
 using DummyIndex = sil::tensor::Covariant<sil::tensor::ScalarIndex>;
 
+using XDualizer = sil::mesher::HalfShiftDualizer<X>;
+using YDualizer = sil::mesher::HalfShiftDualizer<Y>;
+using DDimXDual = sil::mesher::dual_discrete_dimension_t<XDualizer, DDimX>;
+using DDimYDual = sil::mesher::dual_discrete_dimension_t<YDualizer, DDimY>;
+
 int main(int argc, char** argv)
 {
     // Initialize PDI, Kokkos and DDC
@@ -306,6 +311,10 @@ int main(int argc, char** argv)
             ddc::select<Y>(upper_bounds),
             ddc::select<DDimY>(nb_cells)));
     ddc::DiscreteDomain<DDimX, DDimY> mesh_xy(x_dom, y_dom);
+    XDualizer const x_dualizer;
+    YDualizer const y_dualizer;
+    ddc::DiscreteDomain<DDimXDual, DDimY> x_face_dom = x_dualizer(mesh_xy);
+    ddc::DiscreteDomain<DDimX, DDimYDual> y_face_dom = y_dualizer(mesh_xy);
 
     assert(static_cast<std::size_t>(mesh_xy.template extent<DDimX>())
            == static_cast<std::size_t>(mesh_xy.template extent<DDimY>()));
@@ -369,9 +378,11 @@ int main(int argc, char** argv)
     double const y_1 = -0.3;
     double const sigma = .5;
 
-    double const v = 10.;
-    double const k = 1.;
-    double const mass = 1.;
+    double const v = 0.5;
+    double const mass = 100;
+    double const linear_coupling_constant = 0;
+    double const quartic_coupling_constant = 0;
+    double const k = mass * v / std::sqrt(1. - v * v);
 
     ddc::parallel_for_each(
             potential.domain(),
@@ -383,10 +394,9 @@ int main(int argc, char** argv)
                                           * std::exp(
                                                   -((x - x_0) * (x - x_0) + (y - y_0) * (y - y_0))
                                                   / 2. / sigma / sigma)
-                                  + 30. * std::sin(k * (x - x_1))
-                                            * std::exp(
-                                                    -((x - x_1) * (x - x_1) + (y - y_1) * (y - y_1))
-                                                    / 2. / sigma / sigma);
+                                  + std::exp(
+                                          -((x - x_1) * (x - x_1) + (y - y_1) * (y - y_1)) / 2.
+                                          / sigma / sigma);
             });
     auto potential_host
             = ddc::create_mirror_view_and_copy(Kokkos::DefaultHostExecutionSpace(), potential);
@@ -395,17 +405,46 @@ int main(int argc, char** argv)
     ddc::Chunk half_step_potential_alloc(potential_dom, ddc::DeviceAllocator<double>());
     sil::tensor::Tensor half_step_potential(half_step_potential_alloc);
 
-    // Potential gradient
-    [[maybe_unused]] sil::tensor::TensorAccessor<AlphaLow> potential_grad_accessor;
+    // Staggered potential gradients
+    [[maybe_unused]] sil::tensor::TensorAccessor<DummyIndex> scalar_accessor;
+    ddc::DiscreteDomain<DDimXDual, DDimY, DummyIndex>
+            potential_grad_x_dom(x_face_dom, scalar_accessor.domain());
+    ddc::Chunk potential_grad_x_alloc(potential_grad_x_dom, ddc::DeviceAllocator<double>());
+    sil::tensor::Tensor potential_grad_x(potential_grad_x_alloc);
+
+    ddc::DiscreteDomain<DDimX, DDimYDual, DummyIndex>
+            potential_grad_y_dom(y_face_dom, scalar_accessor.domain());
+    ddc::Chunk potential_grad_y_alloc(potential_grad_y_dom, ddc::DeviceAllocator<double>());
+    sil::tensor::Tensor potential_grad_y(potential_grad_y_alloc);
+    auto potential_grad = sil::exterior::make_tensor_form(
+            sil::exterior::component<X>(potential_grad_x),
+            sil::exterior::component<Y>(potential_grad_y));
+
+    // Staggered spatial moments
+    ddc::DiscreteDomain<DDimXDual, DDimY, DummyIndex>
+            spatial_moment_x_dom(x_face_dom, scalar_accessor.domain());
+    ddc::Chunk spatial_moment_x_alloc(spatial_moment_x_dom, ddc::DeviceAllocator<double>());
+    sil::tensor::Tensor spatial_moment_x(spatial_moment_x_alloc);
+
+    ddc::DiscreteDomain<DDimX, DDimYDual, DummyIndex>
+            spatial_moment_y_dom(y_face_dom, scalar_accessor.domain());
+    ddc::Chunk spatial_moment_y_alloc(spatial_moment_y_dom, ddc::DeviceAllocator<double>());
+    sil::tensor::Tensor spatial_moment_y(spatial_moment_y_alloc);
+    auto spatial_moment = sil::exterior::make_tensor_form(
+            sil::exterior::component<X>(spatial_moment_x),
+            sil::exterior::component<Y>(spatial_moment_y));
+
+    auto spatial_moment_x_host = ddc::
+            create_mirror_view_and_copy(Kokkos::DefaultHostExecutionSpace(), spatial_moment_x);
+    auto spatial_moment_y_host = ddc::
+            create_mirror_view_and_copy(Kokkos::DefaultHostExecutionSpace(), spatial_moment_y);
+
+    // Reconstructed spatial moments on nodes, only for export and diagnostics
+    [[maybe_unused]] sil::tensor::TensorAccessor<AlphaLow> spatial_moments_accessor;
     ddc::DiscreteDomain<DDimX, DDimY, AlphaLow>
-            potential_grad_dom(mesh_xy, potential_grad_accessor.domain());
-    ddc::Chunk potential_grad_alloc(potential_grad_dom, ddc::DeviceAllocator<double>());
-    sil::tensor::Tensor potential_grad(potential_grad_alloc);
-
-    // Spatial moments
-    auto& spatial_moments
-            = potential_grad; // We can perform the computations inplace so spatial_moments is just an alias of potential_grad
-
+            spatial_moments_dom(mesh_xy, spatial_moments_accessor.domain());
+    ddc::Chunk spatial_moments_alloc(spatial_moments_dom, ddc::DeviceAllocator<double>());
+    sil::tensor::Tensor spatial_moments(spatial_moments_alloc);
     auto spatial_moments_host = ddc::
             create_mirror_view_and_copy(Kokkos::DefaultHostExecutionSpace(), spatial_moments);
 
@@ -425,15 +464,14 @@ int main(int argc, char** argv)
     ddc::Chunk temporal_moment_alloc(temporal_moment_dom, ddc::DeviceAllocator<double>());
     sil::tensor::Tensor temporal_moment(temporal_moment_alloc);
 
+    double const omega = std::sqrt(k * k + mass * mass);
     ddc::parallel_for_each(
             potential.domain(),
             KOKKOS_LAMBDA(ddc::DiscreteElement<DDimX, DDimY, DummyIndex> elem) {
                 double const x = ddc::coordinate(ddc::DiscreteElement<DDimX>(elem));
                 double const y = ddc::coordinate(ddc::DiscreteElement<DDimY>(elem)) - y_0;
-                // v*dphi/dx of the left wave packet only to get a pure kick along x toward the immobile right one
-                temporal_moment(elem) = -v
-                                        * (-k * std::cos(k * (x - x_0))
-                                           + (x - x_0) / sigma / sigma * std::sin(k * (x - x_0)))
+                temporal_moment(elem) = (-omega * std::cos(k * (x - x_0))
+                                         + v * (x - x_0) / sigma / sigma * std::sin(k * (x - x_0)))
                                         * std::exp(
                                                 -((x - x_0) * (x - x_0) + (y - y_0) * (y - y_0))
                                                 / 2. / sigma / sigma);
@@ -491,36 +529,56 @@ int main(int argc, char** argv)
                 KOKKOS_LAMBDA(ddc::DiscreteElement<DDimX, DDimY, DummyIndex> elem) {
                     half_step_potential(elem)
                             = potential(elem)
-                              - FreeScalarFieldHamiltonian(mass).dH_dpi0(temporal_moment(elem)) * dt
-                                        / 2.;
+                              - FreeScalarFieldHamiltonian(
+                                        mass,
+                                        quartic_coupling_constant,
+                                        linear_coupling_constant)
+                                                .dH_dpi0(temporal_moment(elem))
+                                        * dt / 2.;
                 });
 
-        // Compute the potential gradient
-        sil::exterior::deriv<
-                AlphaLow,
-                DummyIndex>(Kokkos::DefaultExecutionSpace(), potential_grad, half_step_potential);
+        // Compute the staggered potential gradients
+        sil::exterior::deriv(Kokkos::DefaultExecutionSpace(), potential_grad, half_step_potential);
 
         // Compute the spatial moments pi_\alpha by solving dphi/dx^\alpha = dH/dpi_\alpha
         ddc::parallel_for_each(
                 Kokkos::DefaultExecutionSpace(),
-                mesh_xy,
-                KOKKOS_LAMBDA(ddc::DiscreteElement<DDimX, DDimY> elem) {
-                    spatial_moments(elem, ddc::DiscreteElement<AlphaLow>(0))
-                            = FreeScalarFieldHamiltonian(mass).pi1(
-                                    potential_grad(elem, ddc::DiscreteElement<AlphaLow>(0)));
-                    spatial_moments(elem, ddc::DiscreteElement<AlphaLow>(1))
-                            = FreeScalarFieldHamiltonian(mass).pi2(
-                                    potential_grad(elem, ddc::DiscreteElement<AlphaLow>(1)));
+                x_face_dom,
+                KOKKOS_LAMBDA(ddc::DiscreteElement<DDimXDual, DDimY> elem) {
+                    spatial_moment_x(elem, ddc::DiscreteElement<DummyIndex>(0))
+                            = FreeScalarFieldHamiltonian(
+                                      mass,
+                                      quartic_coupling_constant,
+                                      linear_coupling_constant)
+                                      .pi1(potential_grad_x(
+                                              elem,
+                                              ddc::DiscreteElement<DummyIndex>(0)));
+                });
+        ddc::parallel_for_each(
+                Kokkos::DefaultExecutionSpace(),
+                y_face_dom,
+                KOKKOS_LAMBDA(ddc::DiscreteElement<DDimX, DDimYDual> elem) {
+                    spatial_moment_y(elem, ddc::DiscreteElement<DummyIndex>(0))
+                            = FreeScalarFieldHamiltonian(
+                                      mass,
+                                      quartic_coupling_constant,
+                                      linear_coupling_constant)
+                                      .pi2(potential_grad_y(
+                                              elem,
+                                              ddc::DiscreteElement<DummyIndex>(0)));
                 });
-        if (i % nb_iter_between_exports == 0) {
-            ddc::parallel_deepcopy(spatial_moments_host, spatial_moments);
-        }
 
         // Compute the divergence dpi_\alpha/dx^\alpha of the spatial moments, which is the codifferential \delta pi of the spatial moments
-        sil::exterior::codifferential<MetricIndex, AlphaLow, AlphaLow>(
+        ddc::parallel_for_each(
+                Kokkos::DefaultExecutionSpace(),
+                spatial_moments_div.domain(),
+                KOKKOS_LAMBDA(ddc::DiscreteElement<DDimX, DDimY, DummyIndex> elem) {
+                    spatial_moments_div(elem) = 0.;
+                });
+        sil::exterior::codifferential<MetricIndex>(
                 Kokkos::DefaultExecutionSpace(),
                 spatial_moments_div,
-                spatial_moments,
+                spatial_moment,
                 inv_metric);
 
         // Compute dpi_0/dx^0 by solving - dpi_0/dx^0 + dpi_\alpha/dx^\alpha = -dH/dphi and advect pi_0 by a time step dx^0. Also Then, perform the second phi half-advection by solving dphi/dx^0 = -dH/dpi_0
@@ -538,17 +596,29 @@ int main(int argc, char** argv)
                     // Advect temporal moment by half-step, this is what is needed to perform the whole-step potential advection
                     const double half_step_temporal_moment_
                             = temporal_moment_
-                              + (FreeScalarFieldHamiltonian(mass).dH_dphi(half_step_potential_)
-                                 + spatial_moments_div_)
+                              + (spatial_moments_div_
+                                 - FreeScalarFieldHamiltonian(
+                                           mass,
+                                           quartic_coupling_constant,
+                                           linear_coupling_constant)
+                                           .dH_dphi(half_step_potential_))
                                         * dt / 2;
 
                     // Whole-step advection of field state
                     potential(elem)
-                            += FreeScalarFieldHamiltonian(mass).dH_dpi0(half_step_temporal_moment_)
-                               * dt;
+                            -= FreeScalarFieldHamiltonian(
+                                       mass,
+                                       quartic_coupling_constant,
+                                       linear_coupling_constant)
+                                       .dH_dpi0(half_step_temporal_moment_)
+                                       * dt;
                     temporal_moment(elem)
-                            += (FreeScalarFieldHamiltonian(mass).dH_dphi(half_step_potential_)
-                                + spatial_moments_div_)
+                            += (spatial_moments_div_
+                                - FreeScalarFieldHamiltonian(
+                                          mass,
+                                          quartic_coupling_constant,
+                                          linear_coupling_constant)
+                                          .dH_dphi(half_step_potential_))
                                * dt;
                 });
 
@@ -556,6 +626,59 @@ int main(int argc, char** argv)
             ddc::parallel_deepcopy(temporal_moment_host, temporal_moment);
             ddc::parallel_deepcopy(potential_host, potential);
             ddc::parallel_deepcopy(h_spatial_moments_div, spatial_moments_div);
+            ddc::parallel_deepcopy(spatial_moment_x_host, spatial_moment_x);
+            ddc::parallel_deepcopy(spatial_moment_y_host, spatial_moment_y);
+
+            {
+                ddc::DiscreteElement<DDimX, DDimY> const mesh_front = mesh_xy.front();
+                for (std::size_t ix = 0; ix < mesh_xy.template extent<DDimX>(); ++ix) {
+                    for (std::size_t iy = 0; iy < mesh_xy.template extent<DDimY>(); ++iy) {
+                        ddc::DiscreteElement<DDimX, DDimY> const elem
+                                = mesh_front + ddc::DiscreteVector<DDimX, DDimY>(ix, iy);
+                        spatial_moments_host(elem, ddc::DiscreteElement<AlphaLow>(0)) = 0.;
+                        spatial_moments_host(elem, ddc::DiscreteElement<AlphaLow>(1)) = 0.;
+                    }
+                }
+            }
+            {
+                ddc::DiscreteElement<DDimXDual, DDimY> const x_face_front = x_face_dom.front();
+                for (std::size_t ix = 0; ix < x_face_dom.template extent<DDimXDual>(); ++ix) {
+                    for (std::size_t iy = 0; iy < x_face_dom.template extent<DDimY>(); ++iy) {
+                        ddc::DiscreteElement<DDimXDual, DDimY> const elem
+                                = x_face_front + ddc::DiscreteVector<DDimXDual, DDimY>(ix, iy);
+                        ddc::DiscreteElement<DDimX, DDimY> const left_node
+                                = x_dualizer.primal(elem);
+                        double const value
+                                = spatial_moment_x_host(elem, ddc::DiscreteElement<DummyIndex>(0));
+                        spatial_moments_host(left_node, ddc::DiscreteElement<AlphaLow>(0))
+                                += value / 2.;
+                        spatial_moments_host(
+                                left_node + ddc::DiscreteVector<DDimX, DDimY>(1, 0),
+                                ddc::DiscreteElement<AlphaLow>(0))
+                                += value / 2.;
+                    }
+                }
+            }
+            {
+                ddc::DiscreteElement<DDimX, DDimYDual> const y_face_front = y_face_dom.front();
+                for (std::size_t ix = 0; ix < y_face_dom.template extent<DDimX>(); ++ix) {
+                    for (std::size_t iy = 0; iy < y_face_dom.template extent<DDimYDual>(); ++iy) {
+                        ddc::DiscreteElement<DDimX, DDimYDual> const elem
+                                = y_face_front + ddc::DiscreteVector<DDimX, DDimYDual>(ix, iy);
+                        ddc::DiscreteElement<DDimX, DDimY> const lower_node
+                                = y_dualizer.primal(elem);
+                        double const value
+                                = spatial_moment_y_host(elem, ddc::DiscreteElement<DummyIndex>(0));
+                        spatial_moments_host(lower_node, ddc::DiscreteElement<AlphaLow>(1))
+                                += value / 2.;
+                        spatial_moments_host(
+                                lower_node + ddc::DiscreteVector<DDimX, DDimY>(0, 1),
+                                ddc::DiscreteElement<AlphaLow>(1))
+                                += value / 2.;
+                    }
+                }
+            }
+            ddc::parallel_deepcopy(spatial_moments, spatial_moments_host);
 
             ddc::parallel_for_each(
                     Kokkos::DefaultExecutionSpace(),
@@ -566,7 +689,10 @@ int main(int argc, char** argv)
                                    spatial_moments(elem, ddc::DiscreteElement<AlphaLow>(0)),
                                    spatial_moments(elem, ddc::DiscreteElement<AlphaLow>(1))};
                         hamiltonian(elem)
-                                = FreeScalarFieldHamiltonian(mass)
+                                = FreeScalarFieldHamiltonian(
+                                          mass,
+                                          quartic_coupling_constant,
+                                          linear_coupling_constant)
                                           .H(potential(elem, ddc::DiscreteElement<DummyIndex>()),
                                              pi);
                     });
diff --git a/examples/free_scalar_field/free_scalar_field_hamiltonian_functor_generator.py b/examples/free_scalar_field/free_scalar_field_hamiltonian_functor_generator.py
index 772e0e4..0877112 100644
--- a/examples/free_scalar_field/free_scalar_field_hamiltonian_functor_generator.py
+++ b/examples/free_scalar_field/free_scalar_field_hamiltonian_functor_generator.py
@@ -14,17 +14,20 @@
 output_path = Path(sys.argv[2])  # Path of the output file
 output_path.parent.mkdir(parents=True, exist_ok=True)
 
-mass = symbols("mass")
+mass, quartic_coupling_constant, linear_coupling_constant = symbols(
+    "mass quartic_coupling_constant linear_coupling_constant"
+)
 phi = symbols("phi")
 pi = symbols(f"pi0:{N}")
 
 # Minkowski signature (+, -, -, ..., -)
 metric_sign = [-1] + [1] * (N - 1)
 
-# H = 1/2 m^2 phi^2 + 1/2 (-p0^2 + p1^2 - ... - p{N-1}^2)
+# H = linear_coupling_constant * phi + 1/2 m^2 phi^2 + lambda/4 phi^4
+#   + 1/2 (-p0^2 + p1^2 - ... - p{N-1}^2)
 hamiltonian = 0.5 * (
     mass**2 * phi**2 + sum(metric_sign[i] * pi[i] ** 2 for i in range(N))
-)
+) + quartic_coupling_constant * phi**4 / 4 + linear_coupling_constant * phi
 
 # [dH/dphi, dH/dpi0, dH/dpi1, ...]
 hamiltonian_diff = Matrix(
@@ -69,8 +72,17 @@ def preprocess_cxx(expr: str) -> str:
     static constexpr std::size_t N = {N};
 
     const double mass;
-
-    constexpr FreeScalarFieldHamiltonian(double mass_) : mass(mass_) {{}}
+    const double quartic_coupling_constant;
+    const double linear_coupling_constant;
+
+    constexpr FreeScalarFieldHamiltonian(
+            double mass_,
+            double quartic_coupling_constant_,
+            double linear_coupling_constant_)
+        : mass(mass_)
+        , quartic_coupling_constant(quartic_coupling_constant_)
+        , linear_coupling_constant(linear_coupling_constant_)
+    {{}}
 
     constexpr double H(double phi, const std::span<const double, N>& pi) const
     {{
diff --git a/src/similie/exterior/coboundary.hpp b/src/similie/exterior/coboundary.hpp
index b817ec1..9a05e0b 100644
--- a/src/similie/exterior/coboundary.hpp
+++ b/src/similie/exterior/coboundary.hpp
@@ -12,6 +12,7 @@
 #include <similie/misc/portable_stl.hpp>
 #include <similie/misc/specialization.hpp>
 #include <similie/misc/type_seq_conversion.hpp>
+#include <similie/mesher/dualizer.hpp>
 #include <similie/tensor/antisymmetric_tensor.hpp>
 #include <similie/tensor/dummy_index.hpp>
 #include <similie/tensor/tensor_impl.hpp>
@@ -20,6 +21,7 @@
 
 #include "cochain.hpp"
 #include "cosimplex.hpp"
+#include "form.hpp"
 
 
 namespace sil {
@@ -334,6 +336,75 @@ coboundary_tensor_t<TagToAddToCochain, CochainTag, TensorType> deriv(
             TensorType>(exec_space, coboundary_tensor, tensor);
 }
 
+template <
+        class TagToAddToCochain,
+        tensor::TensorIndex CochainTag,
+        misc::Specialization<tensor::Tensor> OutTensorType,
+        misc::Specialization<tensor::Tensor> TensorType,
+        class Dualizer,
+        class ExecSpace>
+OutTensorType deriv(
+        ExecSpace const& exec_space,
+        OutTensorType out_tensor,
+        TensorType tensor,
+        Dualizer const& dualizer)
+{
+    using out_index_t = ddc::type_seq_element_t<
+            0,
+            ddc::to_type_seq_t<typename OutTensorType::indices_domain_t>>;
+    using in_index_t = ddc::type_seq_element_t<
+            0,
+            ddc::to_type_seq_t<typename TensorType::indices_domain_t>>;
+    using in_d_dim_t = mesher::detail::
+            discrete_dimension_for_t<TagToAddToCochain, typename TensorType::non_indices_domain_t>;
+    using out_elem_t = typename OutTensorType::non_indices_domain_t::discrete_element_type;
+    using in_elem_t = typename TensorType::non_indices_domain_t::discrete_element_type;
+    static_assert(out_index_t::rank() == 0);
+    static_assert(in_index_t::rank() == 0);
+    static_assert(!std::is_void_v<in_d_dim_t>);
+
+    ddc::parallel_for_each(
+            "similie_compute_centered_dualized_deriv",
+            exec_space,
+            out_tensor.non_indices_domain(),
+            KOKKOS_LAMBDA(out_elem_t out_elem) {
+                in_elem_t const in_left = dualizer.primal(out_elem);
+                in_elem_t const in_right = in_left + ddc::DiscreteVector<in_d_dim_t>(1);
+                double const dx = static_cast<double>(
+                        ddc::coordinate(ddc::DiscreteElement<in_d_dim_t>(in_right))
+                        - ddc::coordinate(ddc::DiscreteElement<in_d_dim_t>(in_left)));
+                out_tensor.mem(out_elem, ddc::DiscreteElement<out_index_t>(0))
+                        = (tensor.get(in_right, ddc::DiscreteElement<in_index_t>(0))
+                           - tensor.get(in_left, ddc::DiscreteElement<in_index_t>(0)))
+                          / dx;
+            });
+    return out_tensor;
+}
+
+template <
+        class SupportTag,
+        class... Components,
+        misc::Specialization<tensor::Tensor> TensorType,
+        class ExecSpace>
+TensorForm<SupportTag, Components...> deriv(
+        ExecSpace const& exec_space,
+        TensorForm<SupportTag, Components...> out_form,
+        TensorType tensor)
+{
+    using in_index_t = ddc::type_seq_element_t<
+            0,
+            ddc::to_type_seq_t<typename TensorType::indices_domain_t>>;
+    static_assert(in_index_t::rank() == 0);
+
+    (...,
+     sil::exterior::deriv<typename Components::tag, in_index_t>(
+             exec_space,
+             out_form.template component<typename Components::tag>(),
+             tensor,
+             mesher::HalfShiftDualizer<typename Components::tag> {}));
+    return out_form;
+}
+
 } // namespace exterior
 
 } // namespace sil
diff --git a/src/similie/exterior/codifferential.hpp b/src/similie/exterior/codifferential.hpp
index 49a0a2e..3ce2990 100644
--- a/src/similie/exterior/codifferential.hpp
+++ b/src/similie/exterior/codifferential.hpp
@@ -6,6 +6,7 @@
 #include <ddc/ddc.hpp>
 
 #include <similie/exterior/hodge_star.hpp>
+#include <similie/mesher/dualizer.hpp>
 #include <similie/misc/macros.hpp>
 #include <similie/misc/specialization.hpp>
 #include <similie/tensor/character.hpp>
@@ -15,6 +16,7 @@
 
 #include "cochain.hpp"
 #include "cosimplex.hpp"
+#include "form.hpp"
 
 
 namespace sil {
@@ -156,6 +158,125 @@ struct CodifferentialDummyIndexSeq<EndId, T>
     using type = typename CodifferentialDummyIndexSeq_<std::make_index_sequence<EndId>, T>::type;
 };
 
+template <class DDim>
+constexpr bool is_dual_discrete_dimension_v
+        = !std::is_same_v<DDim, mesher::detail::primal_discrete_dimension_t<DDim>>;
+
+template <class DDim, class Domain, class Elem>
+KOKKOS_FUNCTION bool is_on_primal_boundary(Domain const& domain, Elem const& elem)
+{
+    if constexpr (is_dual_discrete_dimension_v<DDim>) {
+        return false;
+    } else {
+        ddc::DiscreteDomain<DDim> const dim_dom(domain);
+        ddc::DiscreteElement<DDim> const dim_elem(elem);
+        return dim_elem.uid() == dim_dom.front().uid() || dim_elem.uid() == dim_dom.back().uid();
+    }
+}
+
+template <
+        class MetricComponentIndex,
+        class TagToRemoveFromCochain,
+        misc::Specialization<tensor::Tensor> MetricType,
+        class Elem>
+KOKKOS_FUNCTION double codifferential_metric_factor(MetricType const& inv_metric, Elem const& elem)
+{
+    if constexpr (misc::Specialization<MetricComponentIndex, tensor::TensorIdentityIndex>) {
+        return 1.;
+    } else {
+        return inv_metric(
+                elem,
+                inv_metric.accessor().template access_element<
+                        TagToRemoveFromCochain,
+                        TagToRemoveFromCochain>());
+    }
+}
+
+template <
+        class AxisTag,
+        misc::Specialization<tensor::Tensor> OutTensorType,
+        misc::Specialization<tensor::Tensor> TensorType>
+KOKKOS_FUNCTION double centered_form_derivative(
+        typename OutTensorType::non_indices_domain_t::discrete_element_type out_elem,
+        TensorType tensor)
+{
+    using out_d_dim_t = mesher::detail::
+            discrete_dimension_for_t<AxisTag, typename OutTensorType::non_indices_domain_t>;
+    using comp_d_dim_t = mesher::detail::
+            discrete_dimension_for_t<AxisTag, typename TensorType::non_indices_domain_t>;
+    using comp_elem_t = typename TensorType::non_indices_domain_t::discrete_element_type;
+    using comp_index_t = ddc::type_seq_element_t<
+            0,
+            ddc::to_type_seq_t<typename TensorType::indices_domain_t>>;
+    static_assert(!std::is_void_v<out_d_dim_t>);
+    static_assert(!std::is_void_v<comp_d_dim_t>);
+    static_assert(comp_index_t::rank() == 0);
+
+    comp_elem_t const aligned_elem = mesher::detail::dualizer_remap_element<comp_elem_t>(out_elem);
+    if constexpr (is_dual_discrete_dimension_v<out_d_dim_t>) {
+        auto const left = aligned_elem;
+        auto const right = left + ddc::DiscreteVector<comp_d_dim_t>(1);
+        double const dx = static_cast<double>(
+                ddc::coordinate(ddc::DiscreteElement<comp_d_dim_t>(right))
+                - ddc::coordinate(ddc::DiscreteElement<comp_d_dim_t>(left)));
+        return (tensor.get(right, ddc::DiscreteElement<comp_index_t>(0))
+                - tensor.get(left, ddc::DiscreteElement<comp_index_t>(0)))
+               / dx;
+    } else {
+        auto const right = aligned_elem;
+        auto const left = right - ddc::DiscreteVector<comp_d_dim_t>(1);
+        double const dx = static_cast<double>(
+                ddc::coordinate(ddc::DiscreteElement<comp_d_dim_t>(right))
+                - ddc::coordinate(ddc::DiscreteElement<comp_d_dim_t>(left)));
+        return (tensor.get(right, ddc::DiscreteElement<comp_index_t>(0))
+                - tensor.get(left, ddc::DiscreteElement<comp_index_t>(0)))
+               / dx;
+    }
+}
+
+template <
+        misc::Specialization<tensor::Tensor> OutTensorType,
+        class SupportTag,
+        class FirstComponent,
+        class SecondComponent,
+        class ExecSpace>
+OutTensorType exterior_derivative_of_tensor_form_2d(
+        ExecSpace const& exec_space,
+        OutTensorType out_tensor,
+        TensorForm<SupportTag, FirstComponent, SecondComponent> tensor_form)
+{
+    using out_index_t = ddc::type_seq_element_t<
+            0,
+            ddc::to_type_seq_t<typename OutTensorType::indices_domain_t>>;
+    using first_out_d_dim_t = mesher::detail::
+            discrete_dimension_for_t<typename FirstComponent::tag, typename OutTensorType::non_indices_domain_t>;
+    using second_out_d_dim_t = mesher::detail::
+            discrete_dimension_for_t<typename SecondComponent::tag, typename OutTensorType::non_indices_domain_t>;
+    static_assert(out_index_t::rank() == 0);
+
+    auto const first_tensor = tensor_form.template component<typename FirstComponent::tag>();
+    auto const second_tensor = tensor_form.template component<typename SecondComponent::tag>();
+
+    ddc::parallel_for_each(
+            "similie_compute_tensor_form_exterior_derivative_2d",
+            exec_space,
+            out_tensor.non_indices_domain(),
+            KOKKOS_LAMBDA(typename OutTensorType::non_indices_domain_t::discrete_element_type elem) {
+                if (is_on_primal_boundary<first_out_d_dim_t>(out_tensor.non_indices_domain(), elem)
+                    || is_on_primal_boundary<second_out_d_dim_t>(out_tensor.non_indices_domain(), elem)) {
+                    return;
+                }
+                out_tensor.mem(elem, ddc::DiscreteElement<out_index_t>(0))
+                        += centered_form_derivative<typename FirstComponent::tag, OutTensorType>(
+                                   elem,
+                                   second_tensor)
+                           - centered_form_derivative<typename SecondComponent::tag, OutTensorType>(
+                                   elem,
+                                   first_tensor);
+            });
+    return out_tensor;
+}
+
 } // namespace detail
 
 template <
@@ -279,6 +400,101 @@ codifferential_tensor_t<TagToRemoveFromCochain, CochainTag, TensorType> codiffer
     return codifferential_tensor;
 }
 
+template <
+        tensor::TensorIndex MetricIndex,
+        class TagToRemoveFromCochain,
+        tensor::TensorIndex CochainTag,
+        misc::Specialization<tensor::Tensor> OutTensorType,
+        misc::Specialization<tensor::Tensor> TensorType,
+        misc::Specialization<tensor::Tensor> MetricType,
+        class Dualizer,
+        class ExecSpace>
+OutTensorType codifferential(
+        ExecSpace const& exec_space,
+        OutTensorType out_tensor,
+        TensorType tensor,
+        MetricType inv_metric,
+        Dualizer const& dualizer)
+{
+    using out_index_t = ddc::type_seq_element_t<
+            0,
+            ddc::to_type_seq_t<typename OutTensorType::indices_domain_t>>;
+    using metric_component_index_t = ddc::type_seq_element_t<
+            0,
+            ddc::to_type_seq_t<typename MetricType::indices_domain_t>>;
+    using out_d_dim_t = mesher::detail::
+            discrete_dimension_for_t<TagToRemoveFromCochain, typename OutTensorType::non_indices_domain_t>;
+    using flux_d_dim_t = mesher::detail::
+            discrete_dimension_for_t<TagToRemoveFromCochain, typename TensorType::non_indices_domain_t>;
+    using out_elem_t = typename OutTensorType::non_indices_domain_t::discrete_element_type;
+    using flux_elem_t = typename TensorType::non_indices_domain_t::discrete_element_type;
+    using flux_index_t = ddc::type_seq_element_t<
+            0,
+            ddc::to_type_seq_t<typename TensorType::indices_domain_t>>;
+    static_assert(out_index_t::rank() == 0);
+    static_assert(!std::is_void_v<out_d_dim_t>);
+    static_assert(!std::is_void_v<flux_d_dim_t>);
+
+    ddc::parallel_for_each(
+            "similie_compute_centered_dualized_codifferential",
+            exec_space,
+            out_tensor.non_indices_domain(),
+            KOKKOS_LAMBDA(out_elem_t out_elem) {
+                ddc::DiscreteDomain<out_d_dim_t> dim_dom(out_tensor.non_indices_domain());
+                ddc::DiscreteElement<out_d_dim_t> dim_elem(out_elem);
+                if (dim_elem.uid() == dim_dom.front().uid() || dim_elem.uid() == dim_dom.back().uid()) {
+                    return;
+                }
+
+                flux_elem_t const right_face = dualizer(out_elem);
+                flux_elem_t const left_face = right_face - ddc::DiscreteVector<flux_d_dim_t>(1);
+                double const dx = static_cast<double>(
+                        ddc::coordinate(ddc::DiscreteElement<flux_d_dim_t>(right_face))
+                        - ddc::coordinate(ddc::DiscreteElement<flux_d_dim_t>(left_face)));
+                double const deriv = (tensor.get(right_face, ddc::DiscreteElement<flux_index_t>(0))
+                                      - tensor.get(left_face, ddc::DiscreteElement<flux_index_t>(0)))
+                                     / dx;
+                double const metric_factor = detail::codifferential_metric_factor<
+                        metric_component_index_t,
+                        TagToRemoveFromCochain>(inv_metric, out_elem);
+                out_tensor.mem(out_elem, ddc::DiscreteElement<out_index_t>(0))
+                        += metric_factor * deriv;
+            });
+    return out_tensor;
+}
+
+template <
+        tensor::TensorIndex MetricIndex,
+        misc::Specialization<tensor::Tensor> OutTensorType,
+        class SupportTag,
+        class FirstComponent,
+        class SecondComponent,
+        misc::Specialization<tensor::Tensor> MetricType,
+        class ExecSpace>
+OutTensorType codifferential(
+        ExecSpace const& exec_space,
+        OutTensorType out_tensor,
+        TensorForm<SupportTag, FirstComponent, SecondComponent> tensor_form,
+        MetricType inv_metric)
+{
+    ddc::Chunk dual_first_alloc(
+            typename SecondComponent::tensor_type::discrete_domain_type(
+                    tensor_form.template component<typename SecondComponent::tag>().domain()),
+            ddc::KokkosAllocator<double, typename ExecSpace::memory_space>());
+    sil::tensor::Tensor dual_first(dual_first_alloc);
+    ddc::Chunk dual_second_alloc(
+            typename FirstComponent::tensor_type::discrete_domain_type(
+                    tensor_form.template component<typename FirstComponent::tag>().domain()),
+            ddc::KokkosAllocator<double, typename ExecSpace::memory_space>());
+    sil::tensor::Tensor dual_second(dual_second_alloc);
+    auto dual_form = make_tensor_form<hodge_dual_support_t<SupportTag>>(
+            component<typename FirstComponent::tag>(dual_first),
+            component<typename SecondComponent::tag>(dual_second));
+
+    sil::exterior::hodge_star(exec_space, dual_form, tensor_form, inv_metric);
+    return detail::exterior_derivative_of_tensor_form_2d(exec_space, out_tensor, dual_form);
+}
+
 } // namespace exterior
 
 } // namespace sil
diff --git a/src/similie/exterior/form.hpp b/src/similie/exterior/form.hpp
index 99b2b38..c190c1e 100644
--- a/src/similie/exterior/form.hpp
+++ b/src/similie/exterior/form.hpp
@@ -3,8 +3,14 @@
 
 #pragma once
 
+#include <tuple>
+#include <type_traits>
+
 #include <ddc/ddc.hpp>
 
+#include <similie/misc/specialization.hpp>
+#include <similie/tensor/tensor_impl.hpp>
+
 #include "cochain.hpp"
 #include "cosimplex.hpp"
 
@@ -12,6 +18,32 @@ namespace sil {
 
 namespace exterior {
 
+struct PrimalSupport
+{
+};
+
+struct DualSupport
+{
+};
+
+template <class SupportTag>
+struct HodgeDualSupport;
+
+template <>
+struct HodgeDualSupport<PrimalSupport>
+{
+    using type = DualSupport;
+};
+
+template <>
+struct HodgeDualSupport<DualSupport>
+{
+    using type = PrimalSupport;
+};
+
+template <class SupportTag>
+using hodge_dual_support_t = typename HodgeDualSupport<SupportTag>::type;
+
 namespace detail {
 
 template <class T, class ElementType, class LayoutStridedPolicy>
@@ -43,6 +75,86 @@ FormWrapper(Head, ElementType, Allocator) -> FormWrapper<Head, ElementType, Layo
 template <class T, class ElementType = double, class LayoutStridedPolicy = Kokkos::LayoutRight>
 using Form = typename detail::FormWrapper<T, ElementType, LayoutStridedPolicy>::type;
 
+template <class Tag, misc::Specialization<tensor::Tensor> TensorType>
+struct FormComponent
+{
+    using tag = Tag;
+    using tensor_type = TensorType;
+
+    TensorType tensor;
+
+    KOKKOS_FUNCTION explicit constexpr FormComponent(TensorType tensor_) noexcept : tensor(tensor_) {}
+};
+
+template <class Tag, misc::Specialization<tensor::Tensor> TensorType>
+KOKKOS_FUNCTION constexpr FormComponent<Tag, TensorType> component(TensorType tensor) noexcept
+{
+    return FormComponent<Tag, TensorType>(tensor);
+}
+
+namespace detail {
+
+template <class Tag, class... Components>
+struct FormComponentByTag;
+
+template <class Tag, class HeadComponent, class... TailComponents>
+struct FormComponentByTag<Tag, HeadComponent, TailComponents...>
+{
+    using type = std::conditional_t<
+            std::is_same_v<Tag, typename HeadComponent::tag>,
+            HeadComponent,
+            typename FormComponentByTag<Tag, TailComponents...>::type>;
+};
+
+template <class Tag>
+struct FormComponentByTag<Tag>
+{
+    using type = void;
+};
+
+template <class... Components>
+constexpr bool has_unique_component_tags_v
+        = ((... && !std::is_same_v<typename Components::tag, typename Components::tag>) || true);
+
+template <class HeadComponent, class... TailComponents>
+constexpr bool has_unique_component_tags_v<HeadComponent, TailComponents...>
+        = ((... && !std::is_same_v<typename HeadComponent::tag, typename TailComponents::tag>))
+          && has_unique_component_tags_v<TailComponents...>;
+
+} // namespace detail
+
+template <class SupportTag, class... Components>
+class TensorForm
+{
+    static_assert(sizeof...(Components) > 0);
+    static_assert(detail::has_unique_component_tags_v<Components...>);
+
+    std::tuple<Components...> m_components;
+
+public:
+    using support_tag = SupportTag;
+    using components_tuple_type = std::tuple<Components...>;
+
+    constexpr explicit TensorForm(Components... components) noexcept
+        : m_components(std::move(components)...)
+    {
+    }
+
+    template <class Tag>
+    KOKKOS_FUNCTION constexpr auto component() const noexcept
+    {
+        using component_t = typename detail::FormComponentByTag<Tag, Components...>::type;
+        static_assert(!std::is_void_v<component_t>);
+        return std::get<component_t>(m_components).tensor;
+    }
+};
+
+template <class SupportTag = PrimalSupport, class... Components>
+constexpr TensorForm<SupportTag, Components...> make_tensor_form(Components... components) noexcept
+{
+    return TensorForm<SupportTag, Components...>(components...);
+}
+
 } // namespace exterior
 
 } // namespace sil
diff --git a/src/similie/exterior/hodge_star.hpp b/src/similie/exterior/hodge_star.hpp
index c3693e1..e60ad3a 100644
--- a/src/similie/exterior/hodge_star.hpp
+++ b/src/similie/exterior/hodge_star.hpp
@@ -7,6 +7,7 @@
 
 #include <similie/misc/factorial.hpp>
 #include <similie/misc/macros.hpp>
+#include <similie/mesher/dualizer.hpp>
 #include <similie/misc/specialization.hpp>
 #include <similie/misc/type_seq_conversion.hpp>
 #include <similie/tensor/antisymmetric_tensor.hpp>
@@ -17,10 +18,157 @@
 #include <similie/tensor/metric.hpp>
 #include <similie/tensor/prime.hpp>
 
+#include "form.hpp"
+
 namespace sil {
 
 namespace exterior {
 
+namespace detail {
+
+template <class Axis1, class Axis2, misc::Specialization<tensor::Tensor> MetricType, class Elem>
+KOKKOS_FUNCTION double inverse_metric_component(MetricType metric, Elem elem)
+{
+    using metric_component_index_t = ddc::type_seq_element_t<
+            0,
+            ddc::to_type_seq_t<typename MetricType::indices_domain_t>>;
+    if constexpr (misc::Specialization<metric_component_index_t, tensor::TensorIdentityIndex>) {
+        return std::is_same_v<Axis1, Axis2> ? 1. : 0.;
+    } else if constexpr (
+            misc::Specialization<metric_component_index_t, tensor::TensorDiagonalIndex>) {
+        if constexpr (std::is_same_v<Axis1, Axis2>) {
+            return metric(
+                    elem,
+                    metric.accessor().template access_element<Axis1, Axis2>());
+        } else {
+            return 0.;
+        }
+    } else {
+        return metric(
+                elem,
+                metric.accessor().template access_element<Axis1, Axis2>());
+    }
+}
+
+template <tensor::TensorIndex MetricIndex, misc::Specialization<tensor::Tensor> MetricType, class Elem>
+KOKKOS_FUNCTION double primal_volume_factor(MetricType metric, Elem elem)
+{
+    using metric_index_1_t = tensor::metric_index_1<tensor::upper_t<MetricIndex>>;
+    using axis1_t = ddc::type_seq_element_t<0, typename metric_index_1_t::type_seq_dimensions>;
+    using axis2_t = ddc::type_seq_element_t<1, typename metric_index_1_t::type_seq_dimensions>;
+    double const gxx = inverse_metric_component<axis1_t, axis1_t>(metric, elem);
+    double const gxy = inverse_metric_component<axis1_t, axis2_t>(metric, elem);
+    double const gyx = inverse_metric_component<axis2_t, axis1_t>(metric, elem);
+    double const gyy = inverse_metric_component<axis2_t, axis2_t>(metric, elem);
+    double const det_inv_metric = gxx * gyy - gxy * gyx;
+    return Kokkos::sqrt(Kokkos::abs(1. / det_inv_metric));
+}
+
+template <class Dom>
+struct IsFullyDualDomain;
+
+template <class... DDims>
+struct IsFullyDualDomain<ddc::DiscreteDomain<DDims...>>
+{
+    static constexpr bool value
+            = (true
+               && ... && !std::is_same_v<DDims, sil::mesher::detail::primal_discrete_dimension_t<DDims>>);
+};
+
+template <class Dom>
+constexpr bool is_fully_dual_domain_v = IsFullyDualDomain<Dom>::value;
+
+template <class OutTensor, class InTensor, class OutIndex, class InIndex>
+struct TensorFormHodgeStarFirstComponentFunctor
+{
+    OutTensor out_tensor;
+    InTensor in_tensor;
+
+    KOKKOS_FUNCTION void operator()(typename InTensor::non_indices_domain_t::discrete_element_type elem) const
+    {
+        out_tensor.mem(elem, ddc::DiscreteElement<OutIndex>(0))
+                = -in_tensor.get(elem, ddc::DiscreteElement<InIndex>(0));
+    }
+};
+
+template <class OutTensor, class InTensor, class OutIndex, class InIndex>
+struct TensorFormHodgeStarSecondComponentFunctor
+{
+    OutTensor out_tensor;
+    InTensor in_tensor;
+
+    KOKKOS_FUNCTION void operator()(typename InTensor::non_indices_domain_t::discrete_element_type elem) const
+    {
+        out_tensor.mem(elem, ddc::DiscreteElement<OutIndex>(0))
+                = in_tensor.get(elem, ddc::DiscreteElement<InIndex>(0));
+    }
+};
+
+template <
+        class FirstTag,
+        class SecondTag,
+        class SupportTag,
+        class... OutComponents,
+        class FirstComponent,
+        class SecondComponent,
+        class ExecSpace>
+void fill_tensor_form_hodge_star_2d(
+        ExecSpace const& exec_space,
+        TensorForm<hodge_dual_support_t<SupportTag>, OutComponents...> out_form,
+        TensorForm<SupportTag, FirstComponent, SecondComponent> in_form)
+{
+    using first_in_tensor_t = typename FirstComponent::tensor_type;
+    using second_in_tensor_t = typename SecondComponent::tensor_type;
+    using out_first_component_t = typename FormComponentByTag<FirstTag, OutComponents...>::type;
+    using out_second_component_t = typename FormComponentByTag<SecondTag, OutComponents...>::type;
+    using first_in_index_t = ddc::type_seq_element_t<
+            0,
+            ddc::to_type_seq_t<typename first_in_tensor_t::indices_domain_t>>;
+    using second_in_index_t = ddc::type_seq_element_t<
+            0,
+            ddc::to_type_seq_t<typename second_in_tensor_t::indices_domain_t>>;
+    using first_out_index_t = ddc::type_seq_element_t<
+            0,
+            ddc::to_type_seq_t<typename second_in_tensor_t::indices_domain_t>>;
+    using second_out_index_t = ddc::type_seq_element_t<
+            0,
+            ddc::to_type_seq_t<typename out_second_component_t::tensor_type::indices_domain_t>>;
+    static_assert(!std::is_void_v<out_first_component_t>);
+    static_assert(!std::is_void_v<out_second_component_t>);
+    static_assert(std::is_same_v<typename out_first_component_t::tensor_type, second_in_tensor_t>);
+    static_assert(std::is_same_v<typename out_second_component_t::tensor_type, first_in_tensor_t>);
+    static_assert(first_in_index_t::rank() == 0);
+    static_assert(second_in_index_t::rank() == 0);
+    static_assert(first_out_index_t::rank() == 0);
+    static_assert(second_out_index_t::rank() == 0);
+
+    auto const in_first = in_form.template component<FirstTag>();
+    auto const in_second = in_form.template component<SecondTag>();
+    auto const out_first = out_form.template component<FirstTag>();
+    auto const out_second = out_form.template component<SecondTag>();
+
+    ddc::parallel_for_each(
+            "similie_compute_tensor_form_hodge_star_first_component",
+            exec_space,
+            out_first.non_indices_domain(),
+            TensorFormHodgeStarFirstComponentFunctor<
+                    decltype(out_first),
+                    decltype(in_second),
+                    first_out_index_t,
+                    second_in_index_t> {out_first, in_second});
+    ddc::parallel_for_each(
+            "similie_compute_tensor_form_hodge_star_second_component",
+            exec_space,
+            out_second.non_indices_domain(),
+            TensorFormHodgeStarSecondComponentFunctor<
+                    decltype(out_second),
+                    decltype(in_first),
+                    second_out_index_t,
+                    first_in_index_t> {out_second, in_first});
+}
+
+} // namespace detail
+
 template <
         misc::Specialization<ddc::detail::TypeSeq> Indices1,
         misc::Specialization<ddc::detail::TypeSeq> Indices2>
@@ -161,6 +309,74 @@ HodgeStarType fill_hodge_star(
     return fill_hodge_star<Indices1, Indices2>(exec_space, hodge_star, metric_det, metric_prod);
 }
 
+template <
+        class SupportTag,
+        class... OutComponents,
+        class FirstComponent,
+        class SecondComponent,
+        class MetricType,
+        class ExecSpace>
+auto hodge_star(
+        ExecSpace const& exec_space,
+        TensorForm<hodge_dual_support_t<SupportTag>, OutComponents...> out_form,
+        TensorForm<SupportTag, FirstComponent, SecondComponent> in_form,
+        MetricType const&) -> TensorForm<hodge_dual_support_t<SupportTag>, OutComponents...>
+{
+    detail::fill_tensor_form_hodge_star_2d<
+            typename FirstComponent::tag,
+            typename SecondComponent::tag>(exec_space, out_form, in_form);
+    return out_form;
+}
+
+template <
+        tensor::TensorIndex MetricIndex,
+        misc::Specialization<tensor::Tensor> OutTensorType,
+        misc::Specialization<tensor::Tensor> InTensorType,
+        misc::Specialization<tensor::Tensor> MetricType,
+        class ExecSpace>
+    requires(
+            ddc::type_seq_element_t<0, ddc::to_type_seq_t<typename OutTensorType::indices_domain_t>>::rank()
+                    == 0
+            && ddc::type_seq_element_t<
+                               0,
+                               ddc::to_type_seq_t<typename InTensorType::indices_domain_t>>::rank()
+                       == 0
+            && OutTensorType::non_indices_domain_t::rank() == 2
+            && InTensorType::non_indices_domain_t::rank() == 2)
+OutTensorType hodge_star(
+        ExecSpace const& exec_space,
+        OutTensorType out_tensor,
+        InTensorType in_tensor,
+        MetricType metric)
+{
+    using out_index_t = ddc::type_seq_element_t<
+            0,
+            ddc::to_type_seq_t<typename OutTensorType::indices_domain_t>>;
+    using in_index_t = ddc::type_seq_element_t<
+            0,
+            ddc::to_type_seq_t<typename InTensorType::indices_domain_t>>;
+    using out_elem_t = typename OutTensorType::non_indices_domain_t::discrete_element_type;
+    using in_elem_t = typename InTensorType::non_indices_domain_t::discrete_element_type;
+    using metric_elem_t = typename MetricType::non_indices_domain_t::discrete_element_type;
+    constexpr bool in_is_fully_dual = detail::is_fully_dual_domain_v<typename InTensorType::non_indices_domain_t>;
+
+    ddc::parallel_for_each(
+            "similie_compute_scalar_tensor_hodge_star_2d",
+            exec_space,
+            out_tensor.non_indices_domain(),
+            KOKKOS_LAMBDA(out_elem_t out_elem) {
+                in_elem_t const in_elem
+                        = sil::mesher::detail::dualizer_remap_element<in_elem_t>(out_elem);
+                metric_elem_t const metric_elem
+                        = sil::mesher::detail::dualizer_remap_element<metric_elem_t>(out_elem);
+                double const volume_factor = detail::primal_volume_factor<MetricIndex>(metric, metric_elem);
+                double const scale = in_is_fully_dual ? 1. / volume_factor : volume_factor;
+                out_tensor.mem(out_elem, ddc::DiscreteElement<out_index_t>(0))
+                        = scale * in_tensor.get(in_elem, ddc::DiscreteElement<in_index_t>(0));
+            });
+    return out_tensor;
+}
+
 } // namespace exterior
 
 } // namespace sil
diff --git a/src/similie/mesher/dualizer.hpp b/src/similie/mesher/dualizer.hpp
new file mode 100644
index 0000000..2d6c001
--- /dev/null
+++ b/src/similie/mesher/dualizer.hpp
@@ -0,0 +1,273 @@
+// SPDX-FileCopyrightText: 2026 Baptiste Legouix
+// SPDX-License-Identifier: MIT
+
+#pragma once
+
+#include <type_traits>
+
+#include <ddc/ddc.hpp>
+
+namespace sil {
+
+namespace mesher {
+
+struct HalfShifted
+{
+};
+
+template <class DDim, class Policy = HalfShifted>
+struct Dual;
+
+template <class CDim>
+struct HalfShiftDualizer;
+
+template <class PrimalDDim>
+struct Dual<PrimalDDim, HalfShifted>
+    : ddc::UniformPointSampling<typename PrimalDDim::continuous_dimension_type>
+{
+    using primal_discrete_dimension_type = PrimalDDim;
+    using dual_policy = HalfShifted;
+};
+
+namespace detail {
+
+template <class CDim, class TypeSeq>
+struct DiscreteDimensionForTypeSeq;
+
+template <class CDim, class DDim, class... Tail>
+struct DiscreteDimensionForTypeSeq<CDim, ddc::detail::TypeSeq<DDim, Tail...>>
+{
+    using type = std::conditional_t<
+            std::is_same_v<typename DDim::continuous_dimension_type, CDim>,
+            DDim,
+            typename DiscreteDimensionForTypeSeq<CDim, ddc::detail::TypeSeq<Tail...>>::type>;
+};
+
+template <class CDim>
+struct DiscreteDimensionForTypeSeq<CDim, ddc::detail::TypeSeq<>>
+{
+    using type = void;
+};
+
+template <class CDim, class T>
+using discrete_dimension_for_t
+        = typename DiscreteDimensionForTypeSeq<CDim, ddc::to_type_seq_t<T>>::type;
+
+template <class Dualizer, class DDim>
+struct DualDiscreteDimension;
+
+template <class CDim, class DDim>
+struct DualDiscreteDimension<HalfShiftDualizer<CDim>, DDim>
+{
+    using type = std::conditional_t<
+            std::is_same_v<typename DDim::continuous_dimension_type, CDim>,
+            Dual<DDim, HalfShifted>,
+            DDim>;
+};
+
+template <class DDim, class = void>
+struct PrimalDiscreteDimension
+{
+    using type = DDim;
+};
+
+template <class DDim>
+struct PrimalDiscreteDimension<DDim, std::void_t<typename DDim::primal_discrete_dimension_type>>
+{
+    using type = typename DDim::primal_discrete_dimension_type;
+};
+
+template <class DDim>
+using primal_discrete_dimension_t = typename PrimalDiscreteDimension<DDim>::type;
+
+template <class Dualizer, class T>
+struct Dualized;
+
+template <class Dualizer, class... DDims>
+struct Dualized<Dualizer, ddc::DiscreteElement<DDims...>>
+{
+    using type = ddc::DiscreteElement<typename DualDiscreteDimension<Dualizer, DDims>::type...>;
+};
+
+template <class Dualizer, class... DDims>
+struct Dualized<Dualizer, ddc::DiscreteVector<DDims...>>
+{
+    using type = ddc::DiscreteVector<typename DualDiscreteDimension<Dualizer, DDims>::type...>;
+};
+
+template <class Dualizer, class... DDims>
+struct Dualized<Dualizer, ddc::DiscreteDomain<DDims...>>
+{
+    using type = ddc::DiscreteDomain<typename DualDiscreteDimension<Dualizer, DDims>::type...>;
+};
+
+template <class Dualizer, class T>
+using dualized_t = typename Dualized<Dualizer, T>::type;
+
+template <class T>
+struct Primalized;
+
+template <class... DDims>
+struct Primalized<ddc::DiscreteElement<DDims...>>
+{
+    using type = ddc::DiscreteElement<primal_discrete_dimension_t<DDims>...>;
+};
+
+template <class... DDims>
+struct Primalized<ddc::DiscreteVector<DDims...>>
+{
+    using type = ddc::DiscreteVector<primal_discrete_dimension_t<DDims>...>;
+};
+
+template <class... DDims>
+struct Primalized<ddc::DiscreteDomain<DDims...>>
+{
+    using type = ddc::DiscreteDomain<primal_discrete_dimension_t<DDims>...>;
+};
+
+template <class T>
+using primalized_t = typename Primalized<T>::type;
+
+template <class TargetDDim, class SourceElem>
+KOKKOS_FUNCTION ddc::DiscreteElement<TargetDDim> dualizer_remap_dimension(SourceElem const& source_elem)
+{
+    using source_d_dim_t
+            = discrete_dimension_for_t<typename TargetDDim::continuous_dimension_type, SourceElem>;
+    static_assert(!std::is_void_v<source_d_dim_t>);
+    return ddc::DiscreteElement<TargetDDim>(ddc::uid<source_d_dim_t>(source_elem));
+}
+
+template <class TargetElem, class SourceElem>
+struct DualizerRemapElement;
+
+template <class... TargetDDims, class SourceElem>
+struct DualizerRemapElement<ddc::DiscreteElement<TargetDDims...>, SourceElem>
+{
+    KOKKOS_FUNCTION static ddc::DiscreteElement<TargetDDims...> apply(SourceElem const& source_elem)
+    {
+        return ddc::DiscreteElement<TargetDDims...>(
+                dualizer_remap_dimension<TargetDDims>(source_elem)...);
+    }
+};
+
+template <class TargetElem, class SourceElem>
+KOKKOS_FUNCTION TargetElem dualizer_remap_element(SourceElem const& source_elem)
+{
+    return DualizerRemapElement<TargetElem, SourceElem>::apply(source_elem);
+}
+
+template <class Dualizer, class DDim>
+auto dualize_subdomain(ddc::DiscreteDomain<DDim> const& subdomain)
+{
+    using dual_d_dim_t = typename DualDiscreteDimension<Dualizer, DDim>::type;
+    if constexpr (std::is_same_v<dual_d_dim_t, DDim>) {
+        return subdomain;
+    } else {
+        static_assert(ddc::is_uniform_point_sampling_v<DDim>);
+        auto const nb_points = subdomain.template extent<DDim>().value();
+        assert(nb_points > 1);
+        auto const lower = ddc::coordinate(subdomain.front())
+                           + ddc::Coordinate<typename DDim::continuous_dimension_type>(
+                                   ddc::step<DDim>() / 2.);
+        auto const upper = ddc::coordinate(subdomain.back())
+                           - ddc::Coordinate<typename DDim::continuous_dimension_type>(
+                                   ddc::step<DDim>() / 2.);
+        return ddc::init_discrete_space<dual_d_dim_t>(dual_d_dim_t::template init<dual_d_dim_t>(
+                lower,
+                upper,
+                ddc::DiscreteVector<dual_d_dim_t>(nb_points - 1)));
+    }
+}
+
+template <class Dualizer, class T>
+struct DualizerApplication;
+
+template <class Dualizer, class... DDims>
+struct DualizerApplication<Dualizer, ddc::DiscreteElement<DDims...>>
+{
+    using input_type = ddc::DiscreteElement<DDims...>;
+    using output_type = dualized_t<Dualizer, input_type>;
+
+    KOKKOS_FUNCTION static output_type apply(input_type const& elem)
+    {
+        return dualizer_remap_element<output_type>(elem);
+    }
+};
+
+template <class Dualizer, class... DDims>
+struct DualizerApplication<Dualizer, ddc::DiscreteVector<DDims...>>
+{
+    using input_type = ddc::DiscreteVector<DDims...>;
+    using output_type = dualized_t<Dualizer, input_type>;
+
+    KOKKOS_FUNCTION static output_type apply(input_type const& vect)
+    {
+        return output_type(
+                ddc::select<DDims>(vect).value()...);
+    }
+};
+
+template <class Dualizer, class... DDims>
+struct DualizerApplication<Dualizer, ddc::DiscreteDomain<DDims...>>
+{
+    using input_type = ddc::DiscreteDomain<DDims...>;
+    using output_type = dualized_t<Dualizer, input_type>;
+
+    static output_type apply(input_type const& dom)
+    {
+        return output_type(dualize_subdomain<Dualizer>(ddc::DiscreteDomain<DDims>(dom))...);
+    }
+};
+
+template <class T>
+struct PrimalizerApplication;
+
+template <class... DDims>
+struct PrimalizerApplication<ddc::DiscreteElement<DDims...>>
+{
+    using input_type = ddc::DiscreteElement<DDims...>;
+    using output_type = primalized_t<input_type>;
+
+    KOKKOS_FUNCTION static output_type apply(input_type const& elem)
+    {
+        return dualizer_remap_element<output_type>(elem);
+    }
+};
+
+template <class... DDims>
+struct PrimalizerApplication<ddc::DiscreteVector<DDims...>>
+{
+    using input_type = ddc::DiscreteVector<DDims...>;
+    using output_type = primalized_t<input_type>;
+
+    KOKKOS_FUNCTION static output_type apply(input_type const& vect)
+    {
+        return output_type(
+                ddc::select<DDims>(vect).value()...);
+    }
+};
+
+} // namespace detail
+
+template <class Dualizer, class DDim>
+using dual_discrete_dimension_t = typename detail::DualDiscreteDimension<Dualizer, DDim>::type;
+
+template <class CDim>
+struct HalfShiftDualizer
+{
+    template <class T>
+    constexpr detail::dualized_t<HalfShiftDualizer<CDim>, T> operator()(T const& obj) const
+    {
+        return detail::DualizerApplication<HalfShiftDualizer<CDim>, T>::apply(obj);
+    }
+
+    template <class T>
+    constexpr detail::primalized_t<T> primal(T const& obj) const
+    {
+        return detail::PrimalizerApplication<T>::apply(obj);
+    }
+};
+
+} // namespace mesher
+
+} // namespace sil
diff --git a/src/similie/similie.hpp b/src/similie/similie.hpp
index 8af73fb..69c3178 100644
--- a/src/similie/similie.hpp
+++ b/src/similie/similie.hpp
@@ -10,6 +10,7 @@ namespace sil {
 
 #include "csr/csr.hpp"
 #include "exterior/exterior.hpp"
+#include "mesher/dualizer.hpp"
 #include "mesher/mesher.hpp"
 #include "tensor/tensor.hpp"
 #include "young_tableau/young_tableau.hpp"
diff --git a/tests/exterior/exterior.cpp b/tests/exterior/exterior.cpp
index c19afd2..87a2578 100644
--- a/tests/exterior/exterior.cpp
+++ b/tests/exterior/exterior.cpp
@@ -41,6 +41,14 @@ struct DDimZ : ddc::UniformPointSampling<Z>
 {
 };
 
+struct DDimX2 : ddc::UniformPointSampling<X>
+{
+};
+
+struct DDimY2 : ddc::UniformPointSampling<Y>
+{
+};
+
 TEST(Chain, Optimization)
 {
     sil::exterior::Chain chain = sil::exterior::
@@ -514,6 +522,144 @@ TEST(Form, Alias)
                     0.);
 }
 
+TEST(Form, TensorFormDeriv)
+{
+    using DummyIndex = sil::tensor::Covariant<sil::tensor::ScalarIndex>;
+    using XDualizer = sil::mesher::HalfShiftDualizer<X>;
+    using YDualizer = sil::mesher::HalfShiftDualizer<Y>;
+    using DDimXDual = sil::mesher::dual_discrete_dimension_t<XDualizer, DDimX>;
+    using DDimYDual = sil::mesher::dual_discrete_dimension_t<YDualizer, DDimY>;
+
+    auto const x_dom = ddc::init_discrete_space<DDimX>(DDimX::init<DDimX>(
+            ddc::Coordinate<X>(0.),
+            ddc::Coordinate<X>(1.),
+            ddc::DiscreteVector<DDimX>(6)));
+    auto const y_dom = ddc::init_discrete_space<DDimY>(DDimY::init<DDimY>(
+            ddc::Coordinate<Y>(0.),
+            ddc::Coordinate<Y>(1.),
+            ddc::DiscreteVector<DDimY>(6)));
+    ddc::DiscreteDomain<DDimX, DDimY> const mesh(x_dom, y_dom);
+    XDualizer const x_dualizer;
+    YDualizer const y_dualizer;
+    ddc::DiscreteDomain<DDimXDual, DDimY> const x_face_dom = x_dualizer(mesh);
+    ddc::DiscreteDomain<DDimX, DDimYDual> const y_face_dom = y_dualizer(mesh);
+
+    [[maybe_unused]] sil::tensor::TensorAccessor<DummyIndex> scalar_accessor;
+    ddc::Chunk scalar_alloc(
+            ddc::DiscreteDomain<DDimX, DDimY, DummyIndex>(mesh, scalar_accessor.domain()),
+            ddc::HostAllocator<double>());
+    sil::tensor::Tensor scalar(scalar_alloc);
+    ddc::host_for_each(mesh, [&](ddc::DiscreteElement<DDimX, DDimY> elem) {
+        double const x = ddc::coordinate(ddc::DiscreteElement<DDimX>(elem));
+        double const y = ddc::coordinate(ddc::DiscreteElement<DDimY>(elem));
+        scalar(elem, ddc::DiscreteElement<DummyIndex>(0)) = x * x + y;
+    });
+
+    ddc::Chunk grad_x_alloc(
+            ddc::DiscreteDomain<DDimXDual, DDimY, DummyIndex>(
+                    x_face_dom,
+                    scalar_accessor.domain()),
+            ddc::HostAllocator<double>());
+    sil::tensor::Tensor grad_x(grad_x_alloc);
+    ddc::Chunk grad_y_alloc(
+            ddc::DiscreteDomain<DDimX, DDimYDual, DummyIndex>(
+                    y_face_dom,
+                    scalar_accessor.domain()),
+            ddc::HostAllocator<double>());
+    sil::tensor::Tensor grad_y(grad_y_alloc);
+    auto grad_form = sil::exterior::make_tensor_form(
+            sil::exterior::component<X>(grad_x),
+            sil::exterior::component<Y>(grad_y));
+
+    sil::exterior::deriv(Kokkos::DefaultHostExecutionSpace(), grad_form, scalar);
+
+    ddc::DiscreteElement<DDimXDual, DDimY> const x_face_elem
+            = x_face_dom.front() + ddc::DiscreteVector<DDimXDual, DDimY>(2, 2);
+    ddc::DiscreteElement<DDimX, DDimYDual> const y_face_elem
+            = y_face_dom.front() + ddc::DiscreteVector<DDimX, DDimYDual>(2, 2);
+    EXPECT_NEAR(grad_x(x_face_elem, ddc::DiscreteElement<DummyIndex>(0)), 1.0, 1e-12);
+    EXPECT_NEAR(grad_y(y_face_elem, ddc::DiscreteElement<DummyIndex>(0)), 1.0, 1e-12);
+}
+
+TEST(Form, TensorFormCodifferential)
+{
+    using MetricIndex = sil::tensor::TensorIdentityIndex<
+            sil::tensor::Covariant<sil::tensor::MetricIndex1<X, Y>>,
+            sil::tensor::Covariant<sil::tensor::MetricIndex2<X, Y>>>;
+    using InverseMetricIndex = sil::tensor::upper_t<MetricIndex>;
+    using DummyIndex = sil::tensor::Covariant<sil::tensor::ScalarIndex>;
+    using XDualizer = sil::mesher::HalfShiftDualizer<X>;
+    using YDualizer = sil::mesher::HalfShiftDualizer<Y>;
+    using DDimXDual = sil::mesher::dual_discrete_dimension_t<XDualizer, DDimX2>;
+    using DDimYDual = sil::mesher::dual_discrete_dimension_t<YDualizer, DDimY2>;
+
+    auto const x_dom = ddc::init_discrete_space<DDimX2>(DDimX2::init<DDimX2>(
+            ddc::Coordinate<X>(0.),
+            ddc::Coordinate<X>(1.),
+            ddc::DiscreteVector<DDimX2>(6)));
+    auto const y_dom = ddc::init_discrete_space<DDimY2>(DDimY2::init<DDimY2>(
+            ddc::Coordinate<Y>(0.),
+            ddc::Coordinate<Y>(1.),
+            ddc::DiscreteVector<DDimY2>(6)));
+    ddc::DiscreteDomain<DDimX2, DDimY2> const mesh(x_dom, y_dom);
+    XDualizer const x_dualizer;
+    YDualizer const y_dualizer;
+    ddc::DiscreteDomain<DDimXDual, DDimY2> const x_face_dom = x_dualizer(mesh);
+    ddc::DiscreteDomain<DDimX2, DDimYDual> const y_face_dom = y_dualizer(mesh);
+
+    [[maybe_unused]] sil::tensor::TensorAccessor<DummyIndex> scalar_accessor;
+    ddc::Chunk grad_x_alloc(
+            ddc::DiscreteDomain<DDimXDual, DDimY2, DummyIndex>(
+                    x_face_dom,
+                    scalar_accessor.domain()),
+            ddc::HostAllocator<double>());
+    sil::tensor::Tensor grad_x(grad_x_alloc);
+    ddc::Chunk grad_y_alloc(
+            ddc::DiscreteDomain<DDimX2, DDimYDual, DummyIndex>(
+                    y_face_dom,
+                    scalar_accessor.domain()),
+            ddc::HostAllocator<double>());
+    sil::tensor::Tensor grad_y(grad_y_alloc);
+
+    ddc::host_for_each(grad_x.domain(), [&](ddc::DiscreteElement<DDimXDual, DDimY2, DummyIndex> elem) {
+        double const x = ddc::coordinate(ddc::DiscreteElement<DDimXDual>(elem));
+        grad_x(elem) = 2. * x;
+    });
+    ddc::host_for_each(grad_y.domain(), [&](ddc::DiscreteElement<DDimX2, DDimYDual, DummyIndex> elem) {
+        grad_y(elem) = 1.;
+    });
+
+    auto form = sil::exterior::make_tensor_form(
+            sil::exterior::component<X>(grad_x),
+            sil::exterior::component<Y>(grad_y));
+
+    [[maybe_unused]] sil::tensor::TensorAccessor<InverseMetricIndex> inv_metric_accessor;
+    ddc::Chunk inv_metric_alloc(
+            ddc::DiscreteDomain<DDimX2, DDimY2, InverseMetricIndex>(
+                    mesh,
+                    inv_metric_accessor.domain()),
+            ddc::HostAllocator<double>());
+    sil::tensor::Tensor inv_metric(inv_metric_alloc);
+
+    ddc::Chunk div_alloc(
+            ddc::DiscreteDomain<DDimX2, DDimY2, DummyIndex>(mesh, scalar_accessor.domain()),
+            ddc::HostAllocator<double>());
+    sil::tensor::Tensor div(div_alloc);
+    ddc::host_for_each(div.domain(), [&](ddc::DiscreteElement<DDimX2, DDimY2, DummyIndex> elem) {
+        div(elem) = 0.;
+    });
+
+    sil::exterior::codifferential<MetricIndex>(
+            Kokkos::DefaultHostExecutionSpace(),
+            div,
+            form,
+            inv_metric);
+
+    ddc::DiscreteElement<DDimX2, DDimY2> const center
+            = mesh.front() + ddc::DiscreteVector<DDimX2, DDimY2>(2, 2);
+    EXPECT_NEAR(div(center, ddc::DiscreteElement<DummyIndex>(0)), 2.0, 1e-12);
+}
+
 TEST(Cochain, Test)
 {
     sil::exterior::Chain
diff --git a/tests/exterior/hodge_star.cpp b/tests/exterior/hodge_star.cpp
index fc8b76f..da345bd 100644
--- a/tests/exterior/hodge_star.cpp
+++ b/tests/exterior/hodge_star.cpp
@@ -7,6 +7,7 @@
 
 #include <gtest/gtest.h>
 #include <similie/exterior/hodge_star.hpp>
+#include <similie/mesher/dualizer.hpp>
 #include <similie/tensor/levi_civita_tensor.hpp>
 #include <similie/tensor/metric.hpp>
 #include <similie/tensor/symmetric_tensor.hpp>
@@ -35,6 +36,22 @@ struct DDimZ : ddc::UniformPointSampling<Z>
 {
 };
 
+struct X2
+{
+};
+
+struct Y2
+{
+};
+
+struct DDimX2 : ddc::UniformPointSampling<X2>
+{
+};
+
+struct DDimY2 : ddc::UniformPointSampling<Y2>
+{
+};
+
 struct Mu : sil::tensor::TensorNaturalIndex<X, Y, Z>
 {
 };
@@ -63,6 +80,8 @@ using HodgeStarDomain = sil::exterior::
 using HodgeStarDomain2 = sil::exterior::
         hodge_star_domain_t<ddc::detail::TypeSeq<RhoUp>, ddc::detail::TypeSeq<MuLow, NuLow>>;
 
+using DummyIndex = sil::tensor::Covariant<sil::tensor::ScalarIndex>;
+
 TEST(HodgeStar, Test)
 {
     ddc::DiscreteDomain<DDimX, DDimY>
@@ -137,3 +156,164 @@ TEST(HodgeStar, Test)
         EXPECT_DOUBLE_EQ(form(elem, form.accessor().access_element<Y, Z>()), 3.);
     });
 }
+
+TEST(HodgeStar, TensorForm1In2D)
+{
+    using XDualizer = sil::mesher::HalfShiftDualizer<X>;
+    using YDualizer = sil::mesher::HalfShiftDualizer<Y>;
+    using DDimXDual = sil::mesher::dual_discrete_dimension_t<XDualizer, DDimX>;
+    using DDimYDual = sil::mesher::dual_discrete_dimension_t<YDualizer, DDimY>;
+
+    auto const x_dom = ddc::init_discrete_space<DDimX>(DDimX::init<DDimX>(
+            ddc::Coordinate<X>(0.),
+            ddc::Coordinate<X>(1.),
+            ddc::DiscreteVector<DDimX>(4)));
+    auto const y_dom = ddc::init_discrete_space<DDimY>(DDimY::init<DDimY>(
+            ddc::Coordinate<Y>(0.),
+            ddc::Coordinate<Y>(1.),
+            ddc::DiscreteVector<DDimY>(4)));
+    ddc::DiscreteDomain<DDimX, DDimY> const mesh(x_dom, y_dom);
+    XDualizer const x_dualizer;
+    YDualizer const y_dualizer;
+    ddc::DiscreteDomain<DDimXDual, DDimY> const x_face_dom = x_dualizer(mesh);
+    ddc::DiscreteDomain<DDimX, DDimYDual> const y_face_dom = y_dualizer(mesh);
+
+    [[maybe_unused]] sil::tensor::TensorAccessor<DummyIndex> scalar_accessor;
+
+    ddc::Chunk primal_x_alloc(
+            ddc::DiscreteDomain<DDimXDual, DDimY, DummyIndex>(x_face_dom, scalar_accessor.domain()),
+            ddc::HostAllocator<double>());
+    sil::tensor::Tensor primal_x(primal_x_alloc);
+    ddc::Chunk primal_y_alloc(
+            ddc::DiscreteDomain<DDimX, DDimYDual, DummyIndex>(y_face_dom, scalar_accessor.domain()),
+            ddc::HostAllocator<double>());
+    sil::tensor::Tensor primal_y(primal_y_alloc);
+
+    ddc::host_for_each(primal_x.domain(), [&](ddc::DiscreteElement<DDimXDual, DDimY, DummyIndex> elem) {
+        primal_x(elem) = 2.;
+    });
+    ddc::host_for_each(primal_y.domain(), [&](ddc::DiscreteElement<DDimX, DDimYDual, DummyIndex> elem) {
+        primal_y(elem) = 3.;
+    });
+
+    auto primal_form = sil::exterior::make_tensor_form(
+            sil::exterior::component<X>(primal_x),
+            sil::exterior::component<Y>(primal_y));
+
+    ddc::Chunk dual_x_alloc(
+            ddc::DiscreteDomain<DDimX, DDimYDual, DummyIndex>(y_face_dom, scalar_accessor.domain()),
+            ddc::HostAllocator<double>());
+    sil::tensor::Tensor dual_x(dual_x_alloc);
+    ddc::Chunk dual_y_alloc(
+            ddc::DiscreteDomain<DDimXDual, DDimY, DummyIndex>(x_face_dom, scalar_accessor.domain()),
+            ddc::HostAllocator<double>());
+    sil::tensor::Tensor dual_y(dual_y_alloc);
+    auto dual_form = sil::exterior::make_tensor_form<sil::exterior::DualSupport>(
+            sil::exterior::component<X>(dual_x),
+            sil::exterior::component<Y>(dual_y));
+
+    sil::exterior::hodge_star(Kokkos::DefaultHostExecutionSpace(), dual_form, primal_form, 0);
+
+    EXPECT_DOUBLE_EQ(
+            dual_x(y_face_dom.front(), ddc::DiscreteElement<DummyIndex>(0)),
+            -3.);
+    EXPECT_DOUBLE_EQ(
+            dual_y(x_face_dom.front(), ddc::DiscreteElement<DummyIndex>(0)),
+            2.);
+
+    ddc::Chunk primal_back_x_alloc(
+            ddc::DiscreteDomain<DDimXDual, DDimY, DummyIndex>(x_face_dom, scalar_accessor.domain()),
+            ddc::HostAllocator<double>());
+    sil::tensor::Tensor primal_back_x(primal_back_x_alloc);
+    ddc::Chunk primal_back_y_alloc(
+            ddc::DiscreteDomain<DDimX, DDimYDual, DummyIndex>(y_face_dom, scalar_accessor.domain()),
+            ddc::HostAllocator<double>());
+    sil::tensor::Tensor primal_back_y(primal_back_y_alloc);
+    auto primal_back = sil::exterior::make_tensor_form(
+            sil::exterior::component<X>(primal_back_x),
+            sil::exterior::component<Y>(primal_back_y));
+
+    sil::exterior::hodge_star(Kokkos::DefaultHostExecutionSpace(), primal_back, dual_form, 0);
+
+    EXPECT_DOUBLE_EQ(
+            primal_back_x(x_face_dom.front(), ddc::DiscreteElement<DummyIndex>(0)),
+            -2.);
+    EXPECT_DOUBLE_EQ(
+            primal_back_y(y_face_dom.front(), ddc::DiscreteElement<DummyIndex>(0)),
+            -3.);
+}
+
+TEST(HodgeStar, ScalarTensor0And2In2D)
+{
+    using ScalarIndex = sil::tensor::Covariant<sil::tensor::ScalarIndex>;
+    using Metric2DIndex = sil::tensor::TensorDiagonalIndex<
+            sil::tensor::Contravariant<sil::tensor::MetricIndex1<X2, Y2>>,
+            sil::tensor::Contravariant<sil::tensor::MetricIndex2<X2, Y2>>>;
+    using XDualizer = sil::mesher::HalfShiftDualizer<X2>;
+    using YDualizer = sil::mesher::HalfShiftDualizer<Y2>;
+    using DDimXDual = sil::mesher::dual_discrete_dimension_t<XDualizer, DDimX2>;
+    using DDimYDual = sil::mesher::dual_discrete_dimension_t<YDualizer, DDimY2>;
+
+    auto const x_dom = ddc::init_discrete_space<DDimX2>(DDimX2::init<DDimX2>(
+            ddc::Coordinate<X2>(0.),
+            ddc::Coordinate<X2>(1.),
+            ddc::DiscreteVector<DDimX2>(4)));
+    auto const y_dom = ddc::init_discrete_space<DDimY2>(DDimY2::init<DDimY2>(
+            ddc::Coordinate<Y2>(0.),
+            ddc::Coordinate<Y2>(1.),
+            ddc::DiscreteVector<DDimY2>(4)));
+    ddc::DiscreteDomain<DDimX2, DDimY2> const mesh(x_dom, y_dom);
+    XDualizer const x_dualizer;
+    YDualizer const y_dualizer;
+    auto const dual_mesh = y_dualizer(x_dualizer(mesh));
+
+    [[maybe_unused]] sil::tensor::TensorAccessor<ScalarIndex> scalar_accessor;
+    ddc::Chunk scalar_alloc(
+            ddc::DiscreteDomain<DDimX2, DDimY2, ScalarIndex>(mesh, scalar_accessor.domain()),
+            ddc::HostAllocator<double>());
+    sil::tensor::Tensor scalar(scalar_alloc);
+    ddc::host_for_each(scalar.domain(), [&](ddc::DiscreteElement<DDimX2, DDimY2, ScalarIndex> elem) {
+        scalar(elem) = 2.;
+    });
+
+    [[maybe_unused]] sil::tensor::TensorAccessor<Metric2DIndex> metric_accessor;
+    ddc::Chunk metric_alloc(
+            ddc::DiscreteDomain<DDimX2, DDimY2, Metric2DIndex>(mesh, metric_accessor.domain()),
+            ddc::HostAllocator<double>());
+    sil::tensor::Tensor metric(metric_alloc);
+    ddc::host_for_each(mesh, [&](ddc::DiscreteElement<DDimX2, DDimY2> elem) {
+        metric(elem, metric.accessor().access_element<X2, X2>()) = 4.;
+        metric(elem, metric.accessor().access_element<Y2, Y2>()) = 9.;
+    });
+
+    ddc::Chunk dual_scalar_alloc(
+            ddc::DiscreteDomain<DDimXDual, DDimYDual, ScalarIndex>(dual_mesh, scalar_accessor.domain()),
+            ddc::HostAllocator<double>());
+    sil::tensor::Tensor dual_scalar(dual_scalar_alloc);
+
+    sil::exterior::hodge_star<Metric2DIndex>(
+            Kokkos::DefaultHostExecutionSpace(),
+            dual_scalar,
+            scalar,
+            metric);
+
+    EXPECT_DOUBLE_EQ(dual_scalar(dual_scalar.domain().front()), 1. / 3.);
+
+    ddc::DiscreteDomain<DDimX2, DDimY2> const primal_submesh(
+            mesh.front(),
+            ddc::DiscreteVector<DDimX2, DDimY2>(
+                    dual_mesh.template extent<DDimXDual>().value(),
+                    dual_mesh.template extent<DDimYDual>().value()));
+    ddc::Chunk scalar_back_alloc(
+            ddc::DiscreteDomain<DDimX2, DDimY2, ScalarIndex>(primal_submesh, scalar_accessor.domain()),
+            ddc::HostAllocator<double>());
+    sil::tensor::Tensor scalar_back(scalar_back_alloc);
+    ddc::parallel_fill(scalar_back, 0.);
+    sil::exterior::hodge_star<Metric2DIndex>(
+            Kokkos::DefaultHostExecutionSpace(),
+            scalar_back,
+            dual_scalar,
+            metric);
+
+    EXPECT_DOUBLE_EQ(scalar_back(scalar_back.domain().front()), 2.);
+}