diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index ecba28f40d2199858b79033aa5c4c9ee3fd6b394..f47c26354426ce62402425cedd3bb432da7c3b86 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -493,7 +493,7 @@ SchemeModule::SchemeModule()
 
   this->_addBuiltinFunction("moleculardiffusion",
                             std::make_shared<BuiltinFunctionEmbedder<std::tuple<
-                              std::shared_ptr<const IDiscreteFunction>,
+                              std::shared_ptr<const IDiscreteFunction>, std::shared_ptr<const IDiscreteFunction>,
                               std::shared_ptr<const IDiscreteFunction>>(const std::shared_ptr<const IDiscreteFunction>&,
                                                                         const std::shared_ptr<const IDiscreteFunction>&,
                                                                         const std::shared_ptr<const IDiscreteFunction>&,
@@ -510,6 +510,7 @@ SchemeModule::SchemeModule()
                                  const std::shared_ptr<const IDiscreteFunction> f,
                                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                    bc_descriptor_list) -> const std::tuple<std::shared_ptr<const IDiscreteFunction>,
+                                                                           std::shared_ptr<const IDiscreteFunction>,
                                                                            std::shared_ptr<const IDiscreteFunction>> {
                                 return VectorDiamondSchemeHandler{alpha, lambdab,           mub, lambda, mu,
                                                                   f,     bc_descriptor_list}
diff --git a/src/scheme/VectorDiamondScheme.cpp b/src/scheme/VectorDiamondScheme.cpp
index cbf198ecdeab61c3d4022cef73d3cd385cd86274..f50accc056d474c7d73287a7166842be20fcc896 100644
--- a/src/scheme/VectorDiamondScheme.cpp
+++ b/src/scheme/VectorDiamondScheme.cpp
@@ -158,8 +158,10 @@ class VectorDiamondSchemeHandler::IVectorDiamondScheme
  public:
   virtual std::shared_ptr<const IDiscreteFunction> getSolution() const     = 0;
   virtual std::shared_ptr<const IDiscreteFunction> getDualSolution() const = 0;
-  virtual std::tuple<std::shared_ptr<const IDiscreteFunction>, std::shared_ptr<const IDiscreteFunction>> apply()
-    const = 0;
+  virtual std::tuple<std::shared_ptr<const IDiscreteFunction>,
+                     std::shared_ptr<const IDiscreteFunction>,
+                     std::shared_ptr<const IDiscreteFunction>>
+  apply() const = 0;
 
   IVectorDiamondScheme()          = default;
   virtual ~IVectorDiamondScheme() = default;
@@ -175,7 +177,7 @@ class VectorDiamondSchemeHandler::VectorDiamondScheme : public VectorDiamondSche
 
   std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>> m_solution;
   std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>> m_dual_solution;
-  std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>> m_fluxes;
+  std::shared_ptr<const DiscreteFunctionP0<Dimension, double>> m_energy_delta;
 
   class DirichletBoundaryCondition
   {
@@ -267,10 +269,12 @@ class VectorDiamondSchemeHandler::VectorDiamondScheme : public VectorDiamondSche
     return m_dual_solution;
   }
 
-  std::tuple<std::shared_ptr<const IDiscreteFunction>, std::shared_ptr<const IDiscreteFunction>>
+  std::tuple<std::shared_ptr<const IDiscreteFunction>,
+             std::shared_ptr<const IDiscreteFunction>,
+             std::shared_ptr<const IDiscreteFunction>>
   apply() const final
   {
-    return {m_solution, m_dual_solution};
+    return {m_solution, m_dual_solution, m_energy_delta};
   }
 
   VectorDiamondScheme(const std::shared_ptr<const MeshType>& mesh,
@@ -945,7 +949,7 @@ class VectorDiamondSchemeHandler::VectorDiamondScheme : public VectorDiamondSche
     }
   }
   // compute the fluxes
-  std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>>
+  void   // std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>
   computeFluxes(const std::shared_ptr<const MeshType>& mesh,
                 const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& alpha,
                 const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& dual_lambdab,
@@ -1200,8 +1204,6 @@ class VectorDiamondSchemeHandler::VectorDiamondScheme : public VectorDiamondSche
         std::shared_ptr mapper =
           DiamondDualConnectivityManager::instance().getConnectivityToDiamondDualConnectivityDataMapper(
             mesh->connectivity());
-        //        m_fluxes = std::make_shared<DiscreteFunctionP0<Dimension, TinyVector<Dimension>>>(diamond_mesh);
-        //        auto& fluxes = *m_fluxes;
 
         CellValue<const double> dual_muj      = dual_mu->cellValues();
         CellValue<const double> dual_lambdaj  = dual_lambda->cellValues();
@@ -1476,15 +1478,34 @@ class VectorDiamondSchemeHandler::VectorDiamondScheme : public VectorDiamondSche
           }
           // exit(0);
         }
+        // Assemble
+        m_energy_delta     = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh);
+        auto& energy_delta = *m_energy_delta;
+        for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+          energy_delta[cell_id] = 0;
+        }
+        double sum_deltae = 0.;
+        // CellValue<double>& deltae = m_energy_delta->cellValues();
+        for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+          for (size_t j = 0; j < face_to_cell_matrix[face_id].size(); j++) {
+            CellId i_id = face_to_cell_matrix[face_id][j];
+            energy_delta[i_id] += flux(face_id, j) / primal_Vj[i_id];
+            sum_deltae += flux(face_id, j);
+          }
+          // exit(0);
+        }
+        std::cout << "sum deltat e " << sum_deltae << "\n";
       }
     } else {
       throw NotImplementedError("not done in 1d");
     }
-    return m_fluxes;
+    //    return m_energy_delta;
   }
 };
 
-std::tuple<std::shared_ptr<const IDiscreteFunction>, std::shared_ptr<const IDiscreteFunction>>
+std::tuple<std::shared_ptr<const IDiscreteFunction>,
+           std::shared_ptr<const IDiscreteFunction>,
+           std::shared_ptr<const IDiscreteFunction>>
 VectorDiamondSchemeHandler::apply() const
 {
   return m_scheme->apply();
diff --git a/src/scheme/VectorDiamondScheme.hpp b/src/scheme/VectorDiamondScheme.hpp
index a0f03a89d5e58cb7582dee9c7659263004d9259d..e2c2e0fef484529ea858c98d9715efb6811be3ce 100644
--- a/src/scheme/VectorDiamondScheme.hpp
+++ b/src/scheme/VectorDiamondScheme.hpp
@@ -39,7 +39,10 @@ class VectorDiamondSchemeHandler
 
   std::shared_ptr<const IDiscreteFunction> dual_solution() const;
 
-  std::tuple<std::shared_ptr<const IDiscreteFunction>, std::shared_ptr<const IDiscreteFunction>> apply() const;
+  std::tuple<std::shared_ptr<const IDiscreteFunction>,
+             std::shared_ptr<const IDiscreteFunction>,
+             std::shared_ptr<const IDiscreteFunction>>
+  apply() const;
 
   VectorDiamondSchemeHandler(
     const std::shared_ptr<const IDiscreteFunction>& alphab,