diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index d494bdec19bd28494d9f2c071b1cb6c4410c9763..e8fbfdcd4980b1533c42ceb2ad407fe203ebcb5b 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/Applications/ApplicationsLib/ProjectData.cpp
@@ -472,7 +472,7 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config,
         else if (type == "THERMO_MECHANICS")
         {
             //! \ogs_file_param{prj__processes__process__THERMO_MECHANICS__dimension}
-            switch (process_config.getConfigParameter<int>("dimension"))
+            switch (_mesh_vec[0]->getDimension())
             {
                 case 2:
                     process = ProcessLib::ThermoMechanics::
diff --git a/ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.cpp b/ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.cpp
index a35ab72aefc3f165ba576efe76ab51a61b786a32..e785b49dc3df2f49110c943e701804d65cccd015 100644
--- a/ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.cpp
+++ b/ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.cpp
@@ -23,12 +23,6 @@ namespace ProcessLib
 {
 namespace ThermoMechanics
 {
-template <int DisplacementDim>
-class ThermoMechanicsProcess;
-
-extern template class ThermoMechanicsProcess<2>;
-extern template class ThermoMechanicsProcess<3>;
-
 template <int DisplacementDim>
 std::unique_ptr<Process> createThermoMechanicsProcess(
     MeshLib::Mesh& mesh,
@@ -182,11 +176,10 @@ std::unique_ptr<Process> createThermoMechanicsProcess(
     ProcessLib::parseSecondaryVariables(config, secondary_variables,
                                         named_function_caller);
 
-    return std::unique_ptr<ThermoMechanicsProcess<DisplacementDim>>{
-        new ThermoMechanicsProcess<DisplacementDim>{
-            mesh, std::move(jacobian_assembler), parameters, integration_order,
-            std::move(process_variables), std::move(process_data),
-            std::move(secondary_variables), std::move(named_function_caller)}};
+    return std::make_unique<ThermoMechanicsProcess<DisplacementDim>>(
+        mesh, std::move(jacobian_assembler), parameters, integration_order,
+        std::move(process_variables), std::move(process_data),
+        std::move(secondary_variables), std::move(named_function_caller));
 }
 
 template std::unique_ptr<Process> createThermoMechanicsProcess<2>(
diff --git a/ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.h b/ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.h
index 5b2849eedd96e7b430cca2782e5b59092be2c4ab..161c0847ee931872bb37fe0059efda6fb0d14273 100644
--- a/ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.h
+++ b/ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.h
@@ -9,7 +9,24 @@
 
 #pragma once
 
-#include "ProcessLib/Process.h"
+#include <memory>
+#include <vector>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+namespace MeshLib
+{
+class Mesh;
+}
+namespace ProcessLib
+{
+class AbstractJacobianAssembler;
+struct ParameterBase;
+class Process;
+class ProcessVariable;
+}
 
 namespace ProcessLib
 {
@@ -24,5 +41,21 @@ std::unique_ptr<Process> createThermoMechanicsProcess(
     unsigned const integration_order,
     BaseLib::ConfigTree const& config);
 
+extern template std::unique_ptr<Process> createThermoMechanicsProcess<2>(
+    MeshLib::Mesh& mesh,
+    std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+    std::vector<ProcessVariable> const& variables,
+    std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+    unsigned const integration_order,
+    BaseLib::ConfigTree const& config);
+
+extern template std::unique_ptr<Process> createThermoMechanicsProcess<3>(
+    MeshLib::Mesh& mesh,
+    std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+    std::vector<ProcessVariable> const& variables,
+    std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+    unsigned const integration_order,
+    BaseLib::ConfigTree const& config);
+
 }  // namespace ThermoMechanics
 }  // namespace ProcessLib
diff --git a/ProcessLib/ThermoMechanics/LocalAssemblerInterface.h b/ProcessLib/ThermoMechanics/LocalAssemblerInterface.h
new file mode 100644
index 0000000000000000000000000000000000000000..d76ebedadf84378d99a3e27ebc8876b7654b10cd
--- /dev/null
+++ b/ProcessLib/ThermoMechanics/LocalAssemblerInterface.h
@@ -0,0 +1,99 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <vector>
+
+#include "NumLib/Extrapolation/ExtrapolatableElement.h"
+#include "ProcessLib/LocalAssemblerInterface.h"
+
+namespace ProcessLib
+{
+namespace ThermoMechanics
+{
+struct ThermoMechanicsLocalAssemblerInterface
+    : public ProcessLib::LocalAssemblerInterface,
+      public NumLib::ExtrapolatableElement
+{
+    virtual std::vector<double> const& getIntPtSigmaXX(
+        const double /*t*/,
+        GlobalVector const& /*current_solution*/,
+        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+        std::vector<double>& cache) const = 0;
+
+    virtual std::vector<double> const& getIntPtSigmaYY(
+        const double /*t*/,
+        GlobalVector const& /*current_solution*/,
+        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+        std::vector<double>& cache) const = 0;
+
+    virtual std::vector<double> const& getIntPtSigmaZZ(
+        const double /*t*/,
+        GlobalVector const& /*current_solution*/,
+        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+        std::vector<double>& cache) const = 0;
+
+    virtual std::vector<double> const& getIntPtSigmaXY(
+        const double /*t*/,
+        GlobalVector const& /*current_solution*/,
+        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+        std::vector<double>& cache) const = 0;
+
+    virtual std::vector<double> const& getIntPtSigmaXZ(
+        const double /*t*/,
+        GlobalVector const& /*current_solution*/,
+        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+        std::vector<double>& cache) const = 0;
+
+    virtual std::vector<double> const& getIntPtSigmaYZ(
+        const double /*t*/,
+        GlobalVector const& /*current_solution*/,
+        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+        std::vector<double>& cache) const = 0;
+
+    virtual std::vector<double> const& getIntPtEpsilonXX(
+        const double /*t*/,
+        GlobalVector const& /*current_solution*/,
+        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+        std::vector<double>& cache) const = 0;
+
+    virtual std::vector<double> const& getIntPtEpsilonYY(
+        const double /*t*/,
+        GlobalVector const& /*current_solution*/,
+        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+        std::vector<double>& cache) const = 0;
+
+    virtual std::vector<double> const& getIntPtEpsilonZZ(
+        const double /*t*/,
+        GlobalVector const& /*current_solution*/,
+        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+        std::vector<double>& cache) const = 0;
+
+    virtual std::vector<double> const& getIntPtEpsilonXY(
+        const double /*t*/,
+        GlobalVector const& /*current_solution*/,
+        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+        std::vector<double>& cache) const = 0;
+
+    virtual std::vector<double> const& getIntPtEpsilonXZ(
+        const double /*t*/,
+        GlobalVector const& /*current_solution*/,
+        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+        std::vector<double>& cache) const = 0;
+
+    virtual std::vector<double> const& getIntPtEpsilonYZ(
+        const double /*t*/,
+        GlobalVector const& /*current_solution*/,
+        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+        std::vector<double>& cache) const = 0;
+};
+
+}  // namespace ThermoMechanics
+}  // namespace ProcessLib
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h b/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h
index f01984566ecc941f2e7a2d07afb17c028e2d9dc5..b5f1abbfbc7ca93b68af676feb1ff7f974f51c13 100644
--- a/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h
@@ -13,26 +13,25 @@
 #include <vector>
 
 #include "MaterialLib/SolidModels/KelvinVector.h"
-#include "MaterialLib/SolidModels/LinearElasticIsotropic.h"
+#include "MaterialLib/SolidModels/MechanicsBase.h"
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "NumLib/Extrapolation/ExtrapolatableElement.h"
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
-#include "NumLib/Function/Interpolation.h"
 #include "ProcessLib/Deformation/BMatrixPolicy.h"
 #include "ProcessLib/Deformation/LinearBMatrix.h"
-#include "ProcessLib/LocalAssemblerInterface.h"
-#include "ProcessLib/LocalAssemblerTraits.h"
 #include "ProcessLib/Parameter/Parameter.h"
 #include "ProcessLib/Utils/InitShapeMatrices.h"
 
+#include "LocalAssemblerInterface.h"
 #include "ThermoMechanicsProcessData.h"
 
 namespace ProcessLib
 {
 namespace ThermoMechanics
 {
-template <typename BMatricesType, typename ShapeMatrixType, int DisplacementDim>
+template <typename BMatricesType, typename ShapeMatricesType,
+          int DisplacementDim>
 struct IntegrationPointData final
 {
     explicit IntegrationPointData(
@@ -43,8 +42,6 @@ struct IntegrationPointData final
     {
     }
 
-    typename ShapeMatrixType::NodalRowVectorType N;
-    typename ShapeMatrixType::GlobalDimNodalMatrixType dNdx;
     typename BMatricesType::KelvinVectorType sigma, sigma_prev;
     typename BMatricesType::KelvinVectorType eps;
     typename BMatricesType::KelvinVectorType eps_m, eps_m_prev;
@@ -55,6 +52,8 @@ struct IntegrationPointData final
         material_state_variables;
 
     double integration_weight;
+    typename ShapeMatricesType::NodalRowVectorType N;
+    typename ShapeMatricesType::GlobalDimNodalMatrixType dNdx;
 
     void pushBackState()
     {
@@ -63,31 +62,6 @@ struct IntegrationPointData final
         material_state_variables->pushBackState();
     }
 
-    static const int kelvin_vector_size =
-        KelvinVectorDimensions<DisplacementDim>::value;
-    using Invariants = MaterialLib::SolidModels::Invariants<kelvin_vector_size>;
-
-    typename BMatricesType::KelvinMatrixType updateConstitutiveRelation(
-        double const t,
-        SpatialPosition const& x_position,
-        double const dt,
-        double const linear_thermal_strain)
-    {
-        // assume isotropic thermal expansion
-        eps_m.noalias() = eps - linear_thermal_strain * Invariants::identity2;
-        auto&& solution = solid_material.integrateStress(
-            t, x_position, dt, eps_m_prev, eps_m, sigma_prev,
-            *material_state_variables);
-
-        if (!solution)
-            OGS_FATAL("Computation of local constitutive relation failed.");
-
-        KelvinMatrixType<DisplacementDim> C;
-        std::tie(sigma, material_state_variables, C) = std::move(*solution);
-
-        return C;
-    }
-
     EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
 };
 
@@ -99,83 +73,6 @@ struct SecondaryData
     std::vector<ShapeMatrixType, Eigen::aligned_allocator<ShapeMatrixType>> N;
 };
 
-struct ThermoMechanicsLocalAssemblerInterface
-    : public ProcessLib::LocalAssemblerInterface,
-      public NumLib::ExtrapolatableElement
-{
-    virtual std::vector<double> const& getIntPtSigmaXX(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& cache) const = 0;
-
-    virtual std::vector<double> const& getIntPtSigmaYY(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& cache) const = 0;
-
-    virtual std::vector<double> const& getIntPtSigmaZZ(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& cache) const = 0;
-
-    virtual std::vector<double> const& getIntPtSigmaXY(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& cache) const = 0;
-
-    virtual std::vector<double> const& getIntPtSigmaXZ(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& cache) const = 0;
-
-    virtual std::vector<double> const& getIntPtSigmaYZ(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& cache) const = 0;
-
-    virtual std::vector<double> const& getIntPtEpsilonXX(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& cache) const = 0;
-
-    virtual std::vector<double> const& getIntPtEpsilonYY(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& cache) const = 0;
-
-    virtual std::vector<double> const& getIntPtEpsilonZZ(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& cache) const = 0;
-
-    virtual std::vector<double> const& getIntPtEpsilonXY(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& cache) const = 0;
-
-    virtual std::vector<double> const& getIntPtEpsilonXZ(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& cache) const = 0;
-
-    virtual std::vector<double> const& getIntPtEpsilonYZ(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& cache) const = 0;
-};
-
 template <typename ShapeFunction, typename IntegrationMethod,
           int DisplacementDim>
 class ThermoMechanicsLocalAssembler
@@ -298,18 +195,18 @@ public:
 
         double const& dt = _process_data.dt;
 
+        unsigned const n_integration_points =
+            _integration_method.getNumberOfPoints();
+
         SpatialPosition x_position;
         x_position.setElementID(_element.getID());
 
-        unsigned const n_integration_points =
-            _integration_method.getNumberOfPoints();
         for (unsigned ip = 0; ip < n_integration_points; ip++)
         {
             x_position.setIntegrationPoint(ip);
             auto const& w = _ip_data[ip].integration_weight;
-
-            auto const& dNdx = _ip_data[ip].dNdx;
             auto const& N = _ip_data[ip].N;
+            auto const& dNdx = _ip_data[ip].dNdx;
 
             auto const x_coord =
                 interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
@@ -318,9 +215,15 @@ public:
                 DisplacementDim, ShapeFunction::NPOINTS,
                 typename BMatricesType::BMatrixType>(dNdx, N, x_coord,
                                                      _is_axially_symmetric);
-            auto const& sigma = _ip_data[ip].sigma;
+
+            auto& sigma = _ip_data[ip].sigma;
+            auto const& sigma_prev = _ip_data[ip].sigma_prev;
             auto& eps = _ip_data[ip].eps;
 
+            auto& eps_m = _ip_data[ip].eps_m;
+            auto& eps_m_prev = _ip_data[ip].eps_m_prev;
+            auto& state = _ip_data[ip].material_state_variables;
+
             double const delta_T =
                 N.dot(T) - _process_data.reference_temperature;
             // calculate thermally induced strain
@@ -334,8 +237,21 @@ public:
             // displacement equation, displacement part
             //
             eps.noalias() = B * u;
-            auto C = _ip_data[ip].updateConstitutiveRelation(
-                t, x_position, dt, linear_thermal_strain);
+
+            using Invariants =
+                MaterialLib::SolidModels::Invariants<kelvin_vector_size>;
+
+            // assume isotropic thermal expansion
+            eps_m.noalias() =
+                eps - linear_thermal_strain * Invariants::identity2;
+            auto&& solution = _ip_data[ip].solid_material.integrateStress(
+                t, x_position, dt, eps_m_prev, eps_m, sigma_prev, *state);
+
+            if (!solution)
+                OGS_FATAL("Computation of local constitutive relation failed.");
+
+            KelvinMatrixType<DisplacementDim> C;
+            std::tie(sigma, state, C) = std::move(*solution);
 
             local_Jac
                 .template block<displacement_size, displacement_size>(
@@ -354,9 +270,6 @@ public:
                        i, i * displacement_size / DisplacementDim)
                     .noalias() = N;
 
-            using Invariants =
-                MaterialLib::SolidModels::Invariants<kelvin_vector_size>;
-
             // calculate real density
             auto const rho_sr =
                 _process_data.reference_solid_density(t, x_position)[0];
@@ -367,6 +280,7 @@ public:
                 .template block<displacement_size, 1>(displacement_index, 0)
                 .noalias() -=
                 (B.transpose() * sigma - N_u.transpose() * rho_s * b) * w;
+
             //
             // displacement equation, temperature part
             //
@@ -551,6 +465,7 @@ private:
 
         return cache;
     }
+
     std::vector<double> const& getIntPtEpsilon(
         std::vector<double>& cache, std::size_t const component) const
     {
@@ -578,8 +493,8 @@ private:
 
     IntegrationMethod _integration_method;
     MeshLib::Element const& _element;
-    bool const _is_axially_symmetric;
     SecondaryData<typename ShapeMatrices::ShapeType> _secondary_data;
+    bool const _is_axially_symmetric;
 
     static const int temperature_index = 0;
     static const int temperature_size = ShapeFunction::NPOINTS;
@@ -590,29 +505,5 @@ private:
         KelvinVectorDimensions<DisplacementDim>::value;
 };
 
-template <typename ShapeFunction, typename IntegrationMethod,
-          unsigned GlobalDim, int DisplacementDim>
-class LocalAssemblerData final
-    : public ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
-                                           DisplacementDim>
-{
-public:
-    LocalAssemblerData(LocalAssemblerData const&) = delete;
-    LocalAssemblerData(LocalAssemblerData&&) = delete;
-
-    LocalAssemblerData(
-        MeshLib::Element const& e,
-        std::size_t const local_matrix_size,
-        bool is_axially_symmetric,
-        unsigned const integration_order,
-        ThermoMechanicsProcessData<DisplacementDim>& process_data)
-        : ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
-                                        DisplacementDim>(
-              e, local_matrix_size, is_axially_symmetric, integration_order,
-              process_data)
-    {
-    }
-};
-
 }  // namespace ThermoMechanics
 }  // namespace ProcessLib
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsProcess-fwd.h b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess-fwd.h
deleted file mode 100644
index fc48ba4a62a5cf5e4d6b861a091e4d34628b9265..0000000000000000000000000000000000000000
--- a/ProcessLib/ThermoMechanics/ThermoMechanicsProcess-fwd.h
+++ /dev/null
@@ -1,15 +0,0 @@
-/**
- * \copyright
- * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#pragma once
-
-#include "ThermoMechanicsProcess.h"
-
-extern template class ProcessLib::ThermoMechanics::ThermoMechanicsProcess<2>;
-extern template class ProcessLib::ThermoMechanics::ThermoMechanicsProcess<3>;
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsProcess-impl.h b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess-impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..25fc9223234b694adda3952887e932355f1f6486
--- /dev/null
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess-impl.h
@@ -0,0 +1,204 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <cassert>
+
+#include "BaseLib/Functional.h"
+#include "ProcessLib/SmallDeformation/CreateLocalAssemblers.h"
+
+#include "ThermoMechanicsFEM.h"
+
+namespace ProcessLib
+{
+namespace ThermoMechanics
+{
+template <int DisplacementDim>
+ThermoMechanicsProcess<DisplacementDim>::ThermoMechanicsProcess(
+    MeshLib::Mesh& mesh,
+    std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+    std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+    unsigned const integration_order,
+    std::vector<std::reference_wrapper<ProcessVariable>>&& process_variables,
+    ThermoMechanicsProcessData<DisplacementDim>&& process_data,
+    SecondaryVariableCollection&& secondary_variables,
+    NumLib::NamedFunctionCaller&& named_function_caller)
+    : Process(mesh, std::move(jacobian_assembler), parameters,
+              integration_order, std::move(process_variables),
+              std::move(secondary_variables), std::move(named_function_caller)),
+      _process_data(std::move(process_data))
+{
+}
+
+template <int DisplacementDim>
+bool ThermoMechanicsProcess<DisplacementDim>::isLinear() const
+{
+    return false;
+}
+
+template <int DisplacementDim>
+void ThermoMechanicsProcess<DisplacementDim>::initializeConcreteProcess(
+    NumLib::LocalToGlobalIndexMap const& dof_table,
+    MeshLib::Mesh const& mesh,
+    unsigned const integration_order)
+{
+    ProcessLib::SmallDeformation::createLocalAssemblers<
+        DisplacementDim, ThermoMechanicsLocalAssembler>(
+        mesh.getElements(), dof_table, _local_assemblers,
+        mesh.isAxiallySymmetric(), integration_order, _process_data);
+
+    // TODO move the two data members somewhere else.
+    // for extrapolation of secondary variables
+    std::vector<MeshLib::MeshSubsets> all_mesh_subsets_single_component;
+    all_mesh_subsets_single_component.emplace_back(
+        _mesh_subset_all_nodes.get());
+    _local_to_global_index_map_single_component.reset(
+        new NumLib::LocalToGlobalIndexMap(
+            std::move(all_mesh_subsets_single_component),
+            // by location order is needed for output
+            NumLib::ComponentOrder::BY_LOCATION));
+
+    Base::_secondary_variables.addSecondaryVariable(
+        "sigma_xx",
+        makeExtrapolator(
+            1, getExtrapolator(), _local_assemblers,
+            &ThermoMechanicsLocalAssemblerInterface::getIntPtSigmaXX));
+
+    Base::_secondary_variables.addSecondaryVariable(
+        "sigma_yy",
+        makeExtrapolator(
+            1, getExtrapolator(), _local_assemblers,
+            &ThermoMechanicsLocalAssemblerInterface::getIntPtSigmaYY));
+
+    Base::_secondary_variables.addSecondaryVariable(
+        "sigma_zz",
+        makeExtrapolator(
+            1, getExtrapolator(), _local_assemblers,
+            &ThermoMechanicsLocalAssemblerInterface::getIntPtSigmaZZ));
+
+    Base::_secondary_variables.addSecondaryVariable(
+        "sigma_xy",
+        makeExtrapolator(
+            1, getExtrapolator(), _local_assemblers,
+            &ThermoMechanicsLocalAssemblerInterface::getIntPtSigmaXY));
+
+    if (DisplacementDim == 3)
+    {
+        Base::_secondary_variables.addSecondaryVariable(
+            "sigma_xz",
+            makeExtrapolator(
+                1, getExtrapolator(), _local_assemblers,
+                &ThermoMechanicsLocalAssemblerInterface::getIntPtSigmaXZ));
+
+        Base::_secondary_variables.addSecondaryVariable(
+            "sigma_yz",
+            makeExtrapolator(
+                1, getExtrapolator(), _local_assemblers,
+                &ThermoMechanicsLocalAssemblerInterface::getIntPtSigmaYZ));
+    }
+    Base::_secondary_variables.addSecondaryVariable(
+        "epsilon_xx",
+        makeExtrapolator(
+            1, getExtrapolator(), _local_assemblers,
+            &ThermoMechanicsLocalAssemblerInterface::getIntPtEpsilonXX));
+
+    Base::_secondary_variables.addSecondaryVariable(
+        "epsilon_yy",
+        makeExtrapolator(
+            1, getExtrapolator(), _local_assemblers,
+            &ThermoMechanicsLocalAssemblerInterface::getIntPtEpsilonYY));
+
+    Base::_secondary_variables.addSecondaryVariable(
+        "epsilon_zz",
+        makeExtrapolator(
+            1, getExtrapolator(), _local_assemblers,
+            &ThermoMechanicsLocalAssemblerInterface::getIntPtEpsilonZZ));
+
+    Base::_secondary_variables.addSecondaryVariable(
+        "epsilon_xy",
+        makeExtrapolator(
+            1, getExtrapolator(), _local_assemblers,
+            &ThermoMechanicsLocalAssemblerInterface::getIntPtEpsilonXY));
+    if (DisplacementDim == 3)
+    {
+        Base::_secondary_variables.addSecondaryVariable(
+            "epsilon_yz",
+            makeExtrapolator(
+                1, getExtrapolator(), _local_assemblers,
+                &ThermoMechanicsLocalAssemblerInterface::getIntPtEpsilonYZ));
+
+        Base::_secondary_variables.addSecondaryVariable(
+            "epsilon_xz",
+            makeExtrapolator(
+                1, getExtrapolator(), _local_assemblers,
+                &ThermoMechanicsLocalAssemblerInterface::getIntPtEpsilonXZ));
+    }
+}
+
+template <int DisplacementDim>
+void ThermoMechanicsProcess<DisplacementDim>::assembleConcreteProcess(
+    const double t, GlobalVector const& x, GlobalMatrix& M, GlobalMatrix& K,
+    GlobalVector& b)
+{
+    DBUG("Assemble ThermoMechanicsProcess.");
+
+    // Call global assembler for each local assembly item.
+    GlobalExecutor::executeMemberDereferenced(
+        _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
+        *_local_to_global_index_map, t, x, M, K, b, _coupling_term);
+}
+
+template <int DisplacementDim>
+void ThermoMechanicsProcess<DisplacementDim>::
+    assembleWithJacobianConcreteProcess(const double t, GlobalVector const& x,
+                                        GlobalVector const& xdot,
+                                        const double dxdot_dx,
+                                        const double dx_dx, GlobalMatrix& M,
+                                        GlobalMatrix& K, GlobalVector& b,
+                                        GlobalMatrix& Jac)
+{
+    DBUG("AssembleJacobian ThermoMechanicsProcess.");
+
+    // Call global assembler for each local assembly item.
+    GlobalExecutor::executeMemberDereferenced(
+        _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
+        _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
+        dx_dx, M, K, b, Jac, _coupling_term);
+}
+
+template <int DisplacementDim>
+void ThermoMechanicsProcess<DisplacementDim>::preTimestepConcreteProcess(
+    GlobalVector const& x, double const t, double const dt)
+{
+    DBUG("PreTimestep ThermoMechanicsProcess.");
+
+    _process_data.dt = dt;
+    _process_data.t = t;
+
+    GlobalExecutor::executeMemberOnDereferenced(
+        &ThermoMechanicsLocalAssemblerInterface::preTimestep, _local_assemblers,
+        *_local_to_global_index_map, x, t, dt);
+}
+
+template <int DisplacementDim>
+void ThermoMechanicsProcess<DisplacementDim>::postTimestepConcreteProcess(
+    GlobalVector const& x)
+{
+    DBUG("PostTimestep ThermoMechanicsProcess.");
+
+    GlobalExecutor::executeMemberOnDereferenced(
+        &ThermoMechanicsLocalAssemblerInterface::postTimestep,
+        _local_assemblers, *_local_to_global_index_map, x);
+}
+
+}  // namespace ThermoMechanics
+}  // namespace ProcessLib
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.cpp b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.cpp
index 55fddd02e3d5aec27f3fbab725e3ab761011e976..9b391011d73a487712f100863f7cc5aecbdca69d 100644
--- a/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.cpp
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.cpp
@@ -8,7 +8,7 @@
  */
 
 #include "ThermoMechanicsProcess.h"
-#include "ThermoMechanicsProcess-fwd.h"
+#include "ThermoMechanicsProcess-impl.h"
 
 namespace ProcessLib
 {
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.h b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.h
index e0d9339821f2e73f32bec08410e0d0bf6578fdd5..85d5116103d78470acca18e62f80ac4dc7b2299e 100644
--- a/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.h
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.h
@@ -9,19 +9,17 @@
 
 #pragma once
 
-#include <cassert>
-
-#include "NumLib/DOF/DOFTableUtil.h"
 #include "ProcessLib/Process.h"
-#include "ProcessLib/SmallDeformation/CreateLocalAssemblers.h"
 
-#include "ThermoMechanicsFEM.h"
+#include "LocalAssemblerInterface.h"
 #include "ThermoMechanicsProcessData.h"
 
 namespace ProcessLib
 {
 namespace ThermoMechanics
 {
+struct ThermoMechanicsLocalAssemblerInterface;
+
 template <int DisplacementDim>
 class ThermoMechanicsProcess final : public Process
 {
@@ -38,181 +36,47 @@ public:
             process_variables,
         ThermoMechanicsProcessData<DisplacementDim>&& process_data,
         SecondaryVariableCollection&& secondary_variables,
-        NumLib::NamedFunctionCaller&& named_function_caller)
-        : Process(mesh, std::move(jacobian_assembler), parameters,
-                  integration_order, std::move(process_variables),
-                  std::move(secondary_variables),
-                  std::move(named_function_caller)),
-          _process_data(std::move(process_data))
-    {
-    }
+        NumLib::NamedFunctionCaller&& named_function_caller);
 
     //! \name ODESystem interface
     //! @{
-
-    bool isLinear() const override { return false; }
+    bool isLinear() const override;
     //! @}
 
 private:
-    using LocalAssemblerInterface = ThermoMechanicsLocalAssemblerInterface;
-
     void initializeConcreteProcess(
         NumLib::LocalToGlobalIndexMap const& dof_table,
         MeshLib::Mesh const& mesh,
-        unsigned const integration_order) override
-    {
-        ProcessLib::SmallDeformation::createLocalAssemblers<
-        DisplacementDim, ThermoMechanicsLocalAssembler>(
-        mesh.getElements(), dof_table, _local_assemblers,
-        mesh.isAxiallySymmetric(), integration_order, _process_data);
-
-        // TODO move the two data members somewhere else.
-        // for extrapolation of secondary variables
-        std::vector<MeshLib::MeshSubsets> all_mesh_subsets_single_component;
-        all_mesh_subsets_single_component.emplace_back(
-            _mesh_subset_all_nodes.get());
-        _local_to_global_index_map_single_component.reset(
-            new NumLib::LocalToGlobalIndexMap(
-                std::move(all_mesh_subsets_single_component),
-                // by location order is needed for output
-                NumLib::ComponentOrder::BY_LOCATION));
-
-        Base::_secondary_variables.addSecondaryVariable(
-            "sigma_xx",
-            makeExtrapolator(
-                1, getExtrapolator(), _local_assemblers,
-                &ThermoMechanicsLocalAssemblerInterface::getIntPtSigmaXX));
-
-        Base::_secondary_variables.addSecondaryVariable(
-            "sigma_yy",
-            makeExtrapolator(
-                1, getExtrapolator(), _local_assemblers,
-                &ThermoMechanicsLocalAssemblerInterface::getIntPtSigmaYY));
-
-        Base::_secondary_variables.addSecondaryVariable(
-            "sigma_zz",
-            makeExtrapolator(
-                1, getExtrapolator(), _local_assemblers,
-                &ThermoMechanicsLocalAssemblerInterface::getIntPtSigmaZZ));
-
-        Base::_secondary_variables.addSecondaryVariable(
-            "sigma_xy",
-            makeExtrapolator(
-                1, getExtrapolator(), _local_assemblers,
-                &ThermoMechanicsLocalAssemblerInterface::getIntPtSigmaXY));
-
-        if (DisplacementDim == 3)
-        {
-            Base::_secondary_variables.addSecondaryVariable(
-                "sigma_xz",
-                makeExtrapolator(
-                    1, getExtrapolator(), _local_assemblers,
-                    &ThermoMechanicsLocalAssemblerInterface::getIntPtSigmaXZ));
-
-            Base::_secondary_variables.addSecondaryVariable(
-                "sigma_yz",
-                makeExtrapolator(
-                    1, getExtrapolator(), _local_assemblers,
-                    &ThermoMechanicsLocalAssemblerInterface::getIntPtSigmaYZ));
-        }
-        Base::_secondary_variables.addSecondaryVariable(
-            "epsilon_xx",
-            makeExtrapolator(
-                1, getExtrapolator(), _local_assemblers,
-                &ThermoMechanicsLocalAssemblerInterface::getIntPtEpsilonXX));
-
-        Base::_secondary_variables.addSecondaryVariable(
-            "epsilon_yy",
-            makeExtrapolator(
-                1, getExtrapolator(), _local_assemblers,
-                &ThermoMechanicsLocalAssemblerInterface::getIntPtEpsilonYY));
-
-        Base::_secondary_variables.addSecondaryVariable(
-            "epsilon_zz",
-            makeExtrapolator(
-                1, getExtrapolator(), _local_assemblers,
-                &ThermoMechanicsLocalAssemblerInterface::getIntPtEpsilonZZ));
-
-        Base::_secondary_variables.addSecondaryVariable(
-            "epsilon_xy",
-            makeExtrapolator(
-                1, getExtrapolator(), _local_assemblers,
-                &ThermoMechanicsLocalAssemblerInterface::getIntPtEpsilonXY));
-        if (DisplacementDim == 3)
-        {
-            Base::_secondary_variables.addSecondaryVariable(
-                "epsilon_yz",
-                makeExtrapolator(1, getExtrapolator(), _local_assemblers,
-                                 &ThermoMechanicsLocalAssemblerInterface::
-                                     getIntPtEpsilonYZ));
-
-            Base::_secondary_variables.addSecondaryVariable(
-                "epsilon_xz",
-                makeExtrapolator(1, getExtrapolator(), _local_assemblers,
-                                 &ThermoMechanicsLocalAssemblerInterface::
-                                     getIntPtEpsilonXZ));
-        }
-    }
+        unsigned const integration_order) override;
 
     void assembleConcreteProcess(
         const double t, GlobalVector const& x, GlobalMatrix& M, GlobalMatrix& K,
-        GlobalVector& b) override
-    {
-        DBUG("Assemble ThermoMechanicsProcess.");
-
-        // Call global assembler for each local assembly item.
-        GlobalExecutor::executeMemberDereferenced(
-            _global_assembler, &VectorMatrixAssembler::assemble,
-            _local_assemblers, *_local_to_global_index_map, t, x, M, K, b,
-            _coupling_term);
-    }
+        GlobalVector& b) override;
 
     void assembleWithJacobianConcreteProcess(
         const double t, GlobalVector const& x, GlobalVector const& xdot,
         const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
-        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override
-    {
-        DBUG("AssembleJacobian ThermoMechanicsProcess.");
-
-        // Call global assembler for each local assembly item.
-        GlobalExecutor::executeMemberDereferenced(
-            _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
-            _local_assemblers, *_local_to_global_index_map, t, x, xdot,
-            dxdot_dx, dx_dx, M, K, b, Jac, _coupling_term);
-    }
+        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override;
 
     void preTimestepConcreteProcess(GlobalVector const& x, double const t,
-                                    double const dt) override
-    {
-        DBUG("PreTimestep ThermoMechanicsProcess.");
-
-        _process_data.dt = dt;
-        _process_data.t = t;
-
-        GlobalExecutor::executeMemberOnDereferenced(
-            &ThermoMechanicsLocalAssemblerInterface::preTimestep,
-            _local_assemblers, *_local_to_global_index_map, x, t, dt);
-    }
+                                    double const dt) override;
 
-    void postTimestepConcreteProcess(GlobalVector const& x) override
-    {
-        DBUG("PostTimestep ThermoMechanicsProcess.");
-
-        GlobalExecutor::executeMemberOnDereferenced(
-            &ThermoMechanicsLocalAssemblerInterface::postTimestep,
-            _local_assemblers, *_local_to_global_index_map, x);
-    }
+    void postTimestepConcreteProcess(GlobalVector const& x) override;
 
 private:
     std::vector<MeshLib::Node*> _base_nodes;
     std::unique_ptr<MeshLib::MeshSubset const> _mesh_subset_base_nodes;
     ThermoMechanicsProcessData<DisplacementDim> _process_data;
 
-    std::vector<std::unique_ptr<LocalAssemblerInterface>> _local_assemblers;
+    std::vector<std::unique_ptr<ThermoMechanicsLocalAssemblerInterface>>
+        _local_assemblers;
 
     std::unique_ptr<NumLib::LocalToGlobalIndexMap>
         _local_to_global_index_map_single_component;
 };
 
+extern template class ThermoMechanicsProcess<2>;
+extern template class ThermoMechanicsProcess<3>;
+
 }  // namespace ThermoMechanics
 }  // namespace ProcessLib
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsProcessData.h b/ProcessLib/ThermoMechanics/ThermoMechanicsProcessData.h
index 5ac1bb4210a7dedd2caa2c0c9bfa9a062b8e6e1b..b59fc749ac87f3668be6d1b646dae93aced478cb 100644
--- a/ProcessLib/ThermoMechanics/ThermoMechanicsProcessData.h
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsProcessData.h
@@ -9,11 +9,14 @@
 
 #pragma once
 
-namespace MeshLib
+namespace MaterialLib
 {
-class Element;
+namespace Solids
+{
+template <int DisplacementDim>
+struct MechanicsBase;
+}
 }
-
 namespace ProcessLib
 {
 namespace ThermoMechanics
@@ -73,8 +76,8 @@ struct ThermoMechanicsProcessData
         thermal_conductivity;  // TODO To be changed as a matrix type variable.
     double const reference_temperature;
     Eigen::Matrix<double, DisplacementDim, 1> const specific_body_force;
-    double dt;
-    double t;
+    double dt = 0;
+    double t = 0;
 
     EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
 };
diff --git a/Tests/Data b/Tests/Data
index 6bbe61eb4d6112f89d5e60ae7967fa806435f5ab..a90755d46a01be845056ea521df0c7f9e203ccfd 160000
--- a/Tests/Data
+++ b/Tests/Data
@@ -1 +1 @@
-Subproject commit 6bbe61eb4d6112f89d5e60ae7967fa806435f5ab
+Subproject commit a90755d46a01be845056ea521df0c7f9e203ccfd