diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index 10258da6e888cbb9a8c5a9834583c641abbe6b03..b414b30e5b71d4c60b731673a7dd55ccb2f9b6e9 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -10,6 +10,7 @@
 #include "Process.h"
 
 #include "NumLib/DOF/ComputeSparsityPattern.h"
+#include "NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h"
 #include "DirichletBc.h"
 #include "NeumannBc.h"
 #include "NeumannBcAssembler.h"
@@ -52,6 +53,9 @@ void Process::initialize()
     DBUG("Compute sparsity pattern");
     computeSparsityPattern();
 
+    DBUG("Initialize the extrapolator");
+    initializeExtrapolator();
+
     initializeConcreteProcess(*_local_to_global_index_map, _mesh,
                               _integration_order);
 
@@ -168,6 +172,48 @@ void Process::constructDofTable()
         std::move(all_mesh_subsets), NumLib::ComponentOrder::BY_LOCATION));
 }
 
+void Process::initializeExtrapolator()
+{
+    NumLib::LocalToGlobalIndexMap const* dof_table_single_component;
+    bool manage_storage;
+
+    if (_local_to_global_index_map->getNumberOfComponents() == 1)
+    {
+        // For single-variable-single-component processes reuse the existing DOF
+        // table.
+        dof_table_single_component = _local_to_global_index_map.get();
+        manage_storage = false;
+    }
+    else
+    {
+        // Otherwise construct a new DOF table.
+        std::vector<std::unique_ptr<MeshLib::MeshSubsets>>
+            all_mesh_subsets_single_component;
+        all_mesh_subsets_single_component.emplace_back(
+            new MeshLib::MeshSubsets(this->_mesh_subset_all_nodes.get()));
+
+        dof_table_single_component = new NumLib::LocalToGlobalIndexMap(
+            std::move(all_mesh_subsets_single_component),
+            // by location order is needed for output
+            NumLib::ComponentOrder::BY_LOCATION);
+        manage_storage = true;
+    }
+
+    MathLib::MatrixSpecifications mat_specs(
+        dof_table_single_component->dofSizeWithoutGhosts(),
+        dof_table_single_component->dofSizeWithoutGhosts(),
+        &dof_table_single_component->getGhostIndices(),
+        nullptr);
+
+    std::unique_ptr<NumLib::Extrapolator> extrapolator(
+        new NumLib::LocalLinearLeastSquaresExtrapolator(
+            mat_specs, *dof_table_single_component));
+
+    // TODO Later on the DOF table can change during the simulation!
+    _extrapolator_data = ExtrapolatorData(
+        std::move(extrapolator), dof_table_single_component, manage_storage);
+}
+
 void Process::setInitialConditions(ProcessVariable const& variable,
                                    int const variable_id,
                                    int const component_id,
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index 72d4009202b733b08518e5116292c273fcb5d3f9..bb9da6ed7015b4ea5414902a91cf99f24328995e 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -14,6 +14,7 @@
 #include "NumLib/ODESolver/ODESystem.h"
 #include "NumLib/ODESolver/TimeDiscretization.h"
 
+#include "ExtrapolatorData.h"
 #include "Parameter.h"
 #include "ProcessOutput.h"
 #include "SecondaryVariable.h"
@@ -85,6 +86,12 @@ public:
         return *_time_discretization;
     }
 
+protected:
+    NumLib::Extrapolator& getExtrapolator()
+    {
+        return *_extrapolator_data.extrapolator;
+    }
+
 private:
     /// Process specific initialization called by initialize().
     virtual void initializeConcreteProcess(
@@ -104,6 +111,8 @@ private:
 
     void constructDofTable();
 
+    void initializeExtrapolator();
+
     /// Sets the initial condition values in the solution vector x for a given
     /// process variable and component.
     void setInitialConditions(ProcessVariable const& variable,
@@ -142,6 +151,8 @@ private:
 
     /// Variables used by this process.
     std::vector<std::reference_wrapper<ProcessVariable>> _process_variables;
+
+    ExtrapolatorData _extrapolator_data;
 };
 
 /// Find process variables in \c variables whose names match the settings under