diff --git a/ProcessLib/ThermoRichardsMechanics/CreateThermoRichardsMechanicsProcess.cpp b/ProcessLib/ThermoRichardsMechanics/CreateThermoRichardsMechanicsProcess.cpp
index ed2c2debfd6e673a69c2a2fec070200908bb1caf..241a5035b0df03e68226e4e5e4f326305b572b22 100644
--- a/ProcessLib/ThermoRichardsMechanics/CreateThermoRichardsMechanicsProcess.cpp
+++ b/ProcessLib/ThermoRichardsMechanics/CreateThermoRichardsMechanicsProcess.cpp
@@ -178,13 +178,20 @@ std::unique_ptr<Process> createThermoRichardsMechanicsProcess(
         mass_lumping = *mass_lumping_ptr;
     }
 
+    const bool use_TaylorHood_elements =
+        variable_p->getShapeFunctionOrder() !=
+                variable_u->getShapeFunctionOrder()
+            ? true
+            : false;
+
     ThermoRichardsMechanicsProcessData<DisplacementDim> process_data{
         materialIDs(mesh),
         std::move(media_map),
         std::move(solid_constitutive_relations),
         initial_stress,
         specific_body_force,
-        mass_lumping};
+        mass_lumping,
+        use_TaylorHood_elements};
 
     SecondaryVariableCollection secondary_variables;
 
diff --git a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.cpp b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.cpp
index dd10b1a117c2e4a6fe820c202a31cf0cde9676c4..1d22580762f49e5c80057840e983540c65340e1a 100644
--- a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.cpp
+++ b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.cpp
@@ -109,12 +109,12 @@ template <int DisplacementDim>
 void ThermoRichardsMechanicsProcess<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());
+    _mesh_subset_all_nodes = std::make_unique<MeshLib::MeshSubset>(
+        _mesh, _mesh.getNodes(), process_data_.use_TaylorHood_elements);
     // 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_);
+    mesh_subset_base_nodes_ = std::make_unique<MeshLib::MeshSubset>(
+        _mesh, base_nodes_, process_data_.use_TaylorHood_elements);
 
     // TODO move the two data members somewhere else.
     // for extrapolation of secondary variables of stress or strain
diff --git a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcessData.h b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcessData.h
index ca97f3846ed5a6e83cb0fa3595f0a515b6049604..d495409e2c47f062101ce079ceab827bea07f3f4 100644
--- a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcessData.h
+++ b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcessData.h
@@ -53,6 +53,8 @@ struct ThermoRichardsMechanicsProcessData
 
     bool const apply_mass_lumping;
 
+    const bool use_TaylorHood_elements;
+
     MeshLib::PropertyVector<double>* element_saturation = nullptr;
     MeshLib::PropertyVector<double>* element_porosity = nullptr;
     MeshLib::PropertyVector<double>* element_stresses = nullptr;