diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index 58be51290e959030b4a0c6177ceee7f8a90d01f9..688784055cec5cc9e93ce20601aa278e0cb3f2a0 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/Applications/ApplicationsLib/ProjectData.cpp
@@ -96,6 +96,9 @@
 #ifdef OGS_BUILD_PROCESS_THERMALTWOPHASEFLOWWITHPP
 #include "ProcessLib/ThermalTwoPhaseFlowWithPP/CreateThermalTwoPhaseFlowWithPPProcess.h"
 #endif
+#ifdef OGS_BUILD_PROCESS_THERMOHYDROMECHANICS
+#include "ProcessLib/ThermoHydroMechanics/CreateThermoHydroMechanicsProcess.h"
+#endif
 #ifdef OGS_BUILD_PROCESS_THERMOMECHANICALPHASEFIELD
 #include "ProcessLib/ThermoMechanicalPhaseField/CreateThermoMechanicalPhaseFieldProcess.h"
 #endif
@@ -741,6 +744,34 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config,
         }
         else
 #endif
+#ifdef OGS_BUILD_PROCESS_THERMOHYDROMECHANICS
+            if (type == "THERMO_HYDRO_MECHANICS")
+        {
+            //! \ogs_file_param{prj__processes__process__THERMO_HYDRO_MECHANICS__dimension}
+            switch (process_config.getConfigParameter<int>("dimension"))
+            {
+                case 2:
+                    process = ProcessLib::ThermoHydroMechanics::
+                        createThermoHydroMechanicsProcess<2>(
+                            *_mesh_vec[0], std::move(jacobian_assembler),
+                            _process_variables, _parameters, integration_order,
+                            process_config);
+                    break;
+                case 3:
+                    process = ProcessLib::ThermoHydroMechanics::
+                        createThermoHydroMechanicsProcess<3>(
+                            *_mesh_vec[0], std::move(jacobian_assembler),
+                            _process_variables, _parameters, integration_order,
+                            process_config);
+                    break;
+                default:
+                    OGS_FATAL(
+                        "THERMO_HYDRO_MECHANICS process does not support given "
+                        "dimension");
+            }
+        }
+        else
+#endif
 #ifdef OGS_BUILD_PROCESS_THERMOMECHANICALPHASEFIELD
             if (type == "THERMO_MECHANICAL_PHASE_FIELD")
         {
diff --git a/CMakeLists.txt b/CMakeLists.txt
index dd233617914b79438a9c68835da08db7ad0750fd..8e5312921ebe1b18d2530291a1387301e2abb6f6 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -118,6 +118,7 @@ set(ProcessesList
     SmallDeformationNonlocal
     TES
     ThermalTwoPhaseFlowWithPP
+    ThermoHydroMechanics
     ThermoMechanicalPhaseField
     ThermoMechanics
     TwoPhaseFlowWithPP
diff --git a/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/c_THERMO_HYDRO_MECHANICS.md b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/c_THERMO_HYDRO_MECHANICS.md
new file mode 100644
index 0000000000000000000000000000000000000000..11e90e8358e39c914965cc6fc2fc121a5febb1f9
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/c_THERMO_HYDRO_MECHANICS.md
@@ -0,0 +1 @@
+Coupled thermo-hydro-mechanical process. It is implemented monolithically.
diff --git a/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/process_variables/i_process_variables.md b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/process_variables/i_process_variables.md
new file mode 100644
index 0000000000000000000000000000000000000000..7b74dddccb1a3ae3ba02e727ca832c803293dd13
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/process_variables/i_process_variables.md
@@ -0,0 +1 @@
+The process variables for displacement, pressure and temperature.
diff --git a/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/process_variables/t_displacement.md b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/process_variables/t_displacement.md
new file mode 100644
index 0000000000000000000000000000000000000000..b6d5f29924d468888547d4a5f9eef957f496c26c
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/process_variables/t_displacement.md
@@ -0,0 +1 @@
+Process variable name for displacement.
diff --git a/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/process_variables/t_pressure.md b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/process_variables/t_pressure.md
new file mode 100644
index 0000000000000000000000000000000000000000..36a6d4b8b7a6860e66353bd83b152a2cae3de678
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/process_variables/t_pressure.md
@@ -0,0 +1 @@
+Process variable name for pressure.
diff --git a/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/process_variables/t_temperature.md b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/process_variables/t_temperature.md
new file mode 100644
index 0000000000000000000000000000000000000000..f6efd419ab89f60030e6f044c02c53fd9cc44763
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/process_variables/t_temperature.md
@@ -0,0 +1 @@
+Process variable name for temperature.
diff --git a/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_biot_coefficient.md b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_biot_coefficient.md
new file mode 100644
index 0000000000000000000000000000000000000000..2db4bbf026703aea93279258976457595f09ae97
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_biot_coefficient.md
@@ -0,0 +1 @@
+Biot coefficient.
diff --git a/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_coupling_scheme.md b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_coupling_scheme.md
new file mode 100644
index 0000000000000000000000000000000000000000..703141a127a49902cdee92ef1400d2f9496e6223
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_coupling_scheme.md
@@ -0,0 +1,2 @@
+An optional input to select the coupling scheme. So far, only the monolithic
+scheme is available for THERMO_HYDRO_MECHANICS, and this input can be omitted.
diff --git a/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_dimension.md b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_dimension.md
new file mode 100644
index 0000000000000000000000000000000000000000..fe666ef74ad120b10ddadd3f2ee29a0b6b8f8dbc
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_dimension.md
@@ -0,0 +1 @@
+Displacement vector dimension. The displacement dimension must be equal to the mesh dimension.
diff --git a/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_fluid_density.md b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_fluid_density.md
new file mode 100644
index 0000000000000000000000000000000000000000..722affa8a460db56881b8daeb0831cea20f5ac62
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_fluid_density.md
@@ -0,0 +1 @@
+\copydoc ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsProcessData::fluid_density
diff --git a/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_fluid_specific_heat_capacity.md b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_fluid_specific_heat_capacity.md
new file mode 100644
index 0000000000000000000000000000000000000000..024bf06abff1e70b276717295a6360e4cbb3ace1
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_fluid_specific_heat_capacity.md
@@ -0,0 +1 @@
+Fluid specific heat capacity is the heat capacity per unit mass of the fluid.
diff --git a/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_fluid_thermal_conductivity.md b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_fluid_thermal_conductivity.md
new file mode 100644
index 0000000000000000000000000000000000000000..75fd9758ed5dcaf750ef28a4568f5f1e09203711
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_fluid_thermal_conductivity.md
@@ -0,0 +1 @@
+Fluid thermal conductivity, the property of fluid to conduct heat.
diff --git a/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_fluid_viscosity.md b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_fluid_viscosity.md
new file mode 100644
index 0000000000000000000000000000000000000000..f49c37cdbfa92dc9d521c138d352ae4508d33a0a
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_fluid_viscosity.md
@@ -0,0 +1 @@
+\copydoc ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsProcessData::fluid_viscosity
diff --git a/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_fluid_volumetric_thermal_expansion_coefficient.md b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_fluid_volumetric_thermal_expansion_coefficient.md
new file mode 100644
index 0000000000000000000000000000000000000000..440f6d0f3143fb2b63caa68149439a5a50675d75
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_fluid_volumetric_thermal_expansion_coefficient.md
@@ -0,0 +1 @@
+Volumetric thermal expansion coefficient of fluid.
diff --git a/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_intrinsic_permeability.md b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_intrinsic_permeability.md
new file mode 100644
index 0000000000000000000000000000000000000000..9559792cd8011dfa8aa950a1d8ffa12ec0d02c9d
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_intrinsic_permeability.md
@@ -0,0 +1 @@
+Permeability of porous medium.
diff --git a/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_porosity.md b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_porosity.md
new file mode 100644
index 0000000000000000000000000000000000000000..418fdcc6bc1a649c41c87583b5d458aff39537a8
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_porosity.md
@@ -0,0 +1 @@
+Porosity of porous medium.
diff --git a/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_reference_temperature.md b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_reference_temperature.md
new file mode 100644
index 0000000000000000000000000000000000000000..6951343d3859fe36e8b0b73e348ebcbf83ca54b6
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_reference_temperature.md
@@ -0,0 +1 @@
+Reference temperature for the thermal part of thermo-hydro-mechanical process.
diff --git a/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_solid_density.md b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_solid_density.md
new file mode 100644
index 0000000000000000000000000000000000000000..197239040b301996ba1eb512365643e2e2880a76
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_solid_density.md
@@ -0,0 +1 @@
+Solid density.
diff --git a/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_solid_linear_thermal_expansion_coefficient.md b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_solid_linear_thermal_expansion_coefficient.md
new file mode 100644
index 0000000000000000000000000000000000000000..38eaff68a74714b0e6f9f31ad19aceb4cae7ba0f
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_solid_linear_thermal_expansion_coefficient.md
@@ -0,0 +1 @@
+Solid linear expansion in one dimension.
diff --git a/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_solid_specific_heat_capacity.md b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_solid_specific_heat_capacity.md
new file mode 100644
index 0000000000000000000000000000000000000000..1696adff39dc64b2f079078467745b99684bd5aa
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_solid_specific_heat_capacity.md
@@ -0,0 +1 @@
+Specific heat capacity, the heat capacity per unit mass of a material.
diff --git a/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_solid_thermal_conductivity.md b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_solid_thermal_conductivity.md
new file mode 100644
index 0000000000000000000000000000000000000000..0d0a26a12127fe7b1ea36565ba9f373e9e0788b2
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_solid_thermal_conductivity.md
@@ -0,0 +1 @@
+Solid thermal conductivity, the property of solid to conduct heat.
diff --git a/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_specific_body_force.md b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_specific_body_force.md
new file mode 100644
index 0000000000000000000000000000000000000000..b6b5fbd2f1431bee0240db9faade96613b41c726
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_specific_body_force.md
@@ -0,0 +1 @@
+Specific body forces that are used to apply gravitational forces, a vector of size equal to the displacement dimension.
\ No newline at end of file
diff --git a/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_specific_storage.md b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_specific_storage.md
new file mode 100644
index 0000000000000000000000000000000000000000..af55c27a7ae0e332abf290dd82c7621c91c6603b
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/THERMO_HYDRO_MECHANICS/t_specific_storage.md
@@ -0,0 +1 @@
+Volumetric average specific storage of porous medium and fluid phases.
diff --git a/ProcessLib/ThermoHydroMechanics/CMakeLists.txt b/ProcessLib/ThermoHydroMechanics/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..4abb37657ef5fad3d72afe1842f0d6cf245b6822
--- /dev/null
+++ b/ProcessLib/ThermoHydroMechanics/CMakeLists.txt
@@ -0,0 +1,12 @@
+APPEND_SOURCE_FILES(SOURCES)
+
+add_library(ThermoHydroMechanics ${SOURCES})
+if(BUILD_SHARED_LIBS)
+    install(TARGETS ThermoHydroMechanics LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR})
+endif()
+target_link_libraries(ThermoHydroMechanics
+    PUBLIC ProcessLib
+    PRIVATE ParameterLib
+)
+
+include(Tests.cmake)
diff --git a/ProcessLib/ThermoHydroMechanics/CreateLocalAssemblers.h b/ProcessLib/ThermoHydroMechanics/CreateLocalAssemblers.h
new file mode 100644
index 0000000000000000000000000000000000000000..84b71bfe230345102e298fd67a8e202295d31d91
--- /dev/null
+++ b/ProcessLib/ThermoHydroMechanics/CreateLocalAssemblers.h
@@ -0,0 +1,88 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2019, 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 <logog/include/logog.hpp>
+
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+
+#include "LocalDataInitializer.h"
+
+namespace ProcessLib
+{
+namespace ThermoHydroMechanics
+{
+namespace detail
+{
+template <int GlobalDim,
+          template <typename, typename, typename, int>
+          class LocalAssemblerImplementation,
+          typename LocalAssemblerInterface, typename... ExtraCtorArgs>
+void createLocalAssemblers(
+    NumLib::LocalToGlobalIndexMap const& dof_table,
+    const unsigned shapefunction_order,
+    std::vector<MeshLib::Element*> const& mesh_elements,
+    std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
+    ExtraCtorArgs&&... extra_ctor_args)
+{
+    // Shape matrices initializer
+    using LocalDataInitializer =
+        LocalDataInitializer<LocalAssemblerInterface,
+                             LocalAssemblerImplementation, GlobalDim,
+                             ExtraCtorArgs...>;
+
+    DBUG("Create local assemblers.");
+    // Populate the vector of local assemblers.
+    local_assemblers.resize(mesh_elements.size());
+
+    LocalDataInitializer initializer(dof_table, shapefunction_order);
+
+    DBUG("Calling local assembler builder for all mesh elements.");
+    GlobalExecutor::transformDereferenced(
+        initializer, mesh_elements, local_assemblers,
+        std::forward<ExtraCtorArgs>(extra_ctor_args)...);
+}
+
+}  // namespace detail
+
+/*! Creates local assemblers for each element of the given \c mesh.
+ *
+ * \tparam LocalAssemblerImplementation the individual local assembler type
+ * \tparam LocalAssemblerInterface the general local assembler interface
+ * \tparam ExtraCtorArgs types of additional constructor arguments.
+ *         Those arguments will be passed to the constructor of
+ *         \c LocalAssemblerImplementation.
+ *
+ * The first two template parameters cannot be deduced from the arguments.
+ * Therefore they always have to be provided manually.
+ */
+template <int GlobalDim,
+          template <typename, typename, typename, int>
+          class LocalAssemblerImplementation,
+          typename LocalAssemblerInterface, typename... ExtraCtorArgs>
+void createLocalAssemblers(
+    const unsigned /*dimension*/,
+    std::vector<MeshLib::Element*> const& mesh_elements,
+    NumLib::LocalToGlobalIndexMap const& dof_table,
+    const unsigned shapefunction_order,
+    std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
+    ExtraCtorArgs&&... extra_ctor_args)
+{
+    DBUG("Create local assemblers.");
+
+    detail::createLocalAssemblers<GlobalDim, LocalAssemblerImplementation>(
+        dof_table, shapefunction_order, mesh_elements, local_assemblers,
+        std::forward<ExtraCtorArgs>(extra_ctor_args)...);
+}
+}  // namespace ThermoHydroMechanics
+
+}  // namespace ProcessLib
diff --git a/ProcessLib/ThermoHydroMechanics/CreateThermoHydroMechanicsProcess.cpp b/ProcessLib/ThermoHydroMechanics/CreateThermoHydroMechanicsProcess.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9339108fb630a013d5b4056b1cf0fa36ecccb0e9
--- /dev/null
+++ b/ProcessLib/ThermoHydroMechanics/CreateThermoHydroMechanicsProcess.cpp
@@ -0,0 +1,312 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "CreateThermoHydroMechanicsProcess.h"
+
+#include <cassert>
+
+#include "MaterialLib/SolidModels/CreateConstitutiveRelation.h"
+#include "MaterialLib/SolidModels/MechanicsBase.h"
+#include "ParameterLib/Utils.h"
+#include "ProcessLib/Output/CreateSecondaryVariables.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
+
+#include "ThermoHydroMechanicsProcess.h"
+#include "ThermoHydroMechanicsProcessData.h"
+
+namespace ProcessLib
+{
+namespace ThermoHydroMechanics
+{
+template <int DisplacementDim>
+std::unique_ptr<Process> createThermoHydroMechanicsProcess(
+    MeshLib::Mesh& mesh,
+    std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+    std::vector<ProcessVariable> const& variables,
+    std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
+    unsigned const integration_order,
+    BaseLib::ConfigTree const& config)
+{
+    //! \ogs_file_param{prj__processes__process__type}
+    config.checkConfigParameter("type", "THERMO_HYDRO_MECHANICS");
+    DBUG("Create ThermoHydroMechanicsProcess.");
+
+    auto const staggered_scheme =
+        //! \ogs_file_param{prj__processes__process__THERMO_HYDRO_MECHANICS__coupling_scheme}
+        config.getConfigParameterOptional<std::string>("coupling_scheme");
+    const bool use_monolithic_scheme =
+        !(staggered_scheme && (*staggered_scheme == "staggered"));
+
+    // Process variable.
+
+    //! \ogs_file_param{prj__processes__process__THERMO_HYDRO_MECHANICS__process_variables}
+    auto const pv_config = config.getConfigSubtree("process_variables");
+
+    ProcessVariable* variable_T;
+    ProcessVariable* variable_p;
+    ProcessVariable* variable_u;
+    std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>
+        process_variables;
+    if (use_monolithic_scheme)  // monolithic scheme.
+    {
+        auto per_process_variables = findProcessVariables(
+            variables, pv_config,
+            {//! \ogs_file_param_special{prj__processes__process__THERMO_HYDRO_MECHANICS__process_variables__temperature}
+             "temperature",
+             //! \ogs_file_param_special{prj__processes__process__THERMO_HYDRO_MECHANICS__process_variables__pressure}
+             "pressure",
+             //! \ogs_file_param_special{prj__processes__process__THERMO_HYDRO_MECHANICS__process_variables__displacement}
+             "displacement"});
+        variable_T = &per_process_variables[0].get();
+        variable_p = &per_process_variables[1].get();
+        variable_u = &per_process_variables[2].get();
+        process_variables.push_back(std::move(per_process_variables));
+    }
+    else  // staggered scheme.
+    {
+        using namespace std::string_literals;
+        for (auto const& variable_name :
+             {"temperature"s, "pressure"s, "displacement"s})
+        {
+            auto per_process_variables =
+                findProcessVariables(variables, pv_config, {variable_name});
+            process_variables.push_back(std::move(per_process_variables));
+        }
+        variable_T = &process_variables[0][0].get();
+        variable_p = &process_variables[1][0].get();
+        variable_u = &process_variables[2][0].get();
+    }
+
+    DBUG("Associate displacement with process variable '%s'.",
+         variable_u->getName().c_str());
+
+    if (variable_u->getNumberOfComponents() != DisplacementDim)
+    {
+        OGS_FATAL(
+            "Number of components of the process variable '%s' is different "
+            "from the displacement dimension: got %d, expected %d",
+            variable_u->getName().c_str(),
+            variable_u->getNumberOfComponents(),
+            DisplacementDim);
+    }
+
+    DBUG("Associate pressure with process variable '%s'.",
+         variable_p->getName().c_str());
+    if (variable_p->getNumberOfComponents() != 1)
+    {
+        OGS_FATAL(
+            "Pressure process variable '%s' is not a scalar variable but has "
+            "%d components.",
+            variable_p->getName().c_str(),
+            variable_p->getNumberOfComponents());
+    }
+
+    DBUG("Associate temperature with process variable '%s'.",
+         variable_T->getName().c_str());
+    if (variable_T->getNumberOfComponents() != 1)
+    {
+        OGS_FATAL(
+            "temperature process variable '%s' is not a scalar variable but "
+            "has %d components.",
+            variable_T->getName().c_str(),
+            variable_T->getNumberOfComponents());
+    }
+
+    auto solid_constitutive_relations =
+        MaterialLib::Solids::createConstitutiveRelations<DisplacementDim>(
+            parameters, config);
+
+    // Intrinsic permeability (only one scalar per element, i.e. the isotropic
+    // case is handled at the moment)
+    auto& intrinsic_permeability = ParameterLib::findParameter<double>(
+        config,
+        //! \ogs_file_param_special{prj__processes__process__THERMO_HYDRO_MECHANICS__intrinsic_permeability}
+        "intrinsic_permeability", parameters, 1);
+
+    DBUG("Use '%s' as intrinsic conductivity parameter.",
+         intrinsic_permeability.name.c_str());
+
+    // Storage coefficient
+    auto& specific_storage = ParameterLib::findParameter<double>(
+        config,
+        //! \ogs_file_param_special{prj__processes__process__THERMO_HYDRO_MECHANICS__specific_storage}
+        "specific_storage", parameters, 1);
+
+    DBUG("Use '%s' as storage coefficient parameter.",
+         specific_storage.name.c_str());
+
+    // Fluid viscosity
+    auto& fluid_viscosity = ParameterLib::findParameter<double>(
+        config,
+        //! \ogs_file_param_special{prj__processes__process__THERMO_HYDRO_MECHANICS__fluid_viscosity}
+        "fluid_viscosity", parameters, 1);
+    DBUG("Use '%s' as fluid viscosity parameter.",
+         fluid_viscosity.name.c_str());
+
+    // Fluid density
+    auto& fluid_density = ParameterLib::findParameter<double>(
+        config,
+        //! \ogs_file_param_special{prj__processes__process__THERMO_HYDRO_MECHANICS__fluid_density}
+        "fluid_density", parameters, 1);
+    DBUG("Use '%s' as fluid density parameter.", fluid_density.name.c_str());
+
+    // Biot coefficient
+    auto& biot_coefficient = ParameterLib::findParameter<double>(
+        config,
+        //! \ogs_file_param_special{prj__processes__process__THERMO_HYDRO_MECHANICS__biot_coefficient}
+        "biot_coefficient", parameters, 1);
+    DBUG("Use '%s' as Biot coefficient parameter.",
+         biot_coefficient.name.c_str());
+
+    // Porosity
+    auto& porosity = ParameterLib::findParameter<double>(
+        config,
+        //! \ogs_file_param_special{prj__processes__process__THERMO_HYDRO_MECHANICS__porosity}
+        "porosity", parameters, 1);
+    DBUG("Use '%s' as porosity parameter.", porosity.name.c_str());
+
+    // Solid density
+    auto& solid_density = ParameterLib::findParameter<double>(
+        config,
+        //! \ogs_file_param_special{prj__processes__process__THERMO_HYDRO_MECHANICS__solid_density}
+        "solid_density", parameters, 1);
+    DBUG("Use '%s' as solid density parameter.", solid_density.name.c_str());
+
+    // linear thermal expansion coefficient for solid
+    auto const& solid_linear_thermal_expansion_coefficient =
+        ParameterLib::findParameter<double>(
+            config,
+            //! \ogs_file_param_special{prj__processes__process__THERMO_HYDRO_MECHANICS__solid_linear_thermal_expansion_coefficient}
+            "solid_linear_thermal_expansion_coefficient", parameters, 1);
+    DBUG("Use '%s' as solid linear thermal expansion coefficient parameter.",
+         solid_linear_thermal_expansion_coefficient.name.c_str());
+
+    // volumetric thermal expansion coefficient for fluid
+    auto const& fluid_volumetric_thermal_expansion_coefficient =
+        ParameterLib::findParameter<double>(
+            config,
+            //! \ogs_file_param_special{prj__processes__process__THERMO_HYDRO_MECHANICS__fluid_volumetric_thermal_expansion_coefficient}
+            "fluid_volumetric_thermal_expansion_coefficient", parameters, 1);
+    DBUG(
+        "Use '%s' as fluid volumetric thermal expansion coefficient "
+        "parameter.",
+        fluid_volumetric_thermal_expansion_coefficient.name.c_str());
+
+    // specific heat capacity for solid
+    auto& solid_specific_heat_capacity = ParameterLib::findParameter<double>(
+        config,
+        //! \ogs_file_param_special{prj__processes__process__THERMO_HYDRO_MECHANICS__solid_specific_heat_capacity}
+        "solid_specific_heat_capacity", parameters, 1);
+    DBUG("Use '%s' as solid specific heat capacity parameter.",
+         solid_specific_heat_capacity.name.c_str());
+
+    // specific heat capacity for fluid
+    auto& fluid_specific_heat_capacity = ParameterLib::findParameter<double>(
+        config,
+        //! \ogs_file_param_special{prj__processes__process__THERMO_HYDRO_MECHANICS__fluid_specific_heat_capacity}
+        "fluid_specific_heat_capacity", parameters, 1);
+    DBUG("Use '%s' as fluid specific heat capacity parameter.",
+         fluid_specific_heat_capacity.name.c_str());
+
+    // thermal conductivity for solid // currently only considers isotropic
+    auto& solid_thermal_conductivity = ParameterLib::findParameter<double>(
+        config,
+        //! \ogs_file_param_special{prj__processes__process__THERMO_HYDRO_MECHANICS__solid_thermal_conductivity}
+        "solid_thermal_conductivity", parameters, 1);
+    DBUG("Use '%s' as solid thermal conductivity parameter.",
+         solid_thermal_conductivity.name.c_str());
+
+    // thermal conductivity for fluid // currently only considers isotropic
+    auto& fluid_thermal_conductivity = ParameterLib::findParameter<double>(
+        config,
+        //! \ogs_file_param_special{prj__processes__process__THERMO_HYDRO_MECHANICS__fluid_thermal_conductivity}
+        "fluid_thermal_conductivity", parameters, 1);
+    DBUG("Use '%s' as fluid thermal conductivity parameter.",
+         fluid_thermal_conductivity.name.c_str());
+
+    // reference temperature
+    auto& reference_temperature = ParameterLib::findParameter<double>(
+        config,
+        //! \ogs_file_param_special{prj__processes__process__THERMO_HYDRO_MECHANICS__reference_temperature}
+        "reference_temperature", parameters, 1);
+    DBUG("Use '%s' as reference temperature parameter.",
+         reference_temperature.name.c_str());
+
+    // Specific body force
+    Eigen::Matrix<double, DisplacementDim, 1> specific_body_force;
+    {
+        std::vector<double> const b =
+            //! \ogs_file_param{prj__processes__process__THERMO_HYDRO_MECHANICS__specific_body_force}
+            config.getConfigParameter<std::vector<double>>(
+                "specific_body_force");
+        if (b.size() != DisplacementDim)
+        {
+            OGS_FATAL(
+                "The size of the specific body force vector does not match the "
+                "displacement dimension. Vector size is %d, displacement "
+                "dimension is %d",
+                b.size(), DisplacementDim);
+        }
+
+        std::copy_n(b.data(), b.size(), specific_body_force.data());
+    }
+
+    ThermoHydroMechanicsProcessData<DisplacementDim> process_data{
+        materialIDs(mesh),
+        std::move(solid_constitutive_relations),
+        intrinsic_permeability,
+        specific_storage,
+        fluid_viscosity,
+        fluid_density,
+        biot_coefficient,
+        porosity,
+        solid_density,
+        solid_linear_thermal_expansion_coefficient,
+        fluid_volumetric_thermal_expansion_coefficient,
+        solid_specific_heat_capacity,
+        fluid_specific_heat_capacity,
+        solid_thermal_conductivity,
+        fluid_thermal_conductivity,
+        reference_temperature,
+        specific_body_force};
+
+    SecondaryVariableCollection secondary_variables;
+
+    NumLib::NamedFunctionCaller named_function_caller(
+        {"ThermoHydroMechanics_displacement"});
+
+    ProcessLib::createSecondaryVariables(config, secondary_variables,
+                                         named_function_caller);
+
+    return std::make_unique<ThermoHydroMechanicsProcess<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),
+        use_monolithic_scheme);
+}
+
+template std::unique_ptr<Process> createThermoHydroMechanicsProcess<2>(
+    MeshLib::Mesh& mesh,
+    std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+    std::vector<ProcessVariable> const& variables,
+    std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
+    unsigned const integration_order,
+    BaseLib::ConfigTree const& config);
+
+template std::unique_ptr<Process> createThermoHydroMechanicsProcess<3>(
+    MeshLib::Mesh& mesh,
+    std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+    std::vector<ProcessVariable> const& variables,
+    std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
+    unsigned const integration_order,
+    BaseLib::ConfigTree const& config);
+
+}  // namespace ThermoHydroMechanics
+}  // namespace ProcessLib
diff --git a/ProcessLib/ThermoHydroMechanics/CreateThermoHydroMechanicsProcess.h b/ProcessLib/ThermoHydroMechanics/CreateThermoHydroMechanicsProcess.h
new file mode 100644
index 0000000000000000000000000000000000000000..59cb38c1bab52ee1af3bb2d348444222d56e977d
--- /dev/null
+++ b/ProcessLib/ThermoHydroMechanics/CreateThermoHydroMechanicsProcess.h
@@ -0,0 +1,64 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2019, 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 <memory>
+#include <vector>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+namespace MeshLib
+{
+class Mesh;
+}
+namespace ParameterLib
+{
+struct ParameterBase;
+}
+namespace ProcessLib
+{
+class AbstractJacobianAssembler;
+class Process;
+class ProcessVariable;
+}  // namespace ProcessLib
+
+namespace ProcessLib
+{
+namespace ThermoHydroMechanics
+{
+template <int DisplacementDim>
+std::unique_ptr<Process> createThermoHydroMechanicsProcess(
+    MeshLib::Mesh& mesh,
+    std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+    std::vector<ProcessVariable> const& variables,
+    std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
+    unsigned const integration_order,
+    BaseLib::ConfigTree const& config);
+
+extern template std::unique_ptr<Process> createThermoHydroMechanicsProcess<2>(
+    MeshLib::Mesh& mesh,
+    std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+    std::vector<ProcessVariable> const& variables,
+    std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
+    unsigned const integration_order,
+    BaseLib::ConfigTree const& config);
+
+extern template std::unique_ptr<Process> createThermoHydroMechanicsProcess<3>(
+    MeshLib::Mesh& mesh,
+    std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+    std::vector<ProcessVariable> const& variables,
+    std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
+    unsigned const integration_order,
+    BaseLib::ConfigTree const& config);
+}  // namespace ThermoHydroMechanics
+}  // namespace ProcessLib
diff --git a/ProcessLib/ThermoHydroMechanics/IntegrationPointData.h b/ProcessLib/ThermoHydroMechanics/IntegrationPointData.h
new file mode 100644
index 0000000000000000000000000000000000000000..8d13b517d39b0e7262f64a90db20670d3b380d70
--- /dev/null
+++ b/ProcessLib/ThermoHydroMechanics/IntegrationPointData.h
@@ -0,0 +1,128 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2019, 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 <memory>
+
+#include "MaterialLib/SolidModels/LinearElasticIsotropic.h"
+#include "MathLib/KelvinVector.h"
+#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "ParameterLib/Parameter.h"
+
+namespace ProcessLib
+{
+namespace ThermoHydroMechanics
+{
+template <typename BMatricesType, typename ShapeMatrixTypeDisplacement,
+          typename ShapeMatricesTypePressure, int DisplacementDim, int NPoints>
+struct IntegrationPointData final
+{
+    explicit IntegrationPointData(
+        MaterialLib::Solids::MechanicsBase<DisplacementDim> const&
+            solid_material)
+        : solid_material(solid_material),
+          material_state_variables(
+              solid_material.createMaterialStateVariables())
+    {
+        // Initialize current time step values
+        static const int kelvin_vector_size =
+            MathLib::KelvinVector::KelvinVectorDimensions<
+                DisplacementDim>::value;
+        sigma_eff.setZero(kelvin_vector_size);
+        eps.setZero(kelvin_vector_size);
+        eps_m.setZero(kelvin_vector_size);
+        eps_m_prev.resize(kelvin_vector_size);
+
+        // Previous time step values are not initialized and are set later.
+        eps_prev.resize(kelvin_vector_size);
+        sigma_eff_prev.resize(kelvin_vector_size);
+    }
+
+    typename ShapeMatrixTypeDisplacement::template MatrixType<
+        DisplacementDim, NPoints * DisplacementDim>
+        N_u_op;
+    typename BMatricesType::KelvinVectorType sigma_eff, sigma_eff_prev;
+    typename BMatricesType::KelvinVectorType eps, eps_prev;
+    typename BMatricesType::KelvinVectorType eps_m, eps_m_prev;
+
+    typename ShapeMatrixTypeDisplacement::NodalRowVectorType N_u;
+    typename ShapeMatrixTypeDisplacement::GlobalDimNodalMatrixType dNdx_u;
+
+    typename ShapeMatricesTypePressure::NodalRowVectorType N_p;
+    typename ShapeMatricesTypePressure::GlobalDimNodalMatrixType dNdx_p;
+
+    MaterialLib::Solids::MechanicsBase<DisplacementDim> const& solid_material;
+    std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
+        DisplacementDim>::MaterialStateVariables>
+        material_state_variables;
+    double integration_weight;
+
+    void pushBackState()
+    {
+        eps_m_prev = eps_m;
+        sigma_eff_prev = sigma_eff;
+        material_state_variables->pushBackState();
+    }
+
+    template <typename DisplacementVectorType>
+    typename BMatricesType::KelvinMatrixType updateConstitutiveRelation(
+        double const t,
+        ParameterLib::SpatialPosition const& x_position,
+        double const dt,
+        DisplacementVectorType const& /*u*/,
+        double const T)
+    {
+        auto&& solution = solid_material.integrateStress(
+            t, x_position, dt, eps_m_prev, eps_m, sigma_eff_prev,
+            *material_state_variables, T);
+
+        if (!solution)
+            OGS_FATAL("Computation of local constitutive relation failed.");
+
+        MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> C;
+        std::tie(sigma_eff, material_state_variables, C) = std::move(*solution);
+
+        return C;
+    }
+
+    template <typename DisplacementVectorType>
+    typename BMatricesType::KelvinMatrixType updateConstitutiveRelationThermal(
+        double const t,
+        ParameterLib::SpatialPosition const& x_position,
+        double const dt,
+        DisplacementVectorType const& /*u*/,
+        double const T,
+        double const thermal_strain)
+    {
+        auto const& identity2 = MathLib::KelvinVector::Invariants<
+            MathLib::KelvinVector::KelvinVectorDimensions<
+                DisplacementDim>::value>::identity2;
+
+        // assume isotropic thermal expansion
+        eps_m.noalias() = eps - thermal_strain * identity2;
+        auto&& solution = solid_material.integrateStress(
+            t, x_position, dt, eps_m_prev, eps_m, sigma_eff_prev,
+            *material_state_variables, T);
+
+        if (!solution)
+            OGS_FATAL("Computation of local constitutive relation failed.");
+
+        MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> C;
+        std::tie(sigma_eff, material_state_variables, C) = std::move(*solution);
+
+        return C;
+    }
+
+    EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
+};
+
+}  // namespace ThermoHydroMechanics
+}  // namespace ProcessLib
diff --git a/ProcessLib/ThermoHydroMechanics/LocalAssemblerInterface.h b/ProcessLib/ThermoHydroMechanics/LocalAssemblerInterface.h
new file mode 100644
index 0000000000000000000000000000000000000000..c528e5ed678b46d0262bfb8abdd8ea8b1e63e2eb
--- /dev/null
+++ b/ProcessLib/ThermoHydroMechanics/LocalAssemblerInterface.h
@@ -0,0 +1,45 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2019, 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 "NumLib/Extrapolation/ExtrapolatableElement.h"
+#include "ProcessLib/LocalAssemblerInterface.h"
+
+namespace ProcessLib
+{
+namespace ThermoHydroMechanics
+{
+struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
+                                 public NumLib::ExtrapolatableElement
+{
+    virtual std::vector<double> getSigma() const = 0;
+
+    virtual std::vector<double> const& getIntPtSigma(
+        const double /*t*/,
+        GlobalVector const& /*current_solution*/,
+        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+        std::vector<double>& cache) const = 0;
+
+    virtual std::vector<double> const& getIntPtEpsilon(
+        const double /*t*/,
+        GlobalVector const& /*current_solution*/,
+        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+        std::vector<double>& cache) const = 0;
+
+    virtual std::vector<double> const& getIntPtDarcyVelocity(
+        const double /*t*/,
+        GlobalVector const& /*current_solution*/,
+        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+        std::vector<double>& cache) const = 0;
+};
+
+}  // namespace ThermoHydroMechanics
+}  // namespace ProcessLib
diff --git a/ProcessLib/ThermoHydroMechanics/LocalDataInitializer.h b/ProcessLib/ThermoHydroMechanics/LocalDataInitializer.h
new file mode 100644
index 0000000000000000000000000000000000000000..f8801c2b3315c5c69a4f8719fb7aeb81870c6069
--- /dev/null
+++ b/ProcessLib/ThermoHydroMechanics/LocalDataInitializer.h
@@ -0,0 +1,290 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2019, 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 <functional>
+#include <memory>
+#include <type_traits>
+#include <typeindex>
+#include <typeinfo>
+#include <unordered_map>
+
+#include "MeshLib/Elements/Elements.h"
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "NumLib/Fem/FiniteElement/LowerDimShapeTable.h"
+#include "NumLib/Fem/Integration/GaussLegendreIntegrationPolicy.h"
+
+#ifndef OGS_MAX_ELEMENT_DIM
+static_assert(false, "The macro OGS_MAX_ELEMENT_DIM is undefined.");
+#endif
+
+#ifndef OGS_MAX_ELEMENT_ORDER
+static_assert(false, "The macro OGS_MAX_ELEMENT_ORDER is undefined.");
+#endif
+
+// The following macros decide which element types will be compiled, i.e.
+// which element types will be available for use in simulations.
+
+#ifdef OGS_ENABLE_ELEMENT_SIMPLEX
+#define ENABLED_ELEMENT_TYPE_SIMPLEX 1u
+#else
+#define ENABLED_ELEMENT_TYPE_SIMPLEX 0u
+#endif
+
+#ifdef OGS_ENABLE_ELEMENT_CUBOID
+#define ENABLED_ELEMENT_TYPE_CUBOID 1u << 1
+#else
+#define ENABLED_ELEMENT_TYPE_CUBOID 0u
+#endif
+
+#ifdef OGS_ENABLE_ELEMENT_PRISM
+#define ENABLED_ELEMENT_TYPE_PRISM 1u << 2
+#else
+#define ENABLED_ELEMENT_TYPE_PRISM 0u
+#endif
+
+#ifdef OGS_ENABLE_ELEMENT_PYRAMID
+#define ENABLED_ELEMENT_TYPE_PYRAMID 1u << 3
+#else
+#define ENABLED_ELEMENT_TYPE_PYRAMID 0u
+#endif
+
+// Dependent element types.
+// Faces of tets, pyramids and prisms are triangles
+#define ENABLED_ELEMENT_TYPE_TRI                                       \
+    ((ENABLED_ELEMENT_TYPE_SIMPLEX) | (ENABLED_ELEMENT_TYPE_PYRAMID) | \
+     (ENABLED_ELEMENT_TYPE_PRISM))
+// Faces of hexes, pyramids and prisms are quads
+#define ENABLED_ELEMENT_TYPE_QUAD                                     \
+    ((ENABLED_ELEMENT_TYPE_CUBOID) | (ENABLED_ELEMENT_TYPE_PYRAMID) | \
+     (ENABLED_ELEMENT_TYPE_PRISM))
+
+// All enabled element types
+#define OGS_ENABLED_ELEMENTS                                          \
+    ((ENABLED_ELEMENT_TYPE_SIMPLEX) | (ENABLED_ELEMENT_TYPE_CUBOID) | \
+     (ENABLED_ELEMENT_TYPE_PYRAMID) | (ENABLED_ELEMENT_TYPE_PRISM))
+
+// Include only what is needed (Well, the conditions are not sharp).
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_SIMPLEX) != 0
+#include "NumLib/Fem/ShapeFunction/ShapeTet10.h"
+#include "NumLib/Fem/ShapeFunction/ShapeTet4.h"
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_TRI) != 0
+#include "NumLib/Fem/ShapeFunction/ShapeTri3.h"
+#include "NumLib/Fem/ShapeFunction/ShapeTri6.h"
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_CUBOID) != 0
+#include "NumLib/Fem/ShapeFunction/ShapeHex20.h"
+#include "NumLib/Fem/ShapeFunction/ShapeHex8.h"
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_QUAD) != 0
+#include "NumLib/Fem/ShapeFunction/ShapeQuad4.h"
+#include "NumLib/Fem/ShapeFunction/ShapeQuad8.h"
+#include "NumLib/Fem/ShapeFunction/ShapeQuad9.h"
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PRISM) != 0
+#include "NumLib/Fem/ShapeFunction/ShapePrism15.h"
+#include "NumLib/Fem/ShapeFunction/ShapePrism6.h"
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PYRAMID) != 0
+#include "NumLib/Fem/ShapeFunction/ShapePyra13.h"
+#include "NumLib/Fem/ShapeFunction/ShapePyra5.h"
+#endif
+
+namespace ProcessLib
+{
+namespace ThermoHydroMechanics
+{
+/// The LocalDataInitializer is a functor creating a local assembler data with
+/// corresponding to the mesh element type shape functions and calling
+/// initialization of the new local assembler data.
+/// For example for MeshLib::Quad a local assembler data with template argument
+/// NumLib::ShapeQuad4 is created.
+///
+/// \attention This is modified version of the ProcessLib::LocalDataInitializer
+/// class which does not include line elements, allows only shapefunction of
+/// order 2.
+template <typename LocalAssemblerInterface,
+          template <typename, typename, typename, int>
+          class ThermoHydroMechanicsLocalAssembler,
+          int GlobalDim, typename... ConstructorArgs>
+class LocalDataInitializer final
+{
+public:
+    using LADataIntfPtr = std::unique_ptr<LocalAssemblerInterface>;
+
+    LocalDataInitializer(NumLib::LocalToGlobalIndexMap const& dof_table,
+                         const unsigned shapefunction_order)
+        : _dof_table(dof_table)
+    {
+        if (shapefunction_order != 2)
+            OGS_FATAL(
+                "The given shape function order %d is not supported.\nOnly "
+                "shape functions of order 2 are supported.",
+                shapefunction_order);
+            // /// Quads and Hexahedra ///////////////////////////////////
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_QUAD) != 0 && \
+    OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 2
+        _builder[std::type_index(typeid(MeshLib::Quad8))] =
+            makeLocalAssemblerBuilder<NumLib::ShapeQuad8>();
+        _builder[std::type_index(typeid(MeshLib::Quad9))] =
+            makeLocalAssemblerBuilder<NumLib::ShapeQuad9>();
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_CUBOID) != 0 && \
+    OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
+        _builder[std::type_index(typeid(MeshLib::Hex20))] =
+            makeLocalAssemblerBuilder<NumLib::ShapeHex20>();
+#endif
+
+        // /// Simplices ////////////////////////////////////////////////
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_TRI) != 0 && \
+    OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 2
+        _builder[std::type_index(typeid(MeshLib::Tri6))] =
+            makeLocalAssemblerBuilder<NumLib::ShapeTri6>();
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_SIMPLEX) != 0 && \
+    OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
+        _builder[std::type_index(typeid(MeshLib::Tet10))] =
+            makeLocalAssemblerBuilder<NumLib::ShapeTet10>();
+#endif
+
+        // /// Prisms ////////////////////////////////////////////////////
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PRISM) != 0 && \
+    OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
+        _builder[std::type_index(typeid(MeshLib::Prism15))] =
+            makeLocalAssemblerBuilder<NumLib::ShapePrism15>();
+#endif
+
+        // /// Pyramids //////////////////////////////////////////////////
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PYRAMID) != 0 && \
+    OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
+        _builder[std::type_index(typeid(MeshLib::Pyramid13))] =
+            makeLocalAssemblerBuilder<NumLib::ShapePyra13>();
+#endif
+    }
+
+    /// Sets the provided \c data_ptr to the newly created local assembler data.
+    ///
+    /// \attention
+    /// The index \c id is not necessarily the mesh item's id. Especially when
+    /// having multiple meshes it will differ from the latter.
+    void operator()(std::size_t const id,
+                    MeshLib::Element const& mesh_item,
+                    LADataIntfPtr& data_ptr,
+                    ConstructorArgs&&... args) const
+    {
+        auto const type_idx = std::type_index(typeid(mesh_item));
+        auto const it = _builder.find(type_idx);
+
+        if (it != _builder.end())
+        {
+            auto const num_local_dof = _dof_table.getNumberOfElementDOF(id);
+            data_ptr = it->second(mesh_item, num_local_dof,
+                                  std::forward<ConstructorArgs>(args)...);
+        }
+        else
+        {
+            OGS_FATAL(
+                "You are trying to build a local assembler for an unknown mesh "
+                "element type (%s)."
+                " Maybe you have disabled this mesh element type in your build "
+                "configuration or this process requires higher order elements.",
+                type_idx.name());
+        }
+    }
+
+private:
+    using LADataBuilder =
+        std::function<LADataIntfPtr(MeshLib::Element const& e,
+                                    std::size_t const local_matrix_size,
+                                    ConstructorArgs&&...)>;
+
+    template <typename ShapeFunctionDisplacement>
+    using IntegrationMethod = typename NumLib::GaussLegendreIntegrationPolicy<
+        typename ShapeFunctionDisplacement::MeshElement>::IntegrationMethod;
+
+    template <typename ShapeFunctionDisplacement,
+              typename ShapeFunctionPressure>
+    using LAData = ThermoHydroMechanicsLocalAssembler<
+        ShapeFunctionDisplacement, ShapeFunctionPressure,
+        IntegrationMethod<ShapeFunctionDisplacement>, GlobalDim>;
+
+    /// A helper forwarding to the correct version of makeLocalAssemblerBuilder
+    /// depending whether the global dimension is less than the shape function's
+    /// dimension or not.
+    template <typename ShapeFunctionDisplacement>
+    static LADataBuilder makeLocalAssemblerBuilder()
+    {
+        return makeLocalAssemblerBuilder<ShapeFunctionDisplacement>(
+            static_cast<std::integral_constant<
+                bool, (GlobalDim >= ShapeFunctionDisplacement::DIM)>*>(
+                nullptr));
+    }
+
+    /// Mapping of element types to local assembler constructors.
+    std::unordered_map<std::type_index, LADataBuilder> _builder;
+
+    NumLib::LocalToGlobalIndexMap const& _dof_table;
+
+    // local assembler builder implementations.
+private:
+    /// Generates a function that creates a new LocalAssembler of type
+    /// LAData<ShapeFunctionDisplacement>. Only functions with shape function's
+    /// dimension less or equal to the global dimension are instantiated, e.g.
+    /// following combinations of shape functions and global dimensions: (Line2,
+    /// 1),
+    /// (Line2, 2), (Line2, 3), (Hex20, 3) but not (Hex20, 2) or (Hex20, 1).
+    template <typename ShapeFunctionDisplacement>
+    static LADataBuilder makeLocalAssemblerBuilder(std::true_type*)
+    {
+        // (Lower order elements = Order(ShapeFunctionDisplacement) - 1).
+        using ShapeFunctionPressure =
+            typename NumLib::LowerDim<ShapeFunctionDisplacement>::type;
+        return [](MeshLib::Element const& e,
+                  std::size_t const local_matrix_size,
+                  ConstructorArgs&&... args) {
+            return LADataIntfPtr{
+                new LAData<ShapeFunctionDisplacement, ShapeFunctionPressure>{
+                    e, local_matrix_size,
+                    std::forward<ConstructorArgs>(args)...}};
+        };
+    }
+
+    /// Returns nullptr for shape functions whose dimensions are less than the
+    /// global dimension.
+    template <typename ShapeFunctionDisplacement>
+    static LADataBuilder makeLocalAssemblerBuilder(std::false_type*)
+    {
+        return nullptr;
+    }
+};
+
+}  // namespace ThermoHydroMechanics
+}  // namespace ProcessLib
+
+#undef ENABLED_ELEMENT_TYPE_SIMPLEX
+#undef ENABLED_ELEMENT_TYPE_CUBOID
+#undef ENABLED_ELEMENT_TYPE_PYRAMID
+#undef ENABLED_ELEMENT_TYPE_PRISM
+#undef ENABLED_ELEMENT_TYPE_TRI
+#undef ENABLED_ELEMENT_TYPE_QUAD
+#undef OGS_ENABLED_ELEMENTS
diff --git a/ProcessLib/ThermoHydroMechanics/Tests.cmake b/ProcessLib/ThermoHydroMechanics/Tests.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..d68df291ac1964d7acca4f37cc278f6a9c04911a
--- /dev/null
+++ b/ProcessLib/ThermoHydroMechanics/Tests.cmake
@@ -0,0 +1,70 @@
+# ThermoHydroMechanics; Small deformation, linear poroelastic, homogeneous
+AddTest(
+    NAME ThermoHydroMechanics_square_1e0
+    PATH ThermoHydroMechanics/Linear/Square_sealed_homogeneous
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS square_1e0.prj
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT OGS_USE_MPI
+    DIFF_DATA
+    expected_square_1e0_pcs_0_ts_10_t_1000.000000.vtu square_1e0_pcs_0_ts_10_t_1000.000000.vtu displacement displacement 1e-8 1e-8
+    expected_square_1e0_pcs_0_ts_10_t_1000.000000.vtu square_1e0_pcs_0_ts_10_t_1000.000000.vtu pressure pressure 1e-8 1e-8
+    expected_square_1e0_pcs_0_ts_10_t_1000.000000.vtu square_1e0_pcs_0_ts_10_t_1000.000000.vtu temperature temperature 1e-8 1e-8
+    expected_square_1e0_pcs_0_ts_10_t_1000.000000.vtu square_1e0_pcs_0_ts_10_t_1000.000000.vtu epsilon epsilon 1e-8 1e-8
+    expected_square_1e0_pcs_0_ts_10_t_1000.000000.vtu square_1e0_pcs_0_ts_10_t_1000.000000.vtu sigma sigma 1e-8 1e-8
+)
+
+# ThermoHydroMechanics; Small deformation, linear poroelastic, sealed, bimaterial
+AddTest(
+    NAME ThermoHydroMechanics_square_1e2_sealed_bimaterial
+    PATH ThermoHydroMechanics/Linear/Beam_sealed_bimaterial
+    RUNTIME 5
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS square_1e2.prj
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT OGS_USE_MPI
+    DIFF_DATA
+    expected_square_1e2_pcs_0_ts_10_t_100.000000.vtu square_1e2_pcs_0_ts_10_t_100.000000.vtu displacement displacement 1e-8 1e-8
+    expected_square_1e2_pcs_0_ts_10_t_100.000000.vtu square_1e2_pcs_0_ts_10_t_100.000000.vtu pressure pressure 1e-8 1e-8
+    expected_square_1e2_pcs_0_ts_10_t_100.000000.vtu square_1e2_pcs_0_ts_10_t_100.000000.vtu temperature temperature 1e-8 1e-8
+    expected_square_1e2_pcs_0_ts_10_t_100.000000.vtu square_1e2_pcs_0_ts_10_t_100.000000.vtu epsilon epsilon 1e-8 1e-8
+    expected_square_1e2_pcs_0_ts_10_t_100.000000.vtu square_1e2_pcs_0_ts_10_t_100.000000.vtu sigma sigma 1e-8 1e-8
+)
+
+# ThermoHydroMechanics; Small deformation, linear poroelastic, unsealed, bimaterial
+AddTest(
+    NAME ThermoHydroMechanics_square_1e2_unsealed_bimaterial
+    PATH ThermoHydroMechanics/Linear/Beam_unsealed_bimaterial
+    RUNTIME 5
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS square_1e2.prj
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT OGS_USE_MPI
+    DIFF_DATA
+    expected_square_1e2_pcs_0_ts_10_t_1000.000000.vtu square_1e2_pcs_0_ts_10_t_1000.000000.vtu displacement displacement 1e-8 1e-8
+    expected_square_1e2_pcs_0_ts_10_t_1000.000000.vtu square_1e2_pcs_0_ts_10_t_1000.000000.vtu pressure pressure 1e-8 1e-8
+    expected_square_1e2_pcs_0_ts_10_t_1000.000000.vtu square_1e2_pcs_0_ts_10_t_1000.000000.vtu temperature temperature 1e-8 1e-8
+    expected_square_1e2_pcs_0_ts_10_t_1000.000000.vtu square_1e2_pcs_0_ts_10_t_1000.000000.vtu epsilon epsilon 1e-8 1e-8
+    expected_square_1e2_pcs_0_ts_10_t_1000.000000.vtu square_1e2_pcs_0_ts_10_t_1000.000000.vtu sigma sigma 1e-8 1e-8
+)
+
+# ThermoHydroMechanics; Small deformation, linear poroelastic, point heat source consolidation
+AddTest(
+    NAME LARGE_ThermoHydroMechanics_square_1e2_point_heat_injection
+    PATH ThermoHydroMechanics/Linear/Point_injection
+    RUNTIME 400
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS square_1e2.prj
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT OGS_USE_MPI
+    DIFF_DATA
+    expected_square_1e0_pcs_0_ts_10_t_50000.000000.vtu square_1e0_pcs_0_ts_10_t_50000.000000.vtu displacement displacement 1e-5 1e-5
+    expected_square_1e0_pcs_0_ts_10_t_50000.000000.vtu square_1e0_pcs_0_ts_10_t_50000.000000.vtu pressure pressure 1e-5 1e-5
+    expected_square_1e0_pcs_0_ts_10_t_50000.000000.vtu square_1e0_pcs_0_ts_10_t_50000.000000.vtu temperature temperature 1e-5 1e-5
+    expected_square_1e0_pcs_0_ts_10_t_50000.000000.vtu square_1e0_pcs_0_ts_10_t_50000.000000.vtu epsilon epsilon 1e-5 1e-5
+    expected_square_1e0_pcs_0_ts_10_t_50000.000000.vtu square_1e0_pcs_0_ts_10_t_50000.000000.vtu sigma sigma 1e-5 1e-5
+)
diff --git a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..aab59804495927c48351932e1a50b026ee35a6c8
--- /dev/null
+++ b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h
@@ -0,0 +1,525 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2019, 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 "ThermoHydroMechanicsFEM.h"
+
+#include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
+#include "MathLib/KelvinVector.h"
+#include "NumLib/Function/Interpolation.h"
+#include "ProcessLib/CoupledSolutionsForStaggeredScheme.h"
+
+namespace ProcessLib
+{
+namespace ThermoHydroMechanics
+{
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+          typename IntegrationMethod, int DisplacementDim>
+ThermoHydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
+                                   ShapeFunctionPressure, IntegrationMethod,
+                                   DisplacementDim>::
+    ThermoHydroMechanicsLocalAssembler(
+        MeshLib::Element const& e,
+        std::size_t const /*local_matrix_size*/,
+        bool const is_axially_symmetric,
+        unsigned const integration_order,
+        ThermoHydroMechanicsProcessData<DisplacementDim>& process_data)
+    : _process_data(process_data),
+      _integration_method(integration_order),
+      _element(e),
+      _is_axially_symmetric(is_axially_symmetric)
+{
+    unsigned const n_integration_points =
+        _integration_method.getNumberOfPoints();
+
+    _ip_data.reserve(n_integration_points);
+    _secondary_data.N_u.resize(n_integration_points);
+
+    auto const shape_matrices_u =
+        initShapeMatrices<ShapeFunctionDisplacement,
+                          ShapeMatricesTypeDisplacement, IntegrationMethod,
+                          DisplacementDim>(e, is_axially_symmetric,
+                                           _integration_method);
+
+    auto const shape_matrices_p =
+        initShapeMatrices<ShapeFunctionPressure, ShapeMatricesTypePressure,
+                          IntegrationMethod, DisplacementDim>(
+            e, is_axially_symmetric, _integration_method);
+
+    auto const& solid_material =
+        MaterialLib::Solids::selectSolidConstitutiveRelation(
+            _process_data.solid_materials, _process_data.material_ids,
+            e.getID());
+
+    for (unsigned ip = 0; ip < n_integration_points; ip++)
+    {
+        _ip_data.emplace_back(solid_material);
+        auto& ip_data = _ip_data[ip];
+        auto const& sm_u = shape_matrices_u[ip];
+        ip_data.integration_weight =
+            _integration_method.getWeightedPoint(ip).getWeight() *
+            sm_u.integralMeasure * sm_u.detJ;
+
+        ip_data.N_u_op = ShapeMatricesTypeDisplacement::template MatrixType<
+            DisplacementDim, displacement_size>::Zero(DisplacementDim,
+                                                      displacement_size);
+        for (int i = 0; i < DisplacementDim; ++i)
+            ip_data.N_u_op
+                .template block<1, displacement_size / DisplacementDim>(
+                    i, i * displacement_size / DisplacementDim)
+                .noalias() = sm_u.N;
+
+        ip_data.N_u = sm_u.N;
+        ip_data.dNdx_u = sm_u.dNdx;
+
+        ip_data.N_p = shape_matrices_p[ip].N;
+        ip_data.dNdx_p = shape_matrices_p[ip].dNdx;
+
+        _secondary_data.N_u[ip] = shape_matrices_u[ip].N;
+    }
+}
+
+// Assembles the local Jacobian matrix. So far, the linearisation of HT part is
+// not considered as that in HT process.
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+          typename IntegrationMethod, int DisplacementDim>
+void ThermoHydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
+                                        ShapeFunctionPressure,
+                                        IntegrationMethod, DisplacementDim>::
+    assembleWithJacobian(double const t, std::vector<double> const& local_x,
+                         std::vector<double> const& local_xdot,
+                         const double /*dxdot_dx*/, const double /*dx_dx*/,
+                         std::vector<double>& /*local_M_data*/,
+                         std::vector<double>& /*local_K_data*/,
+                         std::vector<double>& local_rhs_data,
+                         std::vector<double>& local_Jac_data)
+{
+    assert(local_x.size() ==
+           pressure_size + displacement_size + temperature_size);
+
+    auto T = Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
+        temperature_size> const>(local_x.data() + temperature_index,
+                                 temperature_size);
+
+    auto p = Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
+        pressure_size> const>(local_x.data() + pressure_index, pressure_size);
+
+    auto u =
+        Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
+            displacement_size> const>(local_x.data() + displacement_index,
+                                      displacement_size);
+
+    auto T_dot =
+        Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
+            temperature_size> const>(local_xdot.data() + temperature_index,
+                                     temperature_size);
+
+    auto p_dot =
+        Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
+            pressure_size> const>(local_xdot.data() + pressure_index,
+                                  pressure_size);
+    auto u_dot =
+        Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
+            displacement_size> const>(local_xdot.data() + displacement_index,
+                                      displacement_size);
+
+    auto local_Jac = MathLib::createZeroedMatrix<
+        typename ShapeMatricesTypeDisplacement::template MatrixType<
+            temperature_size + displacement_size + pressure_size,
+            temperature_size + displacement_size + pressure_size>>(
+        local_Jac_data, displacement_size + pressure_size + temperature_size,
+        displacement_size + pressure_size + temperature_size);
+
+    auto local_rhs = MathLib::createZeroedVector<
+        typename ShapeMatricesTypeDisplacement::template VectorType<
+            displacement_size + pressure_size + temperature_size>>(
+        local_rhs_data, displacement_size + pressure_size + temperature_size);
+
+    typename ShapeMatricesTypePressure::NodalMatrixType MTT;
+    MTT.setZero(temperature_size, temperature_size);
+
+    typename ShapeMatricesTypePressure::NodalMatrixType KTT_coeff;
+    KTT_coeff.setZero(temperature_size, temperature_size);
+
+    typename ShapeMatricesTypePressure::NodalMatrixType KTT;
+    KTT.setZero(temperature_size, temperature_size);
+
+    typename ShapeMatricesTypePressure::NodalMatrixType KTp;
+    KTp.setZero(temperature_size, pressure_size);
+
+    typename ShapeMatricesTypePressure::NodalMatrixType KTp_coeff;
+    KTp_coeff.setZero(temperature_size, pressure_size);
+
+    typename ShapeMatricesTypePressure::NodalMatrixType laplace_p;
+    laplace_p.setZero(pressure_size, pressure_size);
+
+    typename ShapeMatricesTypePressure::NodalMatrixType storage_p;
+    storage_p.setZero(pressure_size, pressure_size);
+
+    typename ShapeMatricesTypePressure::NodalMatrixType storage_T;
+    storage_T.setZero(pressure_size, temperature_size);
+
+    typename ShapeMatricesTypeDisplacement::template MatrixType<
+        displacement_size, pressure_size>
+        Kup;
+    Kup.setZero(displacement_size, pressure_size);
+    typename ShapeMatricesTypeDisplacement::template MatrixType<
+        displacement_size, temperature_size>
+        KuT;
+    KuT.setZero(displacement_size, temperature_size);
+
+    double const& dt = _process_data.dt;
+
+    ParameterLib::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& N_u_op = _ip_data[ip].N_u_op;
+
+        auto const& N_u = _ip_data[ip].N_u;
+        auto const& dNdx_u = _ip_data[ip].dNdx_u;
+
+        auto const& N_p = _ip_data[ip].N_p;
+        auto const& dNdx_p = _ip_data[ip].dNdx_p;
+
+        // same shape function for pressure and temperature since they have the
+        // same order
+        auto const& N_T = N_p;
+        auto const& dNdx_T = dNdx_p;
+        auto const T_int_pt = N_T * T;
+
+        auto const x_coord =
+            interpolateXCoordinate<ShapeFunctionDisplacement,
+                                   ShapeMatricesTypeDisplacement>(_element,
+                                                                  N_u);
+        auto const B =
+            LinearBMatrix::computeBMatrix<DisplacementDim,
+                                          ShapeFunctionDisplacement::NPOINTS,
+                                          typename BMatricesType::BMatrixType>(
+                dNdx_u, N_u, x_coord, _is_axially_symmetric);
+
+        auto& eps = _ip_data[ip].eps;
+        auto const& sigma_eff = _ip_data[ip].sigma_eff;
+
+        double const S = _process_data.specific_storage(t, x_position)[0];
+        // TODO: change it like auto const K_over_mu
+        // hydraulicConductivity<GlobalDim>(_process_data.intrinsic_permeability(t,
+        // x_position)/...
+        double const K_over_mu =
+            _process_data.intrinsic_permeability(t, x_position)[0] /
+            _process_data.fluid_viscosity(t, x_position)[0];
+        double const alpha_s =
+            _process_data.solid_linear_thermal_expansion_coefficient(
+                t, x_position)[0];
+        double const beta_f =
+            _process_data.fluid_volumetric_thermal_expansion_coefficient(
+                t, x_position)[0];
+        // TODO: change it like auto const lambda_f =
+        // hydraulicConductivity<GlobalDim>(_process_data.intrinsic_permeability(t,
+        // x_position)/...
+        double const lambda_f =
+            _process_data.fluid_thermal_conductivity(t, x_position)[0];
+        // TODO: change it like auto const lambda_s =
+        // hydraulicConductivity<GlobalDim>(_process_data.intrinsic_permeability(t,
+        // x_position)/...
+        double const lambda_s =
+            _process_data.solid_thermal_conductivity(t, x_position)[0];
+        double const C_f =
+            _process_data.fluid_specific_heat_capacity(t, x_position)[0];
+        double const C_s =
+            _process_data.solid_specific_heat_capacity(t, x_position)[0];
+        double const T0 = _process_data.reference_temperature(t, x_position)[0];
+        auto const alpha = _process_data.biot_coefficient(t, x_position)[0];
+        auto const rho_sr = _process_data.solid_density(t, x_position)[0];
+        auto const rho_fr = _process_data.fluid_density(t, x_position)[0];
+        auto const porosity = _process_data.porosity(t, x_position)[0];
+        auto const& b = _process_data.specific_body_force;
+        auto const& identity2 = MathLib::KelvinVector::Invariants<
+            MathLib::KelvinVector::KelvinVectorDimensions<
+                DisplacementDim>::value>::identity2;
+
+        double const delta_T(T_int_pt - T0);
+        double const thermal_strain = alpha_s * delta_T;
+
+        double const rho_s = rho_sr * (1 - 3 * thermal_strain);
+
+        auto velocity = (-K_over_mu * dNdx_p * p).eval();
+        double const rho_f = rho_fr * (1 - beta_f * delta_T);
+        velocity += K_over_mu * rho_f * b;
+
+        //
+        // displacement equation, displacement part
+        //
+        eps.noalias() = B * u;
+        auto C = _ip_data[ip].updateConstitutiveRelationThermal(
+            t, x_position, dt, u,
+            _process_data.reference_temperature(t, x_position)[0],
+            thermal_strain);
+
+        local_Jac
+            .template block<displacement_size, displacement_size>(
+                displacement_index, displacement_index)
+            .noalias() += B.transpose() * C * B * w;
+
+        auto const rho = rho_s * (1 - porosity) + porosity * rho_f;
+        local_rhs.template segment<displacement_size>(displacement_index)
+            .noalias() -=
+            (B.transpose() * sigma_eff - N_u_op.transpose() * rho * b) * w;
+
+        //
+        // displacement equation, pressure part (K_up)
+        //
+        Kup.noalias() += B.transpose() * alpha * identity2 * N_p * w;
+
+        //
+        // pressure equation, pressure part (K_pp and M_pp).
+        //
+        laplace_p.noalias() += dNdx_p.transpose() * K_over_mu * dNdx_p * w;
+
+        storage_p.noalias() += N_p.transpose() * S * N_p * w;
+        //
+        //  RHS, pressure part
+        //
+        local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
+            dNdx_p.transpose() * rho_f * K_over_mu * b * w;
+        //
+        // pressure equation, temperature part (M_pT)
+        //
+        auto const beta = porosity * beta_f + (1 - porosity) * 3 * alpha_s;
+        storage_T.noalias() += N_T.transpose() * beta * N_T * w;
+
+        //
+        // pressure equation, displacement part.
+        //
+        // Reusing Kup.transpose().
+
+        //
+        // temperature equation, temperature part.
+        //
+        auto const lambda = porosity * lambda_f + (1 - porosity) * lambda_s;
+        KTT.noalias() += (dNdx_T.transpose() * lambda * dNdx_T +
+                          dNdx_T.transpose() * velocity * N_p * rho_f * C_f) *
+                         w;
+        // coefficient matrix which is used for calculating the residual
+        auto const heat_capacity =
+            porosity * C_f * rho_f + (1 - porosity) * C_s * rho_sr;
+        MTT.noalias() += N_T.transpose() * heat_capacity * N_T * w;
+
+        //
+        // temperature equation, pressure part
+        //
+        KTp.noalias() += K_over_mu * rho_f * C_f * N_T.transpose() *
+                         (dNdx_T * T).transpose() * dNdx_T * w;
+    }
+    // temperature equation, temperature part
+    local_Jac
+        .template block<temperature_size, temperature_size>(temperature_index,
+                                                            temperature_index)
+        .noalias() += KTT + MTT / dt;
+
+    // temperature equation, pressure part
+    local_Jac
+        .template block<temperature_size, pressure_size>(temperature_index,
+                                                         pressure_index)
+        .noalias() -= KTp;
+    // displacement equation, temperature part
+    local_Jac
+        .template block<displacement_size, temperature_size>(displacement_index,
+                                                             temperature_index)
+        .noalias() -= KuT;
+
+    // displacement equation, pressure part
+    local_Jac
+        .template block<displacement_size, pressure_size>(displacement_index,
+                                                          pressure_index)
+        .noalias() -= Kup;
+
+    // pressure equation, temperature part.
+    local_Jac
+        .template block<pressure_size, temperature_size>(pressure_index,
+                                                         temperature_index)
+        .noalias() -= storage_T / dt;
+
+    // pressure equation, pressure part.
+    local_Jac
+        .template block<pressure_size, pressure_size>(pressure_index,
+                                                      pressure_index)
+        .noalias() += laplace_p + storage_p / dt;
+
+    // pressure equation, displacement part.
+    local_Jac
+        .template block<pressure_size, displacement_size>(pressure_index,
+                                                          displacement_index)
+        .noalias() += Kup.transpose() / dt;
+
+    // pressure equation (f_p)
+    local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
+        laplace_p * p + storage_p * p_dot - storage_T * T_dot +
+        Kup.transpose() * u_dot;
+
+    // displacement equation (f_u)
+    local_rhs.template segment<displacement_size>(displacement_index)
+        .noalias() += Kup * p;
+
+    // temperature equation (f_T)
+    local_rhs.template segment<temperature_size>(temperature_index).noalias() -=
+        KTT * T + MTT * T_dot;
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+          typename IntegrationMethod, int DisplacementDim>
+std::vector<double> const& ThermoHydroMechanicsLocalAssembler<
+    ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
+    DisplacementDim>::getIntPtDarcyVelocity(const double t,
+                                            GlobalVector const&
+                                                current_solution,
+                                            NumLib::LocalToGlobalIndexMap const&
+                                                dof_table,
+                                            std::vector<double>& cache) const
+{
+    auto const num_intpts = _ip_data.size();
+
+    auto const indices = NumLib::getIndices(_element.getID(), dof_table);
+    assert(!indices.empty());
+    auto const local_x = current_solution.get(indices);
+
+    cache.clear();
+    auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
+        double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
+        cache, DisplacementDim, num_intpts);
+
+    ParameterLib::SpatialPosition pos;
+    pos.setElementID(_element.getID());
+
+    auto p = Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
+        pressure_size> const>(local_x.data() + pressure_index, pressure_size);
+
+    unsigned const n_integration_points =
+        _integration_method.getNumberOfPoints();
+
+    ParameterLib::SpatialPosition x_position;
+    x_position.setElementID(_element.getID());
+    for (unsigned ip = 0; ip < n_integration_points; ip++)
+    {
+        x_position.setIntegrationPoint(ip);
+        double const K_over_mu =
+            _process_data.intrinsic_permeability(t, x_position)[0] /
+            _process_data.fluid_viscosity(t, x_position)[0];
+
+        auto const rho_fr = _process_data.fluid_density(t, x_position)[0];
+        auto const& b = _process_data.specific_body_force;
+
+        // Compute the velocity
+        auto const& dNdx_p = _ip_data[ip].dNdx_p;
+        cache_matrix.col(ip).noalias() =
+            -K_over_mu * dNdx_p * p - K_over_mu * rho_fr * b;
+    }
+
+    return cache;
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+          typename IntegrationMethod, int DisplacementDim>
+void ThermoHydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
+                                        ShapeFunctionPressure,
+                                        IntegrationMethod, DisplacementDim>::
+    postNonLinearSolverConcrete(std::vector<double> const& local_x,
+                                double const t,
+                                bool const use_monolithic_scheme)
+{
+    const int displacement_offset =
+        use_monolithic_scheme ? displacement_index : 0;
+
+    auto u =
+        Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
+            displacement_size> const>(local_x.data() + displacement_offset,
+                                      displacement_size);
+
+    auto T = Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
+        temperature_size> const>(local_x.data() + temperature_index,
+                                 temperature_size);
+
+    double const& dt = _process_data.dt;
+    ParameterLib::SpatialPosition x_position;
+    x_position.setElementID(_element.getID());
+
+    int const n_integration_points = _integration_method.getNumberOfPoints();
+    for (int ip = 0; ip < n_integration_points; ip++)
+    {
+        x_position.setIntegrationPoint(ip);
+        auto const& N_u = _ip_data[ip].N_u;
+        auto const& N_T = _ip_data[ip].N_p;
+        auto const& dNdx_u = _ip_data[ip].dNdx_u;
+
+        auto const x_coord =
+            interpolateXCoordinate<ShapeFunctionDisplacement,
+                                   ShapeMatricesTypeDisplacement>(_element,
+                                                                  N_u);
+        auto const B =
+            LinearBMatrix::computeBMatrix<DisplacementDim,
+                                          ShapeFunctionDisplacement::NPOINTS,
+                                          typename BMatricesType::BMatrixType>(
+                dNdx_u, N_u, x_coord, _is_axially_symmetric);
+
+        double const T0 = _process_data.reference_temperature(t, x_position)[0];
+        double const alpha_s =
+            _process_data.solid_linear_thermal_expansion_coefficient(
+                t, x_position)[0];
+
+        double const T_int_pt = N_T * T;
+
+        double const delta_T(T_int_pt - T0);
+        double const thermal_strain = alpha_s * delta_T;
+
+        auto& eps = _ip_data[ip].eps;
+        eps.noalias() = B * u;
+
+        _ip_data[ip].updateConstitutiveRelationThermal(
+            t, x_position, dt, u,
+            _process_data.reference_temperature(t, x_position)[0],
+            thermal_strain);
+    }
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+          typename IntegrationMethod, int DisplacementDim>
+void ThermoHydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
+                                        ShapeFunctionPressure,
+                                        IntegrationMethod, DisplacementDim>::
+    computeSecondaryVariableConcrete(double const /*t*/,
+                                     std::vector<double> const& local_x)
+{
+    auto p = Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
+        pressure_size> const>(local_x.data() + pressure_index, pressure_size);
+
+    NumLib::interpolateToHigherOrderNodes<
+        ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
+        DisplacementDim>(_element, _is_axially_symmetric, p,
+                         *_process_data.pressure_interpolated);
+
+    auto T = Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
+        pressure_size> const>(local_x.data() + temperature_index,
+                              temperature_size);
+
+    NumLib::interpolateToHigherOrderNodes<
+        ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
+        DisplacementDim>(_element, _is_axially_symmetric, T,
+                         *_process_data.temperature_interpolated);
+}
+
+}  // namespace ThermoHydroMechanics
+}  // namespace ProcessLib
diff --git a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM.h b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM.h
new file mode 100644
index 0000000000000000000000000000000000000000..dc03604ec6a4f69c9ae03775851791a3f8f3cfbf
--- /dev/null
+++ b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM.h
@@ -0,0 +1,262 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2019, 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 <memory>
+#include <vector>
+
+#include "MaterialLib/PhysicalConstant.h"
+#include "MaterialLib/SolidModels/LinearElasticIsotropic.h"
+#include "MathLib/KelvinVector.h"
+#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "NumLib/DOF/DOFTableUtil.h"
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
+#include "ParameterLib/Parameter.h"
+#include "ProcessLib/Deformation/BMatrixPolicy.h"
+#include "ProcessLib/Deformation/LinearBMatrix.h"
+#include "ProcessLib/LocalAssemblerTraits.h"
+#include "ProcessLib/Utils/InitShapeMatrices.h"
+
+#include "IntegrationPointData.h"
+#include "LocalAssemblerInterface.h"
+#include "ThermoHydroMechanicsProcessData.h"
+
+namespace ProcessLib
+{
+namespace ThermoHydroMechanics
+{
+/// Used by for extrapolation of the integration point values. It is ordered
+/// (and stored) by integration points.
+template <typename ShapeMatrixType>
+struct SecondaryData
+{
+    std::vector<ShapeMatrixType, Eigen::aligned_allocator<ShapeMatrixType>> N_u;
+};
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+          typename IntegrationMethod, int DisplacementDim>
+class ThermoHydroMechanicsLocalAssembler : public LocalAssemblerInterface
+{
+public:
+    using ShapeMatricesTypeDisplacement =
+        ShapeMatrixPolicyType<ShapeFunctionDisplacement, DisplacementDim>;
+
+    // Types for pressure.
+    using ShapeMatricesTypePressure =
+        ShapeMatrixPolicyType<ShapeFunctionPressure, DisplacementDim>;
+
+    static int const KelvinVectorSize =
+        MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
+    using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
+
+    ThermoHydroMechanicsLocalAssembler(
+        ThermoHydroMechanicsLocalAssembler const&) = delete;
+    ThermoHydroMechanicsLocalAssembler(ThermoHydroMechanicsLocalAssembler&&) =
+        delete;
+
+    ThermoHydroMechanicsLocalAssembler(
+        MeshLib::Element const& e,
+        std::size_t const /*local_matrix_size*/,
+        bool const is_axially_symmetric,
+        unsigned const integration_order,
+        ThermoHydroMechanicsProcessData<DisplacementDim>& process_data);
+
+    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_rhs_data*/) override
+    {
+        OGS_FATAL(
+            "ThermoHydroMechanicsLocalAssembler: assembly without Jacobian is "
+            "not implemented.");
+    }
+
+    void assembleWithJacobian(double const t,
+                              std::vector<double> const& local_x,
+                              std::vector<double> const& local_xdot,
+                              const double /*dxdot_dx*/, const double /*dx_dx*/,
+                              std::vector<double>& /*local_M_data*/,
+                              std::vector<double>& /*local_K_data*/,
+                              std::vector<double>& local_rhs_data,
+                              std::vector<double>& local_Jac_data) override;
+
+    void preTimestepConcrete(std::vector<double> const& /*local_x*/,
+                             double const /*t*/,
+                             double const /*delta_t*/) override
+    {
+        unsigned const n_integration_points =
+            _integration_method.getNumberOfPoints();
+
+        for (unsigned ip = 0; ip < n_integration_points; ip++)
+        {
+            _ip_data[ip].pushBackState();
+        }
+    }
+
+    void computeSecondaryVariableConcrete(
+        double const t, std::vector<double> const& local_x) override;
+    void postNonLinearSolverConcrete(std::vector<double> const& local_x,
+                                     double const t,
+                                     bool const use_monolithic_scheme) override;
+
+    Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
+        const unsigned integration_point) const override
+    {
+        auto const& N_u = _secondary_data.N_u[integration_point];
+
+        // assumes N is stored contiguously in memory
+        return Eigen::Map<const Eigen::RowVectorXd>(N_u.data(), N_u.size());
+    }
+
+    std::vector<double> const& getIntPtDarcyVelocity(
+        const double t,
+        GlobalVector const& current_solution,
+        NumLib::LocalToGlobalIndexMap const& dof_table,
+        std::vector<double>& cache) const override;
+
+private:
+    std::size_t setSigma(double const* values)
+    {
+        auto const kelvin_vector_size =
+            MathLib::KelvinVector::KelvinVectorDimensions<
+                DisplacementDim>::value;
+        unsigned const n_integration_points =
+            _integration_method.getNumberOfPoints();
+
+        std::vector<double> ip_sigma_values;
+        auto sigma_values =
+            Eigen::Map<Eigen::Matrix<double, kelvin_vector_size, Eigen::Dynamic,
+                                     Eigen::ColMajor> const>(
+                values, kelvin_vector_size, n_integration_points);
+
+        for (unsigned ip = 0; ip < n_integration_points; ++ip)
+        {
+            _ip_data[ip].sigma =
+                MathLib::KelvinVector::symmetricTensorToKelvinVector(
+                    sigma_values.col(ip));
+        }
+
+        return n_integration_points;
+    }
+
+    // TODO (naumov) This method is same as getIntPtSigma but for arguments and
+    // the ordering of the cache_mat.
+    // There should be only one.
+    std::vector<double> getSigma() const override
+    {
+        auto const kelvin_vector_size =
+            MathLib::KelvinVector::KelvinVectorDimensions<
+                DisplacementDim>::value;
+        unsigned const n_integration_points =
+            _integration_method.getNumberOfPoints();
+
+        std::vector<double> ip_sigma_values;
+        auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
+            double, Eigen::Dynamic, kelvin_vector_size, Eigen::RowMajor>>(
+            ip_sigma_values, n_integration_points, kelvin_vector_size);
+
+        for (unsigned ip = 0; ip < n_integration_points; ++ip)
+        {
+            auto const& sigma = _ip_data[ip].sigma_eff;
+            cache_mat.row(ip) =
+                MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma);
+        }
+
+        return ip_sigma_values;
+    }
+
+    std::vector<double> const& getIntPtSigma(
+        const double /*t*/,
+        GlobalVector const& /*current_solution*/,
+        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+        std::vector<double>& cache) const override
+    {
+        static const int kelvin_vector_size =
+            MathLib::KelvinVector::KelvinVectorDimensions<
+                DisplacementDim>::value;
+        unsigned const n_integration_points =
+            _integration_method.getNumberOfPoints();
+
+        cache.clear();
+        auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
+            double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
+            cache, kelvin_vector_size, n_integration_points);
+
+        for (unsigned ip = 0; ip < n_integration_points; ++ip)
+        {
+            auto const& sigma = _ip_data[ip].sigma_eff;
+            cache_mat.col(ip) =
+                MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma);
+        }
+
+        return cache;
+    }
+
+    virtual std::vector<double> const& getIntPtEpsilon(
+        const double /*t*/,
+        GlobalVector const& /*current_solution*/,
+        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+        std::vector<double>& cache) const override
+    {
+        auto const kelvin_vector_size =
+            MathLib::KelvinVector::KelvinVectorDimensions<
+                DisplacementDim>::value;
+        unsigned const n_integration_points =
+            _integration_method.getNumberOfPoints();
+
+        cache.clear();
+        auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
+            double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
+            cache, kelvin_vector_size, n_integration_points);
+
+        for (unsigned ip = 0; ip < n_integration_points; ++ip)
+        {
+            auto const& eps = _ip_data[ip].eps;
+            cache_mat.col(ip) =
+                MathLib::KelvinVector::kelvinVectorToSymmetricTensor(eps);
+        }
+
+        return cache;
+    }
+
+private:
+    ThermoHydroMechanicsProcessData<DisplacementDim>& _process_data;
+
+    using BMatricesType =
+        BMatrixPolicyType<ShapeFunctionDisplacement, DisplacementDim>;
+    using IpData =
+        IntegrationPointData<BMatricesType, ShapeMatricesTypeDisplacement,
+                             ShapeMatricesTypePressure, DisplacementDim,
+                             ShapeFunctionDisplacement::NPOINTS>;
+    std::vector<IpData, Eigen::aligned_allocator<IpData>> _ip_data;
+
+    IntegrationMethod _integration_method;
+    MeshLib::Element const& _element;
+    bool const _is_axially_symmetric;
+    SecondaryData<
+        typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType>
+        _secondary_data;
+
+    // The shape function of pressure has the same form with the shape function
+    // of temperature
+    static const int temperature_index = 0;
+    static const int temperature_size = ShapeFunctionPressure::NPOINTS;
+    static const int pressure_index = ShapeFunctionPressure::NPOINTS;
+    static const int pressure_size = ShapeFunctionPressure::NPOINTS;
+    static const int displacement_index = ShapeFunctionPressure::NPOINTS * 2;
+    static const int displacement_size =
+        ShapeFunctionDisplacement::NPOINTS * DisplacementDim;
+};
+
+}  // namespace ThermoHydroMechanics
+}  // namespace ProcessLib
+
+#include "ThermoHydroMechanicsFEM-impl.h"
diff --git a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsProcess.cpp b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsProcess.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..bf7ba62e719ac51246ed4b387338e90b7ce7ebf6
--- /dev/null
+++ b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsProcess.cpp
@@ -0,0 +1,412 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "ThermoHydroMechanicsProcess.h"
+
+#include <cassert>
+
+#include "MeshLib/Elements/Utils.h"
+#include "NumLib/DOF/ComputeSparsityPattern.h"
+#include "ProcessLib/Process.h"
+#include "ProcessLib/ThermoHydroMechanics/CreateLocalAssemblers.h"
+
+#include "ThermoHydroMechanicsFEM.h"
+#include "ThermoHydroMechanicsProcessData.h"
+
+namespace ProcessLib
+{
+namespace ThermoHydroMechanics
+{
+template <int DisplacementDim>
+ThermoHydroMechanicsProcess<DisplacementDim>::ThermoHydroMechanicsProcess(
+    MeshLib::Mesh& mesh,
+    std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+    std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
+    unsigned const integration_order,
+    std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
+        process_variables,
+    ThermoHydroMechanicsProcessData<DisplacementDim>&& process_data,
+    SecondaryVariableCollection&& secondary_variables,
+    NumLib::NamedFunctionCaller&& named_function_caller,
+    bool const use_monolithic_scheme)
+    : Process(mesh, std::move(jacobian_assembler), parameters,
+              integration_order, std::move(process_variables),
+              std::move(secondary_variables), std::move(named_function_caller),
+              use_monolithic_scheme),
+      _process_data(std::move(process_data))
+{
+    _nodal_forces = MeshLib::getOrCreateMeshProperty<double>(
+        mesh, "NodalForces", MeshLib::MeshItemType::Node, DisplacementDim);
+
+    _hydraulic_flow = MeshLib::getOrCreateMeshProperty<double>(
+        mesh, "HydraulicFlow", MeshLib::MeshItemType::Node, 1);
+}
+
+template <int DisplacementDim>
+bool ThermoHydroMechanicsProcess<DisplacementDim>::isLinear() const
+{
+    return false;
+}
+
+template <int DisplacementDim>
+MathLib::MatrixSpecifications
+ThermoHydroMechanicsProcess<DisplacementDim>::getMatrixSpecifications(
+    const int process_id) const
+{
+    // For the monolithic scheme or the M process (deformation) in the staggered
+    // scheme.
+    if (_use_monolithic_scheme || process_id == 2)
+    {
+        auto const& l = *_local_to_global_index_map;
+        return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
+                &l.getGhostIndices(), &this->_sparsity_pattern};
+    }
+
+    // For staggered scheme and T or H process (pressure).
+    auto const& l = *_local_to_global_index_map_with_base_nodes;
+    return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
+            &l.getGhostIndices(), &_sparsity_pattern_with_linear_element};
+}
+
+template <int DisplacementDim>
+void ThermoHydroMechanicsProcess<DisplacementDim>::constructDofTable()
+{
+    // Create single component dof in every of the mesh's nodes.
+    _mesh_subset_all_nodes =
+        std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh.getNodes());
+    // Create single component dof in the mesh's base nodes.
+    _base_nodes = MeshLib::getBaseNodes(_mesh.getElements());
+    _mesh_subset_base_nodes =
+        std::make_unique<MeshLib::MeshSubset>(_mesh, _base_nodes);
+
+    // TODO move the two data members somewhere else.
+    // for extrapolation of secondary variables of stress or strain
+    std::vector<MeshLib::MeshSubset> all_mesh_subsets_single_component{
+        *_mesh_subset_all_nodes};
+    _local_to_global_index_map_single_component =
+        std::make_unique<NumLib::LocalToGlobalIndexMap>(
+            std::move(all_mesh_subsets_single_component),
+            // by location order is needed for output
+            NumLib::ComponentOrder::BY_LOCATION);
+
+    if (_use_monolithic_scheme)
+    {
+        // For temperature, which is the first
+        std::vector<MeshLib::MeshSubset> all_mesh_subsets{
+            *_mesh_subset_base_nodes};
+
+        // For pressure, which is the second
+        all_mesh_subsets.push_back(*_mesh_subset_base_nodes);
+
+        // For displacement.
+        const int monolithic_process_id = 0;
+        std::generate_n(std::back_inserter(all_mesh_subsets),
+                        getProcessVariables(monolithic_process_id)[2]
+                            .get()
+                            .getNumberOfComponents(),
+                        [&]() { return *_mesh_subset_all_nodes; });
+
+        std::vector<int> const vec_n_components{1, 1, DisplacementDim};
+        _local_to_global_index_map =
+            std::make_unique<NumLib::LocalToGlobalIndexMap>(
+                std::move(all_mesh_subsets), vec_n_components,
+                NumLib::ComponentOrder::BY_LOCATION);
+        assert(_local_to_global_index_map);
+    }
+    else
+    {
+        // For displacement equation.
+        const int process_id = 2;
+        std::vector<MeshLib::MeshSubset> all_mesh_subsets;
+        std::generate_n(
+            std::back_inserter(all_mesh_subsets),
+            getProcessVariables(process_id)[0].get().getNumberOfComponents(),
+            [&]() { return *_mesh_subset_all_nodes; });
+
+        std::vector<int> const vec_n_components{DisplacementDim};
+        _local_to_global_index_map =
+            std::make_unique<NumLib::LocalToGlobalIndexMap>(
+                std::move(all_mesh_subsets), vec_n_components,
+                NumLib::ComponentOrder::BY_LOCATION);
+
+        // For pressure equation or temperature equation.
+        // Collect the mesh subsets with base nodes in a vector.
+        std::vector<MeshLib::MeshSubset> all_mesh_subsets_base_nodes{
+            *_mesh_subset_base_nodes};
+        _local_to_global_index_map_with_base_nodes =
+            std::make_unique<NumLib::LocalToGlobalIndexMap>(
+                std::move(all_mesh_subsets_base_nodes),
+                // by location order is needed for output
+                NumLib::ComponentOrder::BY_LOCATION);
+
+        _sparsity_pattern_with_linear_element = NumLib::computeSparsityPattern(
+            *_local_to_global_index_map_with_base_nodes, _mesh);
+
+        assert(_local_to_global_index_map);
+        assert(_local_to_global_index_map_with_base_nodes);
+    }
+}
+
+template <int DisplacementDim>
+void ThermoHydroMechanicsProcess<DisplacementDim>::initializeConcreteProcess(
+    NumLib::LocalToGlobalIndexMap const& dof_table,
+    MeshLib::Mesh const& mesh,
+    unsigned const integration_order)
+{
+    const int mechanical_process_id = _use_monolithic_scheme ? 0 : 2;
+    const int deformation_variable_id = _use_monolithic_scheme ? 2 : 0;
+    ProcessLib::ThermoHydroMechanics::createLocalAssemblers<
+        DisplacementDim, ThermoHydroMechanicsLocalAssembler>(
+        mesh.getDimension(), mesh.getElements(), dof_table,
+        // use displacement process variable to set shape function order
+        getProcessVariables(mechanical_process_id)[deformation_variable_id]
+            .get()
+            .getShapeFunctionOrder(),
+        _local_assemblers, mesh.isAxiallySymmetric(), integration_order,
+        _process_data);
+
+    _secondary_variables.addSecondaryVariable(
+        "sigma",
+        makeExtrapolator(MathLib::KelvinVector::KelvinVectorType<
+                             DisplacementDim>::RowsAtCompileTime,
+                         getExtrapolator(), _local_assemblers,
+                         &LocalAssemblerInterface::getIntPtSigma));
+
+    _secondary_variables.addSecondaryVariable(
+        "epsilon",
+        makeExtrapolator(MathLib::KelvinVector::KelvinVectorType<
+                             DisplacementDim>::RowsAtCompileTime,
+                         getExtrapolator(), _local_assemblers,
+                         &LocalAssemblerInterface::getIntPtEpsilon));
+
+    _secondary_variables.addSecondaryVariable(
+        "velocity",
+        makeExtrapolator(mesh.getDimension(), getExtrapolator(),
+                         _local_assemblers,
+                         &LocalAssemblerInterface::getIntPtDarcyVelocity));
+
+    _process_data.pressure_interpolated =
+        MeshLib::getOrCreateMeshProperty<double>(
+            const_cast<MeshLib::Mesh&>(mesh), "pressure_interpolated",
+            MeshLib::MeshItemType::Node, 1);
+
+    _process_data.temperature_interpolated =
+        MeshLib::getOrCreateMeshProperty<double>(
+            const_cast<MeshLib::Mesh&>(mesh), "temperature_interpolated",
+            MeshLib::MeshItemType::Node, 1);
+}
+
+template <int DisplacementDim>
+void ThermoHydroMechanicsProcess<
+    DisplacementDim>::initializeBoundaryConditions()
+{
+    if (_use_monolithic_scheme)
+    {
+        const int process_id_of_thermohydromechancs = 0;
+        initializeProcessBoundaryConditionsAndSourceTerms(
+            *_local_to_global_index_map, process_id_of_thermohydromechancs);
+        return;
+    }
+
+    // Staggered scheme:
+    // for the equations of heat transport
+    const int thermal_process_id = 0;
+    initializeProcessBoundaryConditionsAndSourceTerms(
+        *_local_to_global_index_map_with_base_nodes, thermal_process_id);
+
+    // for the equations of mass balance
+    const int hydraulic_process_id = 1;
+    initializeProcessBoundaryConditionsAndSourceTerms(
+        *_local_to_global_index_map_with_base_nodes, hydraulic_process_id);
+
+    // for the equations of deformation.
+    const int mechanical_process_id = 2;
+    initializeProcessBoundaryConditionsAndSourceTerms(
+        *_local_to_global_index_map, mechanical_process_id);
+}
+
+template <int DisplacementDim>
+void ThermoHydroMechanicsProcess<DisplacementDim>::assembleConcreteProcess(
+    const double t, GlobalVector const& x, GlobalMatrix& M, GlobalMatrix& K,
+    GlobalVector& b)
+{
+    DBUG("Assemble the equations for ThermoHydroMechanics");
+
+    // Note: This assembly function is for the Picard nonlinear solver. Since
+    // only the Newton-Raphson method is employed to simulate coupled HM
+    // processes in this class, this function is actually not used so far.
+
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+        dof_table = {std::ref(*_local_to_global_index_map)};
+    // Call global assembler for each local assembly item.
+    GlobalExecutor::executeMemberDereferenced(
+        _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
+        dof_table, t, x, M, K, b, _coupled_solutions);
+}
+
+template <int DisplacementDim>
+void ThermoHydroMechanicsProcess<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)
+{
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+        dof_tables;
+    // For the monolithic scheme
+    if (_use_monolithic_scheme)
+    {
+        DBUG(
+            "Assemble the Jacobian of ThermoHydroMechanics for the monolithic"
+            " scheme.");
+        dof_tables.emplace_back(*_local_to_global_index_map);
+    }
+    else
+    {
+        // For the staggered scheme
+        if (_coupled_solutions->process_id == 0)
+        {
+            DBUG(
+                "Assemble the Jacobian equations of heat transport process in "
+                "ThermoHydroMechanics for the staggered scheme.");
+        }
+        else if (_coupled_solutions->process_id == 1)
+        {
+            DBUG(
+                "Assemble the Jacobian equations of liquid fluid process in "
+                "ThermoHydroMechanics for the staggered scheme.");
+        }
+        else
+        {
+            DBUG(
+                "Assemble the Jacobian equations of mechanical process in "
+                "ThermoHydroMechanics for the staggered scheme.");
+        }
+        dof_tables.emplace_back(*_local_to_global_index_map_with_base_nodes);
+        dof_tables.emplace_back(*_local_to_global_index_map_with_base_nodes);
+        dof_tables.emplace_back(*_local_to_global_index_map);
+    }
+
+    GlobalExecutor::executeMemberDereferenced(
+        _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
+        _local_assemblers, dof_tables, t, x, xdot, dxdot_dx, dx_dx, M, K, b,
+        Jac, _coupled_solutions);
+
+    auto copyRhs = [&](int const variable_id, auto& output_vector) {
+        if (_use_monolithic_scheme)
+        {
+            transformVariableFromGlobalVector(b, variable_id, dof_tables[0],
+                                              output_vector,
+                                              std::negate<double>());
+        }
+        else
+        {
+            transformVariableFromGlobalVector(
+                b, 0, dof_tables[_coupled_solutions->process_id], output_vector,
+                std::negate<double>());
+        }
+    };
+    if (_use_monolithic_scheme || _coupled_solutions->process_id == 1)
+    {
+        copyRhs(0, *_hydraulic_flow);
+    }
+    if (_use_monolithic_scheme || _coupled_solutions->process_id == 2)
+    {
+        copyRhs(1, *_nodal_forces);
+    }
+}
+
+template <int DisplacementDim>
+void ThermoHydroMechanicsProcess<DisplacementDim>::preTimestepConcreteProcess(
+    GlobalVector const& x, double const t, double const dt,
+    const int process_id)
+{
+    DBUG("PreTimestep ThermoHydroMechanicsProcess.");
+
+    _process_data.dt = dt;
+    _process_data.t = t;
+
+    if (hasMechanicalProcess(process_id))
+        GlobalExecutor::executeMemberOnDereferenced(
+            &LocalAssemblerInterface::preTimestep, _local_assemblers,
+            *_local_to_global_index_map, x, t, dt);
+}
+
+template <int DisplacementDim>
+void ThermoHydroMechanicsProcess<DisplacementDim>::postTimestepConcreteProcess(
+    GlobalVector const& x, const double /*t*/, const double /*delta_t*/,
+    const int process_id)
+{
+    DBUG("PostTimestep ThermoHydroMechanicsProcess.");
+    GlobalExecutor::executeMemberOnDereferenced(
+        &LocalAssemblerInterface::postTimestep, _local_assemblers,
+        getDOFTable(process_id), x);
+}
+
+template <int DisplacementDim>
+void ThermoHydroMechanicsProcess<
+    DisplacementDim>::postNonLinearSolverConcreteProcess(GlobalVector const& x,
+                                                         const double t,
+                                                         const int process_id)
+{
+    if (!hasMechanicalProcess(process_id))
+    {
+        return;
+    }
+
+    DBUG("PostNonLinearSolver HydroMechanicsProcess.");
+    // Calculate strain, stress or other internal variables of mechanics.
+    GlobalExecutor::executeMemberOnDereferenced(
+        &LocalAssemblerInterface::postNonLinearSolver, _local_assemblers,
+        getDOFTable(process_id), x, t, _use_monolithic_scheme);
+}
+
+template <int DisplacementDim>
+void ThermoHydroMechanicsProcess<
+    DisplacementDim>::computeSecondaryVariableConcrete(const double t,
+                                                       GlobalVector const& x,
+                                                       const int process_id)
+{
+    DBUG("Compute the secondary variables for HydroMechanicsProcess.");
+    GlobalExecutor::executeMemberOnDereferenced(
+        &LocalAssemblerInterface::computeSecondaryVariable, _local_assemblers,
+        getDOFTable(process_id), t, x, _coupled_solutions);
+}
+
+template <int DisplacementDim>
+std::tuple<NumLib::LocalToGlobalIndexMap*, bool> ThermoHydroMechanicsProcess<
+    DisplacementDim>::getDOFTableForExtrapolatorData() const
+{
+    const bool manage_storage = false;
+    return std::make_tuple(_local_to_global_index_map_single_component.get(),
+                           manage_storage);
+}
+
+template <int DisplacementDim>
+NumLib::LocalToGlobalIndexMap const&
+ThermoHydroMechanicsProcess<DisplacementDim>::getDOFTable(
+    const int process_id) const
+{
+    if (hasMechanicalProcess(process_id))
+    {
+        return *_local_to_global_index_map;
+    }
+
+    // For the equation of pressure
+    return *_local_to_global_index_map_with_base_nodes;
+}
+
+template class ThermoHydroMechanicsProcess<2>;
+template class ThermoHydroMechanicsProcess<3>;
+
+}  // namespace ThermoHydroMechanics
+}  // namespace ProcessLib
diff --git a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsProcess.h b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsProcess.h
new file mode 100644
index 0000000000000000000000000000000000000000..06e715747b012b63942b07fcee6eda1b78657d33
--- /dev/null
+++ b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsProcess.h
@@ -0,0 +1,147 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2019, 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 "LocalAssemblerInterface.h"
+#include "ProcessLib/Process.h"
+#include "ThermoHydroMechanicsProcessData.h"
+
+namespace ProcessLib
+{
+namespace ThermoHydroMechanics
+{
+struct LocalAssemblerInterface;
+
+/// Thermally induced deformation process in linear kinematics
+/// poro-mechanical/biphasic model.
+///
+/// The mixture momentum balance, the mixture mass balance and the mixture
+/// energy balance are solved under fully saturated conditions.
+template <int DisplacementDim>
+class ThermoHydroMechanicsProcess final : public Process
+{
+public:
+    ThermoHydroMechanicsProcess(
+        MeshLib::Mesh& mesh,
+        std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&&
+            jacobian_assembler,
+        std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const&
+            parameters,
+        unsigned const integration_order,
+        std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
+            process_variables,
+        ThermoHydroMechanicsProcessData<DisplacementDim>&& process_data,
+        SecondaryVariableCollection&& secondary_variables,
+        NumLib::NamedFunctionCaller&& named_function_caller,
+        bool const use_monolithic_scheme);
+
+    //! \name ODESystem interface
+    //! @{
+
+    bool isLinear() const override;
+    //! @}
+
+    /**
+     * Get the size and the sparse pattern of the global matrix in order to
+     * create the global matrices and vectors for the system equations of this
+     * process.
+     *
+     * @param process_id Process ID. If the monolithic scheme is applied,
+     *                               process_id = 0. For the staggered scheme,
+     *                               process_id = 0 represents the
+     *                               hydraulic (H) process, while process_id = 1
+     *                               represents the mechanical (M) process.
+     * @return Matrix specifications including size and sparse pattern.
+     */
+    MathLib::MatrixSpecifications getMatrixSpecifications(
+        const int process_id) const override;
+
+private:
+    void constructDofTable() override;
+
+    void initializeConcreteProcess(
+        NumLib::LocalToGlobalIndexMap const& dof_table,
+        MeshLib::Mesh const& mesh,
+        unsigned const integration_order) override;
+
+    void initializeBoundaryConditions() override;
+
+    void assembleConcreteProcess(const double t, GlobalVector const& x,
+                                 GlobalMatrix& M, GlobalMatrix& K,
+                                 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;
+
+    void preTimestepConcreteProcess(GlobalVector const& x, double const t,
+                                    double const dt,
+                                    const int process_id) override;
+
+    void postTimestepConcreteProcess(GlobalVector const& x, const double t,
+                                     const double delta_t,
+                                     int const process_id) override;
+
+    void postNonLinearSolverConcreteProcess(GlobalVector const& x,
+                                            const double t,
+                                            int const process_id) override;
+
+    NumLib::LocalToGlobalIndexMap const& getDOFTable(
+        const int process_id) const override;
+
+private:
+    std::vector<MeshLib::Node*> _base_nodes;
+    std::unique_ptr<MeshLib::MeshSubset const> _mesh_subset_base_nodes;
+    ThermoHydroMechanicsProcessData<DisplacementDim> _process_data;
+
+    std::vector<std::unique_ptr<LocalAssemblerInterface>> _local_assemblers;
+
+    std::unique_ptr<NumLib::LocalToGlobalIndexMap>
+        _local_to_global_index_map_single_component;
+
+    /// Local to global index mapping for base nodes, which is used for linear
+    /// interpolation for pressure in the staggered scheme.
+    std::unique_ptr<NumLib::LocalToGlobalIndexMap>
+        _local_to_global_index_map_with_base_nodes;
+
+    /// Sparsity pattern for the flow equation, and it is initialized only if
+    /// the staggered scheme is used.
+    GlobalSparsityPattern _sparsity_pattern_with_linear_element;
+
+    /// Solutions of the previous time step
+    std::array<std::unique_ptr<GlobalVector>, 2> _xs_previous_timestep;
+
+    void computeSecondaryVariableConcrete(const double t, GlobalVector const& x,
+                                          const int process_id) override;
+    /**
+     * @copydoc ProcessLib::Process::getDOFTableForExtrapolatorData()
+     */
+    std::tuple<NumLib::LocalToGlobalIndexMap*, bool>
+    getDOFTableForExtrapolatorData() const override;
+
+    /// Check whether the process represented by \c process_id is/has
+    /// mechanical process. In the present implementation, the mechanical
+    /// process has process_id == 1 in the staggered scheme.
+    bool hasMechanicalProcess(int const process_id) const
+    {
+        return _use_monolithic_scheme || process_id == 2;
+    }
+
+    MeshLib::PropertyVector<double>* _nodal_forces = nullptr;
+    MeshLib::PropertyVector<double>* _hydraulic_flow = nullptr;
+};
+
+extern template class ThermoHydroMechanicsProcess<2>;
+extern template class ThermoHydroMechanicsProcess<3>;
+
+}  // namespace ThermoHydroMechanics
+}  // namespace ProcessLib
diff --git a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsProcessData.h b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsProcessData.h
new file mode 100644
index 0000000000000000000000000000000000000000..5e3e0ef80e5457b31ecc49fda1378804dab491b5
--- /dev/null
+++ b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsProcessData.h
@@ -0,0 +1,142 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2019, 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 "ParameterLib/Parameter.h"
+
+#include <memory>
+#include <utility>
+
+#include <Eigen/Dense>
+
+namespace MaterialLib
+{
+namespace Solids
+{
+template <int DisplacementDim>
+struct MechanicsBase;
+}
+}  // namespace MaterialLib
+namespace ProcessLib
+{
+namespace ThermoHydroMechanics
+{
+template <int DisplacementDim>
+struct ThermoHydroMechanicsProcessData
+{
+    ThermoHydroMechanicsProcessData(
+        MeshLib::PropertyVector<int> const* const material_ids_,
+        std::map<int,
+                 std::unique_ptr<
+                     MaterialLib::Solids::MechanicsBase<DisplacementDim>>>&&
+            solid_materials_,
+        ParameterLib::Parameter<double> const& intrinsic_permeability_,
+        ParameterLib::Parameter<double> const& specific_storage_,
+        ParameterLib::Parameter<double> const& fluid_viscosity_,
+        ParameterLib::Parameter<double> const& fluid_density_,
+        ParameterLib::Parameter<double> const& biot_coefficient_,
+        ParameterLib::Parameter<double> const& porosity_,
+        ParameterLib::Parameter<double> const& solid_density_,
+        ParameterLib::Parameter<double> const&
+            solid_linear_thermal_expansion_coefficient_,
+        ParameterLib::Parameter<double> const&
+            fluid_volumetric_thermal_expansion_coefficient_,
+        ParameterLib::Parameter<double> const& solid_specific_heat_capacity_,
+        ParameterLib::Parameter<double> const& fluid_specific_heat_capacity_,
+        ParameterLib::Parameter<double> const& solid_thermal_conductivity_,
+        ParameterLib::Parameter<double> const& fluid_thermal_conductivity_,
+        ParameterLib::Parameter<double> const& reference_temperature_,
+        Eigen::Matrix<double, DisplacementDim, 1>
+            specific_body_force_)
+        : material_ids(material_ids_),
+          solid_materials{std::move(solid_materials_)},
+          intrinsic_permeability(intrinsic_permeability_),
+          specific_storage(specific_storage_),
+          fluid_viscosity(fluid_viscosity_),
+          fluid_density(fluid_density_),
+          biot_coefficient(biot_coefficient_),
+          porosity(porosity_),
+          solid_density(solid_density_),
+          solid_linear_thermal_expansion_coefficient(
+              solid_linear_thermal_expansion_coefficient_),
+          fluid_volumetric_thermal_expansion_coefficient(
+              fluid_volumetric_thermal_expansion_coefficient_),
+          solid_specific_heat_capacity(solid_specific_heat_capacity_),
+          solid_thermal_conductivity(solid_thermal_conductivity_),
+          fluid_specific_heat_capacity(fluid_specific_heat_capacity_),
+          fluid_thermal_conductivity(fluid_thermal_conductivity_),
+          reference_temperature(reference_temperature_),
+          specific_body_force(std::move(specific_body_force_))
+    {
+    }
+
+    ThermoHydroMechanicsProcessData(ThermoHydroMechanicsProcessData&& other) =
+        default;
+
+    //! Copies are forbidden.
+    ThermoHydroMechanicsProcessData(ThermoHydroMechanicsProcessData const&) =
+        delete;
+
+    //! Assignments are not needed.
+    void operator=(ThermoHydroMechanicsProcessData const&) = delete;
+
+    //! Assignments are not needed.
+    void operator=(ThermoHydroMechanicsProcessData&&) = delete;
+
+    MeshLib::PropertyVector<int> const* const material_ids;
+
+    /// The constitutive relation for the mechanical part.
+    /// \note Linear elasticity is the only supported one in the moment.
+    std::map<
+        int,
+        std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>>
+        solid_materials;
+    /// Permeability of the solid. A scalar quantity,
+    ///	ParameterLib::Parameter<double>.
+    ParameterLib::Parameter<double> const& intrinsic_permeability;
+    /// Volumetric average specific storage of the solid and fluid phases.
+    /// A scalar quantity, ParameterLib::Parameter<double>.
+    ParameterLib::Parameter<double> const& specific_storage;
+    /// Fluid's viscosity. A scalar quantity, ParameterLib::Parameter<double>.
+    ParameterLib::Parameter<double> const& fluid_viscosity;
+    /// Fluid's density. A scalar quantity, ParameterLib::Parameter<double>.
+    ParameterLib::Parameter<double> const& fluid_density;
+    /// Biot coefficient. A scalar quantity, ParameterLib::Parameter<double>.
+    ParameterLib::Parameter<double> const& biot_coefficient;
+    /// Porosity of the solid. A scalar quantity,
+    /// ParameterLib::Parameter<double>.
+    ParameterLib::Parameter<double> const& porosity;
+    /// Solid's density. A scalar quantity, ParameterLib::Parameter<double>.
+    ParameterLib::Parameter<double> const& solid_density;
+    ParameterLib::Parameter<double> const&
+        solid_linear_thermal_expansion_coefficient;
+    ParameterLib::Parameter<double> const&
+        fluid_volumetric_thermal_expansion_coefficient;
+    ParameterLib::Parameter<double> const& solid_specific_heat_capacity;
+    ParameterLib::Parameter<double> const& solid_thermal_conductivity;
+    ParameterLib::Parameter<double> const& fluid_specific_heat_capacity;
+    ParameterLib::Parameter<double> const& fluid_thermal_conductivity;
+    ParameterLib::Parameter<double> const& reference_temperature;
+    /// Specific body forces applied to solid and fluid.
+    /// It is usually used to apply gravitational forces.
+    /// A vector of displacement dimension's length.
+    Eigen::Matrix<double, DisplacementDim, 1> const specific_body_force;
+    double dt = 0.0;
+    double t = 0.0;
+
+    MeshLib::PropertyVector<double>* pressure_interpolated = nullptr;
+    MeshLib::PropertyVector<double>* temperature_interpolated = nullptr;
+
+    EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
+};
+
+}  // namespace ThermoHydroMechanics
+}  // namespace ProcessLib
diff --git a/Tests/Data/ThermoHydroMechanics/Linear/Beam_sealed_bimaterial/beam_2nd.vtu b/Tests/Data/ThermoHydroMechanics/Linear/Beam_sealed_bimaterial/beam_2nd.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..97046413c5b007a4a85a9a6653957555bacb643e
--- /dev/null
+++ b/Tests/Data/ThermoHydroMechanics/Linear/Beam_sealed_bimaterial/beam_2nd.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:d56b531eb0a81888fe525e17fd08b087bde80b4c15d3b54d58616e6cf771cc05
+size 192434
diff --git a/Tests/Data/ThermoHydroMechanics/Linear/Beam_sealed_bimaterial/expected_square_1e2_pcs_0_ts_10_t_100.000000.vtu b/Tests/Data/ThermoHydroMechanics/Linear/Beam_sealed_bimaterial/expected_square_1e2_pcs_0_ts_10_t_100.000000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..fdd10da3b3a37c01b1bf1004e5a3a38bed2bf314
--- /dev/null
+++ b/Tests/Data/ThermoHydroMechanics/Linear/Beam_sealed_bimaterial/expected_square_1e2_pcs_0_ts_10_t_100.000000.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:ae640efd0d16acc0cfd7175e9d692690fc9b151fecb869a21a83881f2c3842c4
+size 356872
diff --git a/Tests/Data/ThermoHydroMechanics/Linear/Beam_sealed_bimaterial/square_1e2.prj b/Tests/Data/ThermoHydroMechanics/Linear/Beam_sealed_bimaterial/square_1e2.prj
new file mode 100644
index 0000000000000000000000000000000000000000..8f84cb29607644e3aee95655b85ac068742f76e0
--- /dev/null
+++ b/Tests/Data/ThermoHydroMechanics/Linear/Beam_sealed_bimaterial/square_1e2.prj
@@ -0,0 +1,323 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+    <mesh axially_symmetric="true">beam_2nd.vtu</mesh>
+    <geometry>square_1x1.gml</geometry>
+    <processes>
+        <process>
+            <name>THERMO_HYDRO_MECHANICS</name>
+            <type>THERMO_HYDRO_MECHANICS</type>
+            <integration_order>3</integration_order>
+            <dimension>2</dimension>
+            <constitutive_relation>
+                <type>LinearElasticIsotropic</type>
+                <youngs_modulus>E</youngs_modulus>
+                <poissons_ratio>nu</poissons_ratio>
+            </constitutive_relation>
+            <intrinsic_permeability>k</intrinsic_permeability>
+            <specific_storage>S</specific_storage>
+            <biot_coefficient>alpha</biot_coefficient>
+            <porosity>phi</porosity>
+            <solid_density>rho_sr</solid_density>
+            <fluid_density>rho_fr</fluid_density>
+            <fluid_viscosity>mu</fluid_viscosity>
+            <fluid_volumetric_thermal_expansion_coefficient>beta_f</fluid_volumetric_thermal_expansion_coefficient>
+            <solid_linear_thermal_expansion_coefficient>alpha_s</solid_linear_thermal_expansion_coefficient>
+            <solid_specific_heat_capacity>C_s</solid_specific_heat_capacity>
+            <solid_thermal_conductivity>lambda_s</solid_thermal_conductivity>
+            <fluid_specific_heat_capacity>C_f</fluid_specific_heat_capacity>
+            <fluid_thermal_conductivity>lambda_f</fluid_thermal_conductivity>
+            <reference_temperature>T0</reference_temperature>
+            <process_variables>
+                <displacement>displacement</displacement>
+                <pressure>pressure</pressure>
+                <temperature>temperature</temperature>
+            </process_variables>
+            <secondary_variables>
+                <secondary_variable type="static" internal_name="sigma" output_name="sigma"/>
+                <secondary_variable type="static" internal_name="epsilon" output_name="epsilon"/>
+            </secondary_variables>
+            <specific_body_force>0 0</specific_body_force>
+        </process>
+    </processes>
+    <time_loop>
+        <processes>
+            <process ref="THERMO_HYDRO_MECHANICS">
+                <nonlinear_solver>basic_newton</nonlinear_solver>
+                <convergence_criterion>
+                    <type>Residual</type>
+                    <norm_type>NORM2</norm_type>
+                    <abstol>1e-5</abstol>
+                </convergence_criterion>
+                <time_discretization>
+                    <type>BackwardEuler</type>
+                </time_discretization>
+                <time_stepping>
+                    <type>FixedTimeStepping</type>
+                    <t_initial>0</t_initial>
+                    <t_end>100</t_end>
+                    <timesteps>
+                        <pair>
+                            <repeat>100</repeat>
+                            <delta_t>10</delta_t>
+                        </pair>
+                    </timesteps>
+                </time_stepping>
+            </process>
+        </processes>
+        <output>
+            <type>VTK</type>
+            <prefix>square_1e2</prefix>
+            <timesteps>
+                <pair>
+                    <repeat>1</repeat>
+                    <each_steps>10</each_steps>
+                </pair>
+            </timesteps>
+            <variables>
+                <variable>displacement</variable>
+                <variable>pressure</variable>
+                <variable>temperature</variable>
+                <variable>sigma</variable>
+                <variable>epsilon</variable>
+            </variables>
+        </output>
+    </time_loop>
+    <parameters>
+        <parameter>
+            <name>E</name>
+            <type>Group</type>
+            <group_id_property>MaterialIDs</group_id_property>
+            <index_values>
+                <index>0</index>
+                <value>19</value>
+            </index_values>
+            <index_values>
+                <index>1</index>
+                <value>0.95</value>
+            </index_values>
+            <index_values>
+                <index>2</index>
+                <value>9.5</value>
+            </index_values>
+        </parameter>
+        <parameter>
+            <name>alpha_s</name>
+            <type>Group</type>
+            <group_id_property>MaterialIDs</group_id_property>
+            <index_values>
+                <index>0</index>
+                <value>0.23333e-5</value>
+            </index_values>
+            <index_values>
+                <index>1</index>
+                <value>0.23333e-5</value>
+            </index_values>
+            <index_values>
+                <index>2</index>
+                <value>0.46667e-5</value>
+            </index_values>
+        </parameter>
+        <parameter>
+            <name>phi</name>
+            <type>Group</type>
+            <group_id_property>MaterialIDs</group_id_property>
+            <index_values>
+                <index>0</index>
+                <value>0.1</value>
+            </index_values>
+            <index_values>
+                <index>1</index>
+                <value>0</value>
+            </index_values>
+            <index_values>
+                <index>2</index>
+                <value>0.1</value>
+            </index_values>
+        </parameter>
+        <parameter>
+            <name>k</name>
+            <type>Group</type>
+            <group_id_property>MaterialIDs</group_id_property>
+            <index_values>
+                <index>0</index>
+                <value>1e-6</value>
+            </index_values>
+            <index_values>
+                <index>1</index>
+                <value>1e-16</value>
+            </index_values>
+            <index_values>
+                <index>2</index>
+                <value>1e-6</value>
+            </index_values>
+        </parameter>
+        <parameter>
+            <name>alpha</name>
+            <type>Group</type>
+            <group_id_property>MaterialIDs</group_id_property>
+            <index_values>
+                <index>0</index>
+                <value>1</value>
+            </index_values>
+            <index_values>
+                <index>1</index>
+                <value>0</value>
+            </index_values>
+            <index_values>
+                <index>2</index>
+                <value>1</value>
+            </index_values>
+        </parameter>
+        <parameter>
+            <name>nu</name>
+            <type>Constant</type>
+            <value>0.3</value>
+        </parameter>
+        <parameter>
+            <name>S</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+        <parameter>
+            <name>rho_sr</name>
+            <type>Constant</type>
+            <value>1.75e-6</value>
+        </parameter>
+        <parameter>
+            <name>beta_f</name>
+            <type>Constant</type>
+            <value>2.07e-4</value>
+        </parameter>
+        <parameter>
+            <name>lambda_s</name>
+            <type>Constant</type>
+            <value>2.5e-3</value>
+        </parameter>
+        <parameter>
+            <name>C_s</name>
+            <type>Constant</type>
+            <value>1714.0</value>
+        </parameter>
+        <parameter>
+            <name>T0</name>
+            <type>Constant</type>
+            <value>283.15</value>
+        </parameter>
+        <parameter>
+            <name>displacement0</name>
+            <type>Constant</type>
+            <values>0 0</values>
+        </parameter>
+        <parameter>
+            <name>pressure_ic</name>
+            <type>Constant</type>
+            <values>0</values>
+        </parameter>
+        <parameter>
+            <name>dirichlet0</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+        <parameter>
+            <name>temperature_ic</name>
+            <type>Constant</type>
+            <value>283.15</value>
+        </parameter>
+        <parameter>
+            <name>temperature_bc_left</name>
+            <type>Constant</type>
+            <value>363.15</value>
+        </parameter>
+        <parameter>
+            <name>rho_fr</name>
+            <type>Constant</type>
+            <value>1e-6</value>
+        </parameter>
+        <parameter>
+            <name>mu</name>
+            <type>Constant</type>
+            <value>1.278e-9</value>
+        </parameter>
+        <parameter>
+            <name>C_f</name>
+            <type>Constant</type>
+            <value>4186</value>
+        </parameter>
+        <parameter>
+            <name>lambda_f</name>
+            <type>Constant</type>
+            <value>0.56e-3</value>
+        </parameter>
+    </parameters>
+    <process_variables>
+        <process_variable>
+            <name>displacement</name>
+            <components>2</components>
+            <order>2</order>
+            <initial_condition>displacement0</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <geometrical_set>square_1x1_geometry</geometrical_set>
+                    <geometry>left</geometry>
+                    <type>Dirichlet</type>
+                    <component>0</component>
+                    <parameter>dirichlet0</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>square_1x1_geometry</geometrical_set>
+                    <geometry>bottom</geometry>
+                    <type>Dirichlet</type>
+                    <component>1</component>
+                    <parameter>dirichlet0</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+        <process_variable>
+            <name>pressure</name>
+            <components>1</components>
+            <order>1</order>
+            <initial_condition>pressure_ic</initial_condition>
+            <boundary_conditions>
+            </boundary_conditions>
+        </process_variable>
+        <process_variable>
+            <name>temperature</name>
+            <components>1</components>
+            <order>1</order>
+            <initial_condition>temperature_ic</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <geometrical_set>square_1x1_geometry</geometrical_set>
+                    <geometry>left</geometry>
+                    <type>Dirichlet</type>
+                    <component>0</component>
+                    <parameter>temperature_bc_left</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+    </process_variables>
+    <nonlinear_solvers>
+        <nonlinear_solver>
+            <name>basic_newton</name>
+            <type>Newton</type>
+            <max_iter>50</max_iter>
+            <linear_solver>general_linear_solver</linear_solver>
+        </nonlinear_solver>
+    </nonlinear_solvers>
+    <linear_solvers>
+        <linear_solver>
+            <name>general_linear_solver</name>
+            <lis>-i bicgstab -p ilu -tol 1e-16 -maxiter 10000</lis>
+            <eigen>
+                <solver_type>BiCGSTAB</solver_type>
+                <precon_type>ILUT</precon_type>
+                <max_iteration_step>10000</max_iteration_step>
+                <error_tolerance>1e-16</error_tolerance>
+            </eigen>
+            <petsc>
+                <prefix>sd</prefix>
+                <parameters>-sd_ksp_type cg -sd_pc_type bjacobi -sd_ksp_rtol 1e-16 -sd_ksp_max_it 10000</parameters>
+            </petsc>
+        </linear_solver>
+    </linear_solvers>
+</OpenGeoSysProject>
diff --git a/Tests/Data/ThermoHydroMechanics/Linear/Beam_sealed_bimaterial/square_1x1.gml b/Tests/Data/ThermoHydroMechanics/Linear/Beam_sealed_bimaterial/square_1x1.gml
new file mode 100644
index 0000000000000000000000000000000000000000..472140a5202012a22af1b5875a2f1f72dc82be15
--- /dev/null
+++ b/Tests/Data/ThermoHydroMechanics/Linear/Beam_sealed_bimaterial/square_1x1.gml
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:cc39b7e701db9413add64e21ebdca8d5192de5996d4b00b43848fcedad979b0e
+size 1133
diff --git a/Tests/Data/ThermoHydroMechanics/Linear/Beam_unsealed_bimaterial/beam_2nd.vtu b/Tests/Data/ThermoHydroMechanics/Linear/Beam_unsealed_bimaterial/beam_2nd.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..97046413c5b007a4a85a9a6653957555bacb643e
--- /dev/null
+++ b/Tests/Data/ThermoHydroMechanics/Linear/Beam_unsealed_bimaterial/beam_2nd.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:d56b531eb0a81888fe525e17fd08b087bde80b4c15d3b54d58616e6cf771cc05
+size 192434
diff --git a/Tests/Data/ThermoHydroMechanics/Linear/Beam_unsealed_bimaterial/expected_square_1e2_pcs_0_ts_10_t_1000.000000.vtu b/Tests/Data/ThermoHydroMechanics/Linear/Beam_unsealed_bimaterial/expected_square_1e2_pcs_0_ts_10_t_1000.000000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..e01d0d066e501c91a1105467b00601ef82bb1493
--- /dev/null
+++ b/Tests/Data/ThermoHydroMechanics/Linear/Beam_unsealed_bimaterial/expected_square_1e2_pcs_0_ts_10_t_1000.000000.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:0e79fa4019cf033b3f93f65cb56b37a0666e7a89ceaf510b9e9f17159c28a82a
+size 335560
diff --git a/Tests/Data/ThermoHydroMechanics/Linear/Beam_unsealed_bimaterial/square_1e2.prj b/Tests/Data/ThermoHydroMechanics/Linear/Beam_unsealed_bimaterial/square_1e2.prj
new file mode 100644
index 0000000000000000000000000000000000000000..e5ae4709a9984e82ccae571959b37a592f6abaa6
--- /dev/null
+++ b/Tests/Data/ThermoHydroMechanics/Linear/Beam_unsealed_bimaterial/square_1e2.prj
@@ -0,0 +1,328 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+    <mesh axially_symmetric="true">beam_2nd.vtu</mesh>
+    <geometry>square_1x1.gml</geometry>
+    <processes>
+        <process>
+            <name>THERMO_HYDRO_MECHANICS</name>
+            <type>THERMO_HYDRO_MECHANICS</type>
+            <integration_order>3</integration_order>
+            <dimension>2</dimension>
+            <constitutive_relation>
+                <type>LinearElasticIsotropic</type>
+                <youngs_modulus>E</youngs_modulus>
+                <poissons_ratio>nu</poissons_ratio>
+            </constitutive_relation>
+            <intrinsic_permeability>k</intrinsic_permeability>
+            <specific_storage>S</specific_storage>
+            <biot_coefficient>alpha</biot_coefficient>
+            <porosity>phi</porosity>
+            <solid_density>rho_sr</solid_density>
+            <fluid_density>rho_fr</fluid_density>
+            <fluid_viscosity>mu</fluid_viscosity>
+            <fluid_volumetric_thermal_expansion_coefficient>beta_f</fluid_volumetric_thermal_expansion_coefficient>
+            <solid_linear_thermal_expansion_coefficient>alpha_s</solid_linear_thermal_expansion_coefficient>
+            <solid_specific_heat_capacity>C_s</solid_specific_heat_capacity>
+            <solid_thermal_conductivity>lambda_s</solid_thermal_conductivity>
+            <fluid_specific_heat_capacity>C_f</fluid_specific_heat_capacity>
+            <fluid_thermal_conductivity>lambda_f</fluid_thermal_conductivity>
+            <reference_temperature>T0</reference_temperature>
+            <process_variables>
+                <displacement>displacement</displacement>
+                <pressure>pressure</pressure>
+                <temperature>temperature</temperature>
+            </process_variables>
+            <secondary_variables>
+                <secondary_variable type="static" internal_name="sigma" output_name="sigma"/>
+                <secondary_variable type="static" internal_name="epsilon" output_name="epsilon"/>
+            </secondary_variables>
+            <specific_body_force>0 0</specific_body_force>
+        </process>
+    </processes>
+    <time_loop>
+        <processes>
+            <process ref="THERMO_HYDRO_MECHANICS">
+                <nonlinear_solver>basic_newton</nonlinear_solver>
+                <convergence_criterion>
+                    <type>Residual</type>
+                    <norm_type>NORM2</norm_type>
+                    <abstol>1e-5</abstol>
+                </convergence_criterion>
+                <time_discretization>
+                    <type>BackwardEuler</type>
+                </time_discretization>
+                <time_stepping>
+                    <type>FixedTimeStepping</type>
+                    <t_initial>0</t_initial>
+                    <t_end>1000</t_end>
+                    <timesteps>
+                        <pair>
+                            <repeat>10</repeat>
+                            <delta_t>100</delta_t>
+                        </pair>
+                    </timesteps>
+                </time_stepping>
+            </process>
+        </processes>
+        <output>
+            <type>VTK</type>
+            <prefix>square_1e2</prefix>
+            <timesteps>
+                <pair>
+                    <repeat>1</repeat>
+                    <each_steps>10</each_steps>
+                </pair>
+            </timesteps>
+            <variables>
+                <variable>displacement</variable>
+                <variable>pressure</variable>
+                <variable>temperature</variable>
+                <variable>sigma</variable>
+                <variable>epsilon</variable>
+            </variables>
+        </output>
+    </time_loop>
+    <parameters>
+        <parameter>
+            <name>E</name>
+            <type>Group</type>
+            <group_id_property>MaterialIDs</group_id_property>
+            <index_values>
+                <index>0</index>
+                <value>20</value>
+            </index_values>
+            <index_values>
+                <index>1</index>
+                <value>0.1</value>
+            </index_values>
+            <index_values>
+                <index>2</index>
+                <value>10</value>
+            </index_values>
+        </parameter>
+        <parameter>
+            <name>alpha_s</name>
+            <type>Group</type>
+            <group_id_property>MaterialIDs</group_id_property>
+            <index_values>
+                <index>0</index>
+                <value>0.2333333333e-5</value>
+            </index_values>
+            <index_values>
+                <index>1</index>
+                <value>0.2333333333e-5</value>
+            </index_values>
+            <index_values>
+                <index>2</index>
+                <value>0.4666666667e-5</value>
+            </index_values>
+        </parameter>
+        <parameter>
+            <name>phi</name>
+            <type>Group</type>
+            <group_id_property>MaterialIDs</group_id_property>
+            <index_values>
+                <index>0</index>
+                <value>0.1</value>
+            </index_values>
+            <index_values>
+                <index>1</index>
+                <value>0</value>
+            </index_values>
+            <index_values>
+                <index>2</index>
+                <value>0.1</value>
+            </index_values>
+        </parameter>
+        <parameter>
+            <name>k</name>
+            <type>Group</type>
+            <group_id_property>MaterialIDs</group_id_property>
+            <index_values>
+                <index>0</index>
+                <value>1e-6</value>
+            </index_values>
+            <index_values>
+                <index>1</index>
+                <value>1e-6</value>
+            </index_values>
+            <index_values>
+                <index>2</index>
+                <value>1e-6</value>
+            </index_values>
+        </parameter>
+        <parameter>
+            <name>alpha</name>
+            <type>Group</type>
+            <group_id_property>MaterialIDs</group_id_property>
+            <index_values>
+                <index>0</index>
+                <value>1</value>
+            </index_values>
+            <index_values>
+                <index>1</index>
+                <value>0</value>
+            </index_values>
+            <index_values>
+                <index>2</index>
+                <value>1</value>
+            </index_values>
+        </parameter>
+        <parameter>
+            <name>nu</name>
+            <type>Constant</type>
+            <value>0.3</value>
+        </parameter>
+        <parameter>
+            <name>S</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+        <parameter>
+            <name>rho_sr</name>
+            <type>Constant</type>
+            <value>1.75e-6</value>
+        </parameter>
+        <parameter>
+            <name>beta_f</name>
+            <type>Constant</type>
+            <value>2.07e-4</value>
+        </parameter>
+        <parameter>
+            <name>lambda_s</name>
+            <type>Constant</type>
+            <value>2.5e-3</value>
+        </parameter>
+        <parameter>
+            <name>C_s</name>
+            <type>Constant</type>
+            <value>1714.0</value>
+        </parameter>
+        <parameter>
+            <name>T0</name>
+            <type>Constant</type>
+            <value>283.15</value>
+        </parameter>
+        <parameter>
+            <name>displacement0</name>
+            <type>Constant</type>
+            <values>0 0</values>
+        </parameter>
+        <parameter>
+            <name>pressure_ic</name>
+            <type>Constant</type>
+            <values>0</values>
+        </parameter>
+        <parameter>
+            <name>dirichlet0</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+        <parameter>
+            <name>Neumann0</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+        <parameter>
+            <name>temperature_ic</name>
+            <type>Constant</type>
+            <value>283.15</value>
+        </parameter>
+        <parameter>
+            <name>temperature_bc_left</name>
+            <type>Constant</type>
+            <value>363.15</value>
+        </parameter>
+        <parameter>
+            <name>rho_fr</name>
+            <type>Constant</type>
+            <value>1e-6</value>
+        </parameter>
+        <parameter>
+            <name>mu</name>
+            <type>Constant</type>
+            <value>1.278e-9</value>
+        </parameter>
+        <parameter>
+            <name>C_f</name>
+            <type>Constant</type>
+            <value>4186</value>
+        </parameter>
+        <parameter>
+            <name>lambda_f</name>
+            <type>Constant</type>
+            <value>0.56</value>
+        </parameter>
+    </parameters>
+    <process_variables>
+        <process_variable>
+            <name>displacement</name>
+            <components>2</components>
+            <order>2</order>
+            <initial_condition>displacement0</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <geometrical_set>square_1x1_geometry</geometrical_set>
+                    <geometry>left</geometry>
+                    <type>Dirichlet</type>
+                    <component>0</component>
+                    <parameter>dirichlet0</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>square_1x1_geometry</geometrical_set>
+                    <geometry>bottom</geometry>
+                    <type>Dirichlet</type>
+                    <component>1</component>
+                    <parameter>dirichlet0</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+        <process_variable>
+            <name>pressure</name>
+            <components>1</components>
+            <order>1</order>
+            <initial_condition>pressure_ic</initial_condition>
+            <boundary_conditions>
+            </boundary_conditions>
+        </process_variable>
+        <process_variable>
+            <name>temperature</name>
+            <components>1</components>
+            <order>1</order>
+            <initial_condition>temperature_ic</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <geometrical_set>square_1x1_geometry</geometrical_set>
+                    <geometry>left</geometry>
+                    <type>Dirichlet</type>
+                    <component>0</component>
+                    <parameter>temperature_bc_left</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+    </process_variables>
+    <nonlinear_solvers>
+        <nonlinear_solver>
+            <name>basic_newton</name>
+            <type>Newton</type>
+            <max_iter>50</max_iter>
+            <linear_solver>general_linear_solver</linear_solver>
+        </nonlinear_solver>
+    </nonlinear_solvers>
+    <linear_solvers>
+        <linear_solver>
+            <name>general_linear_solver</name>
+            <lis>-i bicgstab -p ilu -tol 1e-16 -maxiter 10000</lis>
+            <eigen>
+                <solver_type>BiCGSTAB</solver_type>
+                <precon_type>ILUT</precon_type>
+                <max_iteration_step>10000</max_iteration_step>
+                <error_tolerance>1e-16</error_tolerance>
+            </eigen>
+            <petsc>
+                <prefix>sd</prefix>
+                <parameters>-sd_ksp_type cg -sd_pc_type bjacobi -sd_ksp_rtol 1e-16 -sd_ksp_max_it 10000</parameters>
+            </petsc>
+        </linear_solver>
+    </linear_solvers>
+</OpenGeoSysProject>
diff --git a/Tests/Data/ThermoHydroMechanics/Linear/Beam_unsealed_bimaterial/square_1x1.gml b/Tests/Data/ThermoHydroMechanics/Linear/Beam_unsealed_bimaterial/square_1x1.gml
new file mode 100644
index 0000000000000000000000000000000000000000..472140a5202012a22af1b5875a2f1f72dc82be15
--- /dev/null
+++ b/Tests/Data/ThermoHydroMechanics/Linear/Beam_unsealed_bimaterial/square_1x1.gml
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:cc39b7e701db9413add64e21ebdca8d5192de5996d4b00b43848fcedad979b0e
+size 1133
diff --git a/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/expected_square_1e0_pcs_0_ts_10_t_1000.000000.vtu b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/expected_square_1e0_pcs_0_ts_10_t_1000.000000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..04d4e90ea8364ffed808ac8e361e32b70da8aa4e
--- /dev/null
+++ b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/expected_square_1e0_pcs_0_ts_10_t_1000.000000.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:04bb17ef4747b70c8f3df046cbba73bd16c02077dc271973fb9bcc43616cf3e9
+size 457220
diff --git a/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/expected_square_1e0_pcs_0_ts_10_t_50000.000000.vtu b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/expected_square_1e0_pcs_0_ts_10_t_50000.000000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..405dda331fd71fe77e76a1cef3da3ae6f2d8307d
--- /dev/null
+++ b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/expected_square_1e0_pcs_0_ts_10_t_50000.000000.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:56ed62cdfe5f19514b193eac3919f66c94647d9a0c354d9b31d0ac61253ce164
+size 460562
diff --git a/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/quarter_002_2nd.vtu b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/quarter_002_2nd.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..b2f12bf3f2d9a122a55b0fe1e16a07cd0ca4e100
--- /dev/null
+++ b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/quarter_002_2nd.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:875fd6c188f4a5b59b9b341c7cc908d105c2d8155c35ffe0cd1985e68eb0de04
+size 185692
diff --git a/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/square_1e2.prj b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/square_1e2.prj
new file mode 100644
index 0000000000000000000000000000000000000000..159678752a13ff78d5403d48d92e63a752933384
--- /dev/null
+++ b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/square_1e2.prj
@@ -0,0 +1,293 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+    <mesh axially_symmetric="true">quarter_002_2nd.vtu</mesh>
+    <geometry>square_1x1.gml</geometry>
+    <processes>
+        <process>
+            <name>THERMO_HYDRO_MECHANICS</name>
+            <type>THERMO_HYDRO_MECHANICS</type>
+            <integration_order>4</integration_order>
+            <dimension>2</dimension>
+            <constitutive_relation>
+                <type>LinearElasticIsotropic</type>
+                <youngs_modulus>E</youngs_modulus>
+                <poissons_ratio>nu</poissons_ratio>
+            </constitutive_relation>
+            <intrinsic_permeability>k</intrinsic_permeability>
+            <specific_storage>S</specific_storage>
+            <biot_coefficient>alpha</biot_coefficient>
+            <porosity>phi</porosity>
+            <solid_density>rho_sr</solid_density>
+            <fluid_density>rho_fr</fluid_density>
+            <fluid_viscosity>mu</fluid_viscosity>
+            <fluid_volumetric_thermal_expansion_coefficient>beta_f</fluid_volumetric_thermal_expansion_coefficient>
+            <solid_linear_thermal_expansion_coefficient>alpha_s</solid_linear_thermal_expansion_coefficient>
+            <solid_specific_heat_capacity>C_s</solid_specific_heat_capacity>
+            <solid_thermal_conductivity>lambda_s</solid_thermal_conductivity>
+            <fluid_specific_heat_capacity>C_f</fluid_specific_heat_capacity>
+            <fluid_thermal_conductivity>lambda_f</fluid_thermal_conductivity>
+            <reference_temperature>T0</reference_temperature>
+            <process_variables>
+                <displacement>displacement</displacement>
+                <pressure>pressure</pressure>
+                <temperature>temperature</temperature>
+            </process_variables>
+            <secondary_variables>
+                <secondary_variable type="static" internal_name="sigma" output_name="sigma"/>
+                <secondary_variable type="static" internal_name="epsilon" output_name="epsilon"/>
+            </secondary_variables>
+            <specific_body_force>0 0</specific_body_force>
+        </process>
+    </processes>
+    <time_loop>
+        <processes>
+            <process ref="THERMO_HYDRO_MECHANICS">
+                <nonlinear_solver>basic_newton</nonlinear_solver>
+                <convergence_criterion>
+                    <type>PerComponentDeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <abstols>1e-5 1e-5 1e-5 1e-5</abstols>
+                </convergence_criterion>
+                <time_discretization>
+                    <type>BackwardEuler</type>
+                </time_discretization>
+                <time_stepping>
+                    <type>FixedTimeStepping</type>
+                    <t_initial>0</t_initial>
+                    <t_end>50000</t_end>
+                    <timesteps>
+                        <pair>
+                            <repeat>10</repeat>
+                            <delta_t>5000</delta_t>
+                        </pair>
+                    </timesteps>
+                </time_stepping>
+            </process>
+        </processes>
+        <output>
+            <type>VTK</type>
+            <prefix>square_1e0</prefix>
+            <timesteps>
+                <pair>
+                    <repeat>1</repeat>
+                    <each_steps>10</each_steps>
+                </pair>
+            </timesteps>
+            <variables>
+                <variable>displacement</variable>
+                <variable>pressure</variable>
+                <variable>temperature</variable>
+                <variable>sigma</variable>
+                <variable>epsilon</variable>
+            </variables>
+        </output>
+    </time_loop>
+    <parameters>
+        <parameter>
+            <name>E</name>
+            <type>Constant</type>
+            <value>5000000000</value>
+        </parameter>
+        <parameter>
+            <name>nu</name>
+            <type>Constant</type>
+            <value>0.3</value>
+        </parameter>
+        <parameter>
+            <name>k</name>
+            <type>Constant</type>
+            <value>2e-20</value>
+        </parameter>
+        <parameter>
+            <name>S</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+        <parameter>
+            <name>alpha</name>
+            <type>Constant</type>
+            <value>1</value>
+        </parameter>
+        <parameter>
+            <name>phi</name>
+            <type>Constant</type>
+            <value>0.16</value>
+        </parameter>
+        <parameter>
+            <name>rho_sr</name>
+            <type>Constant</type>
+            <value>2290</value>
+        </parameter>
+        <parameter>
+            <name>alpha_s</name>
+            <type>Constant</type>
+            <value>1.5e-5</value>
+        </parameter>
+        <parameter>
+            <name>beta_f</name>
+            <type>Constant</type>
+            <value>4e-4</value>
+        </parameter>
+        <parameter>
+            <name>lambda_s</name>
+            <type>Constant</type>
+            <value>1.838</value>
+        </parameter>
+        <parameter>
+            <name>C_s</name>
+            <type>Constant</type>
+            <value>917.654</value>
+        </parameter>
+        <parameter>
+            <name>T0</name>
+            <type>Constant</type>
+            <value>273.15</value>
+        </parameter>
+        <parameter>
+            <name>displacement0</name>
+            <type>Constant</type>
+            <values>0 0</values>
+        </parameter>
+        <parameter>
+            <name>pressure_ic</name>
+            <type>Constant</type>
+            <values>0</values>
+        </parameter>
+        <parameter>
+            <name>dirichlet0</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+        <parameter>
+            <name>Neumann0</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+        <parameter>
+            <name>temperature_ic</name>
+            <type>Constant</type>
+            <value>273.15</value>
+        </parameter>
+        <parameter>
+            <name>pressure_bc_left</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+        <parameter>
+            <name>temperature_bc_left</name>
+            <type>Constant</type>
+            <value>273.15</value>
+        </parameter>
+        <parameter>
+            <name>temperature_source_term</name>
+            <type>Constant</type>
+            <value>150</value>
+        </parameter>
+        <parameter>
+            <name>rho_fr</name>
+            <type>Constant</type>
+            <value>999.1</value>
+        </parameter>
+        <parameter>
+            <name>mu</name>
+            <type>Constant</type>
+            <value>1e-3</value>
+        </parameter>
+        <parameter>
+            <name>C_f</name>
+            <type>Constant</type>
+            <value>4280</value>
+        </parameter>
+        <parameter>
+            <name>lambda_f</name>
+            <type>Constant</type>
+            <value>0.6</value>
+        </parameter>
+    </parameters>
+    <process_variables>
+        <process_variable>
+            <name>displacement</name>
+            <components>2</components>
+            <order>2</order>
+            <initial_condition>displacement0</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <geometrical_set>square_1x1_geometry</geometrical_set>
+                    <geometry>bottom</geometry>
+                    <type>Dirichlet</type>
+                    <component>1</component>
+                    <parameter>dirichlet0</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>square_1x1_geometry</geometrical_set>
+                    <geometry>left</geometry>
+                    <type>Dirichlet</type>
+                    <component>0</component>
+                    <parameter>dirichlet0</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+        <process_variable>
+            <name>pressure</name>
+            <components>1</components>
+            <order>1</order>
+            <initial_condition>pressure_ic</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <geometrical_set>square_1x1_geometry</geometrical_set>
+                    <geometry>out</geometry>
+                    <type>Dirichlet</type>
+                    <component>0</component>
+                    <parameter>pressure_bc_left</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+        <process_variable>
+            <name>temperature</name>
+            <components>1</components>
+            <order>1</order>
+            <initial_condition>temperature_ic</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <geometrical_set>square_1x1_geometry</geometrical_set>
+                    <geometry>out</geometry>
+                    <type>Dirichlet</type>
+                    <component>0</component>
+                    <parameter>temperature_bc_left</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+            <source_terms>
+                <source_term>
+                    <geometrical_set>square_1x1_geometry</geometrical_set>
+                    <geometry>center</geometry>
+                    <type>Nodal</type>
+                    <parameter>temperature_source_term</parameter>
+                </source_term>
+            </source_terms>
+        </process_variable>
+    </process_variables>
+    <nonlinear_solvers>
+        <nonlinear_solver>
+            <name>basic_newton</name>
+            <type>Newton</type>
+            <max_iter>50</max_iter>
+            <linear_solver>general_linear_solver</linear_solver>
+        </nonlinear_solver>
+    </nonlinear_solvers>
+    <linear_solvers>
+        <linear_solver>
+            <name>general_linear_solver</name>
+            <lis>-i bicgstab -p ilu -tol 1e-16 -maxiter 10000</lis>
+            <eigen>
+                <solver_type>BiCGSTAB</solver_type>
+                <precon_type>ILUT</precon_type>
+                <max_iteration_step>10000</max_iteration_step>
+                <error_tolerance>1e-8</error_tolerance>
+            </eigen>
+            <petsc>
+                <prefix>sd</prefix>
+                <parameters>-sd_ksp_type cg -sd_pc_type bjacobi -sd_ksp_rtol 1e-16 -sd_ksp_max_it 10000</parameters>
+            </petsc>
+        </linear_solver>
+    </linear_solvers>
+</OpenGeoSysProject>
diff --git a/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/square_1x1.gml b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/square_1x1.gml
new file mode 100644
index 0000000000000000000000000000000000000000..e06ad003f98d662df50868f579cf42e09bda24ec
--- /dev/null
+++ b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/square_1x1.gml
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:fe58d69d8a678f35e866f0dcadec37aad6c42cbfa12cacc4ec5473494f822e38
+size 2598
diff --git a/Tests/Data/ThermoHydroMechanics/Linear/Square_sealed_homogeneous/expected_square_1e0_pcs_0_ts_10_t_1000.000000.vtu b/Tests/Data/ThermoHydroMechanics/Linear/Square_sealed_homogeneous/expected_square_1e0_pcs_0_ts_10_t_1000.000000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..a862716ac7a3be8ee9bd2b857e9802540fa8e86d
--- /dev/null
+++ b/Tests/Data/ThermoHydroMechanics/Linear/Square_sealed_homogeneous/expected_square_1e0_pcs_0_ts_10_t_1000.000000.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:dcd83e50c02cceec44c7447cf98973fd03eb0faded7a7136670e526bb6466efa
+size 3775
diff --git a/Tests/Data/ThermoHydroMechanics/Linear/Square_sealed_homogeneous/square_1e0.prj b/Tests/Data/ThermoHydroMechanics/Linear/Square_sealed_homogeneous/square_1e0.prj
new file mode 100644
index 0000000000000000000000000000000000000000..727f1e8b899af85de36ec0957eac26889d461303
--- /dev/null
+++ b/Tests/Data/ThermoHydroMechanics/Linear/Square_sealed_homogeneous/square_1e0.prj
@@ -0,0 +1,263 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+    <mesh axially_symmetric="true">square_1x1_quad2_1e0.vtu</mesh>
+    <geometry>square_1x1.gml</geometry>
+    <processes>
+        <process>
+            <name>THERMO_HYDRO_MECHANICS</name>
+            <type>THERMO_HYDRO_MECHANICS</type>
+            <integration_order>4</integration_order>
+            <dimension>2</dimension>
+            <constitutive_relation>
+                <type>LinearElasticIsotropic</type>
+                <youngs_modulus>E</youngs_modulus>
+                <poissons_ratio>nu</poissons_ratio>
+            </constitutive_relation>
+            <intrinsic_permeability>k</intrinsic_permeability>
+            <specific_storage>S</specific_storage>
+            <biot_coefficient>alpha</biot_coefficient>
+            <porosity>phi</porosity>
+            <solid_density>rho_sr</solid_density>
+            <fluid_density>rho_fr</fluid_density>
+            <fluid_viscosity>mu</fluid_viscosity>
+            <fluid_volumetric_thermal_expansion_coefficient>beta_f</fluid_volumetric_thermal_expansion_coefficient>
+            <solid_linear_thermal_expansion_coefficient>alpha_s</solid_linear_thermal_expansion_coefficient>
+            <solid_specific_heat_capacity>C_s</solid_specific_heat_capacity>
+            <solid_thermal_conductivity>lambda_s</solid_thermal_conductivity>
+            <fluid_specific_heat_capacity>C_f</fluid_specific_heat_capacity>
+            <fluid_thermal_conductivity>lambda_f</fluid_thermal_conductivity>
+            <reference_temperature>T0</reference_temperature>
+            <process_variables>
+                <temperature>temperature</temperature>
+                <pressure>pressure</pressure>
+                <displacement>displacement</displacement>
+            </process_variables>
+            <secondary_variables>
+                <secondary_variable type="static" internal_name="sigma" output_name="sigma"/>
+                <secondary_variable type="static" internal_name="epsilon" output_name="epsilon"/>
+            </secondary_variables>
+            <specific_body_force>0 0</specific_body_force>
+        </process>
+    </processes>
+    <time_loop>
+        <processes>
+            <process ref="THERMO_HYDRO_MECHANICS">
+                <nonlinear_solver>basic_newton</nonlinear_solver>
+                <convergence_criterion>
+                    <type>PerComponentDeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <abstols>1e-7 1e-5 1e-5 1e-5</abstols>
+                </convergence_criterion>
+                <time_discretization>
+                    <type>BackwardEuler</type>
+                </time_discretization>
+                <time_stepping>
+                    <type>FixedTimeStepping</type>
+                    <t_initial>0</t_initial>
+                    <t_end>1000</t_end>
+                    <timesteps>
+                        <pair>
+                            <repeat>10</repeat>
+                            <delta_t>100</delta_t>
+                        </pair>
+                    </timesteps>
+                </time_stepping>
+            </process>
+        </processes>
+        <output>
+            <type>VTK</type>
+            <prefix>square_1e0</prefix>
+            <timesteps>
+                <pair>
+                    <repeat>1</repeat>
+                    <each_steps>10</each_steps>
+                </pair>
+            </timesteps>
+            <variables>
+                <variable>displacement</variable>
+                <variable>pressure</variable>
+                <variable>temperature</variable>
+                <variable>sigma</variable>
+                <variable>epsilon</variable>
+            </variables>
+        </output>
+    </time_loop>
+    <parameters>
+        <parameter>
+            <name>E</name>
+            <type>Constant</type>
+            <value>21</value>
+        </parameter>
+        <parameter>
+            <name>nu</name>
+            <type>Constant</type>
+            <value>0.3</value>
+        </parameter>
+        <parameter>
+            <name>k</name>
+            <type>Constant</type>
+            <value>1e-6</value>
+        </parameter>
+        <parameter>
+            <name>S</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+        <parameter>
+            <name>alpha</name>
+            <type>Constant</type>
+            <value>1</value>
+        </parameter>
+        <parameter>
+            <name>phi</name>
+            <type>Constant</type>
+            <value>0.4</value>
+        </parameter>
+        <parameter>
+            <name>rho_sr</name>
+            <type>Constant</type>
+            <value>1.75e-6</value>
+        </parameter>
+        <parameter>
+            <name>alpha_s</name>
+            <type>Constant</type>
+            <value>0.7e-5</value>
+        </parameter>
+        <parameter>
+            <name>beta_f</name>
+            <type>Constant</type>
+            <value>2.07e-4</value>
+        </parameter>
+        <parameter>
+            <name>lambda_s</name>
+            <type>Constant</type>
+            <value>0.9e-3</value>
+        </parameter>
+        <parameter>
+            <name>C_s</name>
+            <type>Constant</type>
+            <value>1714.0</value>
+        </parameter>
+        <parameter>
+            <name>T0</name>
+            <type>Constant</type>
+            <value>273.15</value>
+        </parameter>
+        <parameter>
+            <name>displacement0</name>
+            <type>Constant</type>
+            <values>0 0</values>
+        </parameter>
+        <parameter>
+            <name>pressure_ic</name>
+            <type>Constant</type>
+            <values>0</values>
+        </parameter>
+        <parameter>
+            <name>dirichlet0</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+        <parameter>
+            <name>temperature_ic</name>
+            <type>Constant</type>
+            <value>273.15</value>
+        </parameter>
+        <parameter>
+            <name>temperature_bc_left</name>
+            <type>Constant</type>
+            <value>353.15</value>
+        </parameter>
+        <parameter>
+            <name>rho_fr</name>
+            <type>Constant</type>
+            <value>1e-6</value>
+        </parameter>
+        <parameter>
+            <name>mu</name>
+            <type>Constant</type>
+            <value>1.278e-9</value>
+        </parameter>
+        <parameter>
+            <name>C_f</name>
+            <type>Constant</type>
+            <value>4186</value>
+        </parameter>
+        <parameter>
+            <name>lambda_f</name>
+            <type>Constant</type>
+            <value>0.56</value>
+        </parameter>
+    </parameters>
+    <process_variables>
+        <process_variable>
+            <name>displacement</name>
+            <components>2</components>
+            <order>2</order>
+            <initial_condition>displacement0</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <geometrical_set>square_1x1_geometry</geometrical_set>
+                    <geometry>left</geometry>
+                    <type>Dirichlet</type>
+                    <component>0</component>
+                    <parameter>dirichlet0</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>square_1x1_geometry</geometrical_set>
+                    <geometry>bottom</geometry>
+                    <type>Dirichlet</type>
+                    <component>1</component>
+                    <parameter>dirichlet0</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+        <process_variable>
+            <name>pressure</name>
+            <components>1</components>
+            <order>1</order>
+            <initial_condition>pressure_ic</initial_condition>
+            <boundary_conditions>
+            </boundary_conditions>
+        </process_variable>
+        <process_variable>
+            <name>temperature</name>
+            <components>1</components>
+            <order>1</order>
+            <initial_condition>temperature_ic</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <geometrical_set>square_1x1_geometry</geometrical_set>
+                    <geometry>left</geometry>
+                    <type>Dirichlet</type>
+                    <component>0</component>
+                    <parameter>temperature_bc_left</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+    </process_variables>
+    <nonlinear_solvers>
+        <nonlinear_solver>
+            <name>basic_newton</name>
+            <type>Newton</type>
+            <max_iter>50</max_iter>
+            <linear_solver>general_linear_solver</linear_solver>
+        </nonlinear_solver>
+    </nonlinear_solvers>
+    <linear_solvers>
+        <linear_solver>
+            <name>general_linear_solver</name>
+            <lis>-i bicgstab -p ilu -tol 1e-16 -maxiter 10000</lis>
+            <eigen>
+                <solver_type>BiCGSTAB</solver_type>
+                <precon_type>ILUT</precon_type>
+                <max_iteration_step>10000</max_iteration_step>
+                <error_tolerance>1e-16</error_tolerance>
+            </eigen>
+            <petsc>
+                <prefix>sd</prefix>
+                <parameters>-sd_ksp_type cg -sd_pc_type bjacobi -sd_ksp_rtol 1e-16 -sd_ksp_max_it 10000</parameters>
+            </petsc>
+        </linear_solver>
+    </linear_solvers>
+</OpenGeoSysProject>
diff --git a/Tests/Data/ThermoHydroMechanics/Linear/Square_sealed_homogeneous/square_1x1.gml b/Tests/Data/ThermoHydroMechanics/Linear/Square_sealed_homogeneous/square_1x1.gml
new file mode 100644
index 0000000000000000000000000000000000000000..d47737cc4098963de6a933f2677a99eb5ed34345
--- /dev/null
+++ b/Tests/Data/ThermoHydroMechanics/Linear/Square_sealed_homogeneous/square_1x1.gml
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:6080489c93ce68c1a51c292698ecb7fef60a1d0cc1a3655c90988aa9a7d3b03d
+size 1125
diff --git a/Tests/Data/ThermoHydroMechanics/Linear/Square_sealed_homogeneous/square_1x1_quad2_1e0.vtu b/Tests/Data/ThermoHydroMechanics/Linear/Square_sealed_homogeneous/square_1x1_quad2_1e0.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..832750be45ef5ae5dd500ef6e30bb6eba9468224
--- /dev/null
+++ b/Tests/Data/ThermoHydroMechanics/Linear/Square_sealed_homogeneous/square_1x1_quad2_1e0.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:63e0885eded57de4ed24c0146d5ba3caa54916f2614c696a6c34d3936e52426b
+size 1632