diff --git a/ProcessLib/TES/TESLocalAssembler-impl.h b/ProcessLib/TES/TESLocalAssembler-impl.h
index 331f67b63523159a21b2fcf18215b9d90e4061e7..391b29307367829b5a3552c200e0c983a3a04696 100644
--- a/ProcessLib/TES/TESLocalAssembler-impl.h
+++ b/ProcessLib/TES/TESLocalAssembler-impl.h
@@ -194,10 +194,47 @@ template <typename ShapeFunction_, typename IntegrationMethod_,
           unsigned GlobalDim>
 std::vector<double> const& TESLocalAssembler<
     ShapeFunction_, IntegrationMethod_,
-    GlobalDim>::getIntegrationPointValues(TESIntPtVariables const var,
-                                          std::vector<double>& cache) const
+    GlobalDim>::getIntPtSolidDensity(std::vector<double>& /*cache*/) const
 {
-    return _d.getIntegrationPointValues(var, cache);
+    return _d.getData().solid_density;
+}
+
+template <typename ShapeFunction_, typename IntegrationMethod_,
+          unsigned GlobalDim>
+std::vector<double> const& TESLocalAssembler<
+    ShapeFunction_, IntegrationMethod_,
+    GlobalDim>::getIntPtReactionRate(std::vector<double>& /*cache*/) const
+{
+    return _d.getData().reaction_rate;
+}
+
+template <typename ShapeFunction_, typename IntegrationMethod_,
+          unsigned GlobalDim>
+std::vector<double> const& TESLocalAssembler<
+    ShapeFunction_, IntegrationMethod_,
+    GlobalDim>::getIntPtDarcyVelocityX(std::vector<double>& /*cache*/) const
+{
+    return _d.getData().velocity[0];
+}
+
+template <typename ShapeFunction_, typename IntegrationMethod_,
+          unsigned GlobalDim>
+std::vector<double> const& TESLocalAssembler<
+    ShapeFunction_, IntegrationMethod_,
+    GlobalDim>::getIntPtDarcyVelocityY(std::vector<double>& /*cache*/) const
+{
+    assert(_d.getData().velocity.size() >= 2);
+    return _d.getData().velocity[1];
+}
+
+template <typename ShapeFunction_, typename IntegrationMethod_,
+          unsigned GlobalDim>
+std::vector<double> const& TESLocalAssembler<
+    ShapeFunction_, IntegrationMethod_,
+    GlobalDim>::getIntPtDarcyVelocityZ(std::vector<double>& /*cache*/) const
+{
+    assert(_d.getData().velocity.size() >= 3);
+    return _d.getData().velocity[2];
 }
 
 template <typename ShapeFunction_, typename IntegrationMethod_,
diff --git a/ProcessLib/TES/TESLocalAssembler.h b/ProcessLib/TES/TESLocalAssembler.h
index 37ce25c30459fb2f945a8eac62ec0086c0589fb6..5509c0cacf4a32a075975c1727e4d1a450e80d2f 100644
--- a/ProcessLib/TES/TESLocalAssembler.h
+++ b/ProcessLib/TES/TESLocalAssembler.h
@@ -17,7 +17,7 @@
 #include "TESAssemblyParams.h"
 #include "TESLocalAssemblerInner-fwd.h"
 
-#include "NumLib/Extrapolation/Extrapolator.h"
+#include "NumLib/Extrapolation/ExtrapolatableElement.h"
 
 namespace ProcessLib
 {
@@ -25,13 +25,28 @@ namespace TES
 {
 class TESLocalAssemblerInterface
     : public ProcessLib::LocalAssemblerInterface,
-      public NumLib::Extrapolatable<TESIntPtVariables>
+      public NumLib::ExtrapolatableElement
 {
 public:
     virtual ~TESLocalAssemblerInterface() = default;
 
     virtual bool checkBounds(std::vector<double> const& local_x,
                              std::vector<double> const& local_x_prev_ts) = 0;
+
+    virtual std::vector<double> const& getIntPtSolidDensity(
+        std::vector<double>& /*cache*/) const = 0;
+
+    virtual std::vector<double> const& getIntPtReactionRate(
+        std::vector<double>& /*cache*/) const = 0;
+
+    virtual std::vector<double> const& getIntPtDarcyVelocityX(
+        std::vector<double>& /*cache*/) const = 0;
+
+    virtual std::vector<double> const& getIntPtDarcyVelocityY(
+        std::vector<double>& /*cache*/) const = 0;
+
+    virtual std::vector<double> const& getIntPtDarcyVelocityZ(
+        std::vector<double>& /*cache*/) const = 0;
 };
 
 template <typename ShapeFunction_, typename IntegrationMethod_,
@@ -66,9 +81,20 @@ public:
     bool checkBounds(std::vector<double> const& local_x,
                      std::vector<double> const& local_x_prev_ts) override;
 
-    std::vector<double> const& getIntegrationPointValues(
-        TESIntPtVariables const var, std::vector<double>& cache) const override;
+    std::vector<double> const& getIntPtSolidDensity(
+        std::vector<double>& /*cache*/) const override;
+
+    std::vector<double> const& getIntPtReactionRate(
+        std::vector<double>& /*cache*/) const override;
+
+    std::vector<double> const& getIntPtDarcyVelocityX(
+        std::vector<double>& /*cache*/) const override;
+
+    std::vector<double> const& getIntPtDarcyVelocityY(
+        std::vector<double>& /*cache*/) const override;
 
+    std::vector<double> const& getIntPtDarcyVelocityZ(
+        std::vector<double>& /*cache*/) const override;
 private:
     std::vector<ShapeMatrices> _shape_matrices;
 
diff --git a/ProcessLib/TES/TESLocalAssemblerInner-impl.h b/ProcessLib/TES/TESLocalAssemblerInner-impl.h
index d587920df981851367fd1f4e838c8fb2412595d5..b0d63166cd60ad493793a48d369cc3de3a234d52 100644
--- a/ProcessLib/TES/TESLocalAssemblerInner-impl.h
+++ b/ProcessLib/TES/TESLocalAssemblerInner-impl.h
@@ -218,57 +218,6 @@ void TESLocalAssemblerInner<Traits>::preEachAssembleIntegrationPoint(
     _d.rho_GR = fluid_density(_d.p, _d.T, _d.vapour_mass_fraction);
 }
 
-template <typename Traits>
-std::vector<double> const&
-TESLocalAssemblerInner<Traits>::getIntegrationPointValues(
-    TESIntPtVariables const var, std::vector<double>& cache) const
-{
-    switch (var)
-    {
-        case TESIntPtVariables::REACTION_RATE:
-            return _d.reaction_rate;
-        case TESIntPtVariables::SOLID_DENSITY:
-            return _d.solid_density;
-        case TESIntPtVariables::VELOCITY_X:
-            return _d.velocity[0];
-        case TESIntPtVariables::VELOCITY_Y:
-            assert(_d.velocity.size() >= 2);
-            return _d.velocity[1];
-        case TESIntPtVariables::VELOCITY_Z:
-            assert(_d.velocity.size() >= 3);
-            return _d.velocity[2];
-
-        case TESIntPtVariables::LOADING:
-        {
-            auto& Cs = cache;
-            Cs.clear();
-            Cs.reserve(_d.solid_density.size());
-
-            for (auto rho_SR : _d.solid_density)
-            {
-                Cs.push_back(rho_SR / _d.ap.rho_SR_dry - 1.0);
-            }
-
-            return Cs;
-        }
-
-        // TODO that's an element value, ain't it?
-        case TESIntPtVariables::REACTION_DAMPING_FACTOR:
-        {
-            auto const num_integration_points = _d.solid_density.size();
-            auto& alphas = cache;
-            alphas.clear();
-            alphas.resize(num_integration_points,
-                          _d.reaction_adaptor->getReactionDampingFactor());
-
-            return alphas;
-        }
-    }
-
-    cache.clear();
-    return cache;
-}
-
 template <typename Traits>
 void TESLocalAssemblerInner<Traits>::assembleIntegrationPoint(
     unsigned integration_point,
diff --git a/ProcessLib/TES/TESLocalAssemblerInner.h b/ProcessLib/TES/TESLocalAssemblerInner.h
index 3e36cb4718059fb80bb4ffd989008ced34fcb05d..77d4638fc507c22d697a7cc89dee6f0bc5b39587 100644
--- a/ProcessLib/TES/TESLocalAssemblerInner.h
+++ b/ProcessLib/TES/TESLocalAssemblerInner.h
@@ -23,17 +23,6 @@ namespace ProcessLib
 {
 namespace TES
 {
-enum class TESIntPtVariables : unsigned
-{
-    SOLID_DENSITY,
-    REACTION_RATE,
-    VELOCITY_X,
-    VELOCITY_Y,
-    VELOCITY_Z,
-    LOADING,
-    REACTION_DAMPING_FACTOR
-};
-
 template <typename Traits>
 class TESLocalAssemblerInner
 {
@@ -56,9 +45,6 @@ public:
 
     void preEachAssemble();
 
-    std::vector<double> const& getIntegrationPointValues(
-        TESIntPtVariables const var, std::vector<double>& cache) const;
-
     // TODO better encapsulation
     AssemblyParams const& getAssemblyParameters() const { return _d.ap; }
     TESFEMReactionAdaptor const& getReactionAdaptor() const
diff --git a/ProcessLib/TES/TESProcess.cpp b/ProcessLib/TES/TESProcess.cpp
index b35eb3c189645aa58fc0ca75b8874d2bc1bbbd69..a73f5571f1adb4ce193d3c98a6cfcc39c73415b2 100644
--- a/ProcessLib/TES/TESProcess.cpp
+++ b/ProcessLib/TES/TESProcess.cpp
@@ -204,22 +204,27 @@ void TESProcess::initializeConcreteProcess(
                                                         std::move(fcts));
     };
     auto makeEx =
-        [&](TESIntPtVariables var) -> SecondaryVariableFunctions {
-        return ProcessLib::makeExtrapolator(var, *_extrapolator,
-                                            _local_assemblers);
+        [&](std::vector<double> const& (TESLocalAssemblerInterface::*method)(
+            std::vector<double>&)const) -> SecondaryVariableFunctions {
+        return ProcessLib::makeExtrapolator(*_extrapolator, _local_assemblers,
+                                            method);
     };
 
-    add2nd("solid_density", 1, makeEx(TESIntPtVariables::SOLID_DENSITY));
-    add2nd("reaction_rate", 1, makeEx(TESIntPtVariables::REACTION_RATE));
-    add2nd("velocity_x", 1, makeEx(TESIntPtVariables::VELOCITY_X));
+    // add2nd("solid_density", 1, makeEx(TESIntPtVariables::SOLID_DENSITY));
+    add2nd("reaction_rate", 1,
+           makeEx(&TESLocalAssemblerInterface::getIntPtReactionRate));
+    add2nd("velocity_x", 1,
+           makeEx(&TESLocalAssemblerInterface::getIntPtDarcyVelocityX));
     if (mesh.getDimension() >= 2)
-        add2nd("velocity_y", 1, makeEx(TESIntPtVariables::VELOCITY_Y));
+        add2nd("velocity_y", 1,
+               makeEx(&TESLocalAssemblerInterface::getIntPtDarcyVelocityY));
     if (mesh.getDimension() >= 3)
-        add2nd("velocity_z", 1, makeEx(TESIntPtVariables::VELOCITY_Z));
+        add2nd("velocity_z", 1,
+               makeEx(&TESLocalAssemblerInterface::getIntPtDarcyVelocityZ));
 
-    add2nd("loading", 1, makeEx(TESIntPtVariables::LOADING));
+    /*add2nd("loading", 1, makeEx(TESIntPtVariables::LOADING));
     add2nd("reaction_damping_factor", 1,
-           makeEx(TESIntPtVariables::REACTION_DAMPING_FACTOR));
+           makeEx(TESIntPtVariables::REACTION_DAMPING_FACTOR));*/
 
     namespace PH = std::placeholders;
     using Self = TESProcess;
diff --git a/ProcessLib/TES/TESProcess.h b/ProcessLib/TES/TESProcess.h
index 987f84be056ade216626e93881447609a1217989..3f670b68db058379135ca06411e0cebf042109d1 100644
--- a/ProcessLib/TES/TESProcess.h
+++ b/ProcessLib/TES/TESProcess.h
@@ -53,10 +53,9 @@ public:
     bool isLinear() const override { return false; }
 private:
     using ExtrapolatorInterface =
-        NumLib::Extrapolator<TESIntPtVariables, TESLocalAssemblerInterface>;
+        NumLib::Extrapolator;
     using ExtrapolatorImplementation =
-        NumLib::LocalLinearLeastSquaresExtrapolator<
-            TESIntPtVariables, TESLocalAssemblerInterface>;
+        NumLib::LocalLinearLeastSquaresExtrapolator;
 
     void initializeConcreteProcess(
         NumLib::LocalToGlobalIndexMap const& dof_table,