diff --git a/ProcessLib/HeatConduction/HeatConductionFEM.h b/ProcessLib/HeatConduction/HeatConductionFEM.h
index 738ee4a8c46d883fb7602b1534ac69c494a065e8..cd412cac74d00427cbe38a1d0dea425989ae485c 100644
--- a/ProcessLib/HeatConduction/HeatConductionFEM.h
+++ b/ProcessLib/HeatConduction/HeatConductionFEM.h
@@ -80,10 +80,10 @@ public:
         (void)local_matrix_size;
     }
 
-    void assemble(double const t, std::vector<double> const& local_x,
+    void assemble_HeatConduction(double const t, std::vector<double> const& local_x,
                   std::vector<double>& local_M_data,
                   std::vector<double>& local_K_data,
-                  std::vector<double>& /*local_b_data*/) override
+                  std::vector<double>& /*local_b_data*/)
     {
         auto const local_matrix_size = local_x.size();
         // This assertion is valid only if all nodal d.o.f. use the same shape
@@ -118,6 +118,25 @@ public:
         }
     }
 
+    void assemble(double const t, std::vector<double> const& local_x,
+                  std::vector<double>& local_M_data,
+                  std::vector<double>& local_K_data,
+                  std::vector<double>& local_b_data) override
+    {
+        assemble_HeatConduction(t, local_x, local_M_data, local_K_data,
+                                local_b_data);
+    }
+
+    void coupling_assemble(double const t, std::vector<double> const& local_x,
+                           std::vector<double>& local_M_data,
+                           std::vector<double>& local_K_data,
+                           std::vector<double>& local_b_data,
+                           LocalCouplingTerm const& coupled_term) override
+    {
+        assemble_HeatConduction(t, local_x, local_M_data, local_K_data,
+                                local_b_data);
+    }
+
     void computeSecondaryVariableConcrete(
                                     const double t,
                                     std::vector<double> const& local_x) override
diff --git a/ProcessLib/HeatConduction/HeatConductionProcess.cpp b/ProcessLib/HeatConduction/HeatConductionProcess.cpp
index 1f3782e6aa5ffa629dd229f5bbe66f1d1af3be18..012fb88bbc8960daa42839f5005df0d05851de8b 100644
--- a/ProcessLib/HeatConduction/HeatConductionProcess.cpp
+++ b/ProcessLib/HeatConduction/HeatConductionProcess.cpp
@@ -33,6 +33,22 @@ HeatConductionProcess::HeatConductionProcess(
 {
 }
 
+void HeatConductionProcess::preTimestepConcreteProcess(GlobalVector const& x,
+                                            const double /*t*/,
+                                            const double /*delta_t*/)
+{
+    if (!_x_previous_timestep)
+    {
+        _x_previous_timestep =
+            MathLib::MatrixVectorTraits<GlobalVector>::newInstance(x);
+    }
+    else
+    {
+        auto x0 = _x_previous_timestep.get();
+        MathLib::LinAlg::copy(x, *x0);
+    }
+}
+
 void HeatConductionProcess::initializeConcreteProcess(
     NumLib::LocalToGlobalIndexMap const& dof_table,
     MeshLib::Mesh const& mesh,
diff --git a/ProcessLib/HeatConduction/HeatConductionProcess.h b/ProcessLib/HeatConduction/HeatConductionProcess.h
index c6eaf6b103ae2c926d1b0f258ca772399e365cab..06f54405591188eba3e57320842261030b2906bc 100644
--- a/ProcessLib/HeatConduction/HeatConductionProcess.h
+++ b/ProcessLib/HeatConduction/HeatConductionProcess.h
@@ -41,6 +41,15 @@ public:
     void computeSecondaryVariableConcrete(double const t,
                                           GlobalVector const& x) override;
 
+    void preTimestepConcreteProcess(GlobalVector const& x, const double t,
+                                    const double delta_t) override;
+
+    // Get the solution of the previous time step.
+    virtual GlobalVector* getPreviousTimeStepSolution() const override
+    {
+        return _x_previous_timestep.get();
+    }
+
 private:
     void initializeConcreteProcess(
         NumLib::LocalToGlobalIndexMap const& dof_table,
@@ -63,6 +72,9 @@ private:
 
     std::vector<std::unique_ptr<HeatConductionLocalAssemblerInterface>>
         _local_assemblers;
+
+    /// Solution of the previous time step
+    std::unique_ptr<GlobalVector> _x_previous_timestep = nullptr;
 };
 
 }  // namespace HeatConduction