diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index d1fa4bf7ce816e054d36a216cf66fc1084a4c6c7..2e72192ea23d93067549ab19ef78316264e89704 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -639,17 +639,16 @@ SchemeModule::SchemeModule()
                               }));
 
   this->_addBuiltinFunction("hybridheat",
-                            std::make_shared<BuiltinFunctionEmbedder<std::tuple<
-                              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>&,
-                                                                        const std::vector<std::shared_ptr<
-                                                                          const IBoundaryConditionDescriptor>>&,
-                                                                        const std::vector<std::shared_ptr<
-                                                                          const IBoundaryConditionDescriptor>>&,
-                                                                        const std::shared_ptr<const IDiscreteFunction>&,
-                                                                        const double&, const double&, const double&)>>(
+                            std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
+                              const IDiscreteFunction>(const std::shared_ptr<const IDiscreteFunction>&,
+                                                       const std::shared_ptr<const IDiscreteFunction>&,
+                                                       const std::shared_ptr<const IDiscreteFunction>&,
+                                                       const std::vector<
+                                                         std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+                                                       const std::vector<
+                                                         std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+                                                       const std::shared_ptr<const IDiscreteFunction>&, const double&,
+                                                       const double&, const double&)>>(
 
                               [](const std::shared_ptr<const IDiscreteFunction>& kappa,
                                  const std::shared_ptr<const IDiscreteFunction>& f,
@@ -660,8 +659,7 @@ SchemeModule::SchemeModule()
                                    bc_prev_descriptor_list,
                                  const std::shared_ptr<const IDiscreteFunction>& P, const double& dt,
                                  const double& cn_coeff,
-                                 const double& lambda) -> const std::tuple<std::shared_ptr<const IDiscreteFunction>,
-                                                                           std::shared_ptr<const IDiscreteFunction>> {
+                                 const double& lambda) -> const std::shared_ptr<const IDiscreteFunction> {
                                 return ScalarHybridSchemeHandler{kappa,
                                                                  f,
                                                                  f_prev,
diff --git a/src/scheme/ScalarHybridScheme.cpp b/src/scheme/ScalarHybridScheme.cpp
index 187d452b65ef2d56d95c7c1d72e25c78df12d1d4..f0c4414e97208ed5fb253ca15151ea25358a7e23 100644
--- a/src/scheme/ScalarHybridScheme.cpp
+++ b/src/scheme/ScalarHybridScheme.cpp
@@ -6,8 +6,7 @@
 class ScalarHybridSchemeHandler::IScalarHybridScheme
 {
  public:
-  virtual std::tuple<std::shared_ptr<const IDiscreteFunction>, std::shared_ptr<const IDiscreteFunction>> getSolution()
-    const = 0;
+  virtual std::shared_ptr<const IDiscreteFunction> getSolution() const = 0;
 
   // virtual std::shared_ptr<const IDiscreteFunction> getSolution() const = 0;
 
@@ -24,7 +23,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
   using MeshDataType     = MeshData<Dimension>;
 
   std::shared_ptr<const DiscreteFunctionP0<Dimension, double>> m_cell_temperature;
-  std::shared_ptr<const DiscreteFunctionP0<Dimension, double>> m_node_temperature;
   //  std::shared_ptr<const DiscreteFunctionP0<Dimension, double>> m_solution;
 
   class DirichletBoundaryCondition
@@ -153,10 +151,10 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
   };
 
  public:
-  std::tuple<std::shared_ptr<const IDiscreteFunction>, std::shared_ptr<const IDiscreteFunction>>
+  std::shared_ptr<const IDiscreteFunction>
   getSolution() const
   {
-    return {m_cell_temperature, m_node_temperature};
+    return {m_cell_temperature};
   }
   // std::shared_ptr<const IDiscreteFunction>
   // getSolution() const final
@@ -337,11 +335,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
 
     MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
 
-    std::shared_ptr dual_mesh = DualMeshManager::instance().getMedianDualMesh(*mesh);
-
-    std::shared_ptr mapper =
-      DualConnectivityManager::instance().getPrimalToMedianDualConnectivityDataMapper(mesh->connectivity());
-
     const NodeValue<const TinyVector<Dimension>>& xr = mesh->xr();
 
     const FaceValue<const TinyVector<Dimension>>& xl = mesh_data.xl();
@@ -1496,23 +1489,13 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
       return compute_primal_node_temperature;
     }();
 
-    const CellValue<const double> dual_node_temperature = [=] {
-      CellValue<double> my_dual_temperature{dual_mesh->connectivity()};
-      mapper->toDualCell(primal_node_temperature, my_dual_temperature);
-      return my_dual_temperature;
-    }();
-
     // m_solution     = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh);
     m_cell_temperature = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh);
-    m_node_temperature = std::make_shared<DiscreteFunctionP0<Dimension, double>>(dual_mesh);
     // auto& solution = *m_solution;
     auto& cell_temperature = *m_cell_temperature;
-    auto& node_temperature = *m_node_temperature;
+
     parallel_for(
       mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { cell_temperature[cell_id] = T[cell_id]; });
-    parallel_for(
-      dual_mesh->numberOfCells(),
-      PUGS_LAMBDA(CellId cell_id) { node_temperature[cell_id] = dual_node_temperature[cell_id]; });
   }
 };
 
@@ -1521,7 +1504,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
 // {
 //   return m_scheme->getSolution();
 // }
-std::tuple<std::shared_ptr<const IDiscreteFunction>, std::shared_ptr<const IDiscreteFunction>>
+std::shared_ptr<const IDiscreteFunction>
 ScalarHybridSchemeHandler::solution() const
 {
   return m_scheme->getSolution();
@@ -1567,7 +1550,19 @@ ScalarHybridSchemeHandler::ScalarHybridSchemeHandler(
     break;
   }
   case 3: {
-    throw NormalError("The scheme is not implemented in dimension 3");
+    using MeshType                   = Mesh<Connectivity<3>>;
+    using DiscreteTensorFunctionType = DiscreteFunctionP0<3, TinyMatrix<3>>;
+    using DiscreteScalarFunctionType = DiscreteFunctionP0<3, double>;
+
+    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+
+    m_scheme =
+      std::make_unique<ScalarHybridScheme<3>>(mesh, std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k),
+                                              std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
+                                              std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f_prev),
+                                              bc_descriptor_list, bc_prev_descriptor_list,
+                                              std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(P), dt,
+                                              cn_coeff, lambda);
     break;
   }
   default: {
diff --git a/src/scheme/ScalarHybridScheme.hpp b/src/scheme/ScalarHybridScheme.hpp
index b110887d7bcce04d4456f8185bb6201ed983a8ae..9de2c37cd382a3692aab3bec81438468a8edffb3 100644
--- a/src/scheme/ScalarHybridScheme.hpp
+++ b/src/scheme/ScalarHybridScheme.hpp
@@ -40,7 +40,7 @@ class ScalarHybridSchemeHandler
   std::unique_ptr<IScalarHybridScheme> m_scheme;
 
   // std::shared_ptr<const IDiscreteFunction> solution() const;
-  std::tuple<std::shared_ptr<const IDiscreteFunction>, std::shared_ptr<const IDiscreteFunction>> solution() const;
+  std::shared_ptr<const IDiscreteFunction> solution() const;
 
   ScalarHybridSchemeHandler(
     const std::shared_ptr<const IDiscreteFunction>& cell_k,