diff --git a/MaterialLib/SolidModels/CreateConstitutiveRelation.cpp b/MaterialLib/SolidModels/CreateConstitutiveRelation.cpp
index 41a0002acdbd0eb18a813955d03381eb91587770..1f8ac0c61d0c0189ecd38d12ff4e2d9304c4e537 100644
--- a/MaterialLib/SolidModels/CreateConstitutiveRelation.cpp
+++ b/MaterialLib/SolidModels/CreateConstitutiveRelation.cpp
@@ -33,47 +33,83 @@ createConstitutiveRelation(
     std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
     BaseLib::ConfigTree const& config)
 {
-    auto const constitutive_relation_config =
-        //! \ogs_file_param{material__solid__constitutive_relation}
-        config.getConfigSubtree("constitutive_relation");
-
     auto const type =
         //! \ogs_file_param{material__solid__constitutive_relation__type}
-        constitutive_relation_config.peekConfigParameter<std::string>("type");
+        config.peekConfigParameter<std::string>("type");
 
     if (type == "Ehlers")
     {
         return MaterialLib::Solids::Ehlers::createEhlers<DisplacementDim>(
-            parameters, constitutive_relation_config);
+            parameters, config);
     }
     if (type == "LinearElasticIsotropic")
     {
         const bool skip_type_checking = false;
         return MaterialLib::Solids::createLinearElasticIsotropic<
-            DisplacementDim>(
-            parameters, constitutive_relation_config, skip_type_checking);
+            DisplacementDim>(parameters, config, skip_type_checking);
     }
     if (type == "Lubby2")
     {
         return MaterialLib::Solids::Lubby2::createLubby2<DisplacementDim>(
-            parameters, constitutive_relation_config);
+            parameters, config);
     }
     if (type == "CreepBGRa")
     {
         return MaterialLib::Solids::Creep::createCreepBGRa<DisplacementDim>(
-            parameters, constitutive_relation_config);
+            parameters, config);
     }
     OGS_FATAL("Cannot construct constitutive relation of given type \'%s\'.",
               type.c_str());
 }
 
-template std::unique_ptr<MaterialLib::Solids::MechanicsBase<2>>
-createConstitutiveRelation<2>(
+template <int DisplacementDim>
+std::map<int,
+         std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>>
+createConstitutiveRelations(
+    std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+    BaseLib::ConfigTree const& config)
+{
+    auto const constitutive_relation_configs =
+        //! \ogs_file_param{material__solid__constitutive_relation}
+        config.getConfigSubtreeList("constitutive_relation");
+
+    std::map<int, std::unique_ptr<MechanicsBase<DisplacementDim>>>
+        constitutive_relations;
+
+    for (auto const& constitutive_relation_config :
+         constitutive_relation_configs)
+    {
+        int const material_id =
+            constitutive_relation_config.getConfigAttribute<int>("id", 0);
+
+        if (constitutive_relations.find(material_id) !=
+            constitutive_relations.end())
+        {
+            OGS_FATAL(
+                "Multiple constitutive relations were specified for the same "
+                "material id %d. Keep in mind, that if no material id is "
+                "specified, it is assumed to be 0 by default.",
+                material_id);
+        }
+
+        constitutive_relations.emplace(
+            material_id,
+            createConstitutiveRelation<DisplacementDim>(
+                parameters, constitutive_relation_config));
+    }
+
+    DBUG("Found %d constitutive relations.", constitutive_relations.size());
+
+    return constitutive_relations;
+}
+
+template std::map<int, std::unique_ptr<MaterialLib::Solids::MechanicsBase<2>>>
+createConstitutiveRelations<2>(
     std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
     BaseLib::ConfigTree const& config);
 
-template std::unique_ptr<MaterialLib::Solids::MechanicsBase<3>>
-createConstitutiveRelation<3>(
+template std::map<int, std::unique_ptr<MaterialLib::Solids::MechanicsBase<3>>>
+createConstitutiveRelations<3>(
     std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
     BaseLib::ConfigTree const& config);
 }
diff --git a/MaterialLib/SolidModels/CreateConstitutiveRelation.h b/MaterialLib/SolidModels/CreateConstitutiveRelation.h
index 15c43ebf81a20513d59894b0c8be6b0a036049b6..e1e4bba9421f49941a6d9e2bc41c1cfda87a96ea 100644
--- a/MaterialLib/SolidModels/CreateConstitutiveRelation.h
+++ b/MaterialLib/SolidModels/CreateConstitutiveRelation.h
@@ -11,6 +11,7 @@
 
 #pragma once
 
+#include <map>
 #include <memory>
 #include <vector>
 
@@ -32,18 +33,21 @@ template <int DisplacementDim>
 struct MechanicsBase;
 
 template <int DisplacementDim>
-std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>
-createConstitutiveRelation(
+std::map<int,
+         std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>>
+createConstitutiveRelations(
     std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
     BaseLib::ConfigTree const& config);
 
-extern template std::unique_ptr<MaterialLib::Solids::MechanicsBase<2>>
-createConstitutiveRelation<2>(
+extern template std::map<int,
+                         std::unique_ptr<MaterialLib::Solids::MechanicsBase<2>>>
+createConstitutiveRelations<2>(
     std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
     BaseLib::ConfigTree const& config);
 
-extern template std::unique_ptr<MaterialLib::Solids::MechanicsBase<3>>
-createConstitutiveRelation<3>(
+extern template std::map<int,
+                         std::unique_ptr<MaterialLib::Solids::MechanicsBase<3>>>
+createConstitutiveRelations<3>(
     std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
     BaseLib::ConfigTree const& config);
 }
diff --git a/MaterialLib/SolidModels/CreepBGRa.cpp b/MaterialLib/SolidModels/CreepBGRa.cpp
index 78da877423ef33e63e2acc59b751cdf6f2612ad0..b59f6d85341dbb7558204c82f198616ba8243aa6 100644
--- a/MaterialLib/SolidModels/CreepBGRa.cpp
+++ b/MaterialLib/SolidModels/CreepBGRa.cpp
@@ -39,7 +39,7 @@ CreepBGRa<DisplacementDim>::integrateStress(
     KelvinVector const& sigma_prev,
     typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
     /*material_state_variables*/,
-    double const T)
+    double const T) const
 {
     using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
 
diff --git a/MaterialLib/SolidModels/CreepBGRa.h b/MaterialLib/SolidModels/CreepBGRa.h
index ee29d0e23c4d75249d39d2f34302534aea1f5772..affcbd6091f589ecd05d0c562c7e20e68262c516 100644
--- a/MaterialLib/SolidModels/CreepBGRa.h
+++ b/MaterialLib/SolidModels/CreepBGRa.h
@@ -55,7 +55,7 @@ public:
 
     std::unique_ptr<
         typename MechanicsBase<DisplacementDim>::MaterialStateVariables>
-    createMaterialStateVariables() override
+    createMaterialStateVariables() const override
     {
         return LinearElasticIsotropic<
             DisplacementDim>::createMaterialStateVariables();
@@ -75,17 +75,17 @@ public:
     {
     }
 
-    boost::optional<
-        std::tuple<KelvinVector, std::unique_ptr<typename MechanicsBase<
-                                     DisplacementDim>::MaterialStateVariables>,
-                   KelvinMatrix>>
+    boost::optional<std::tuple<KelvinVector,
+                               std::unique_ptr<typename MechanicsBase<
+                                   DisplacementDim>::MaterialStateVariables>,
+                               KelvinMatrix>>
     integrateStress(
         double const t, ProcessLib::SpatialPosition const& x, double const dt,
         KelvinVector const& eps_prev, KelvinVector const& eps,
         KelvinVector const& sigma_prev,
         typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
             material_state_variables,
-        double const T) override;
+        double const T) const override;
 
     ConstitutiveModel getConstitutiveModel() const override
     {
diff --git a/MaterialLib/SolidModels/Ehlers.cpp b/MaterialLib/SolidModels/Ehlers.cpp
index 23e0f9aef5a9ce77aeae9f40376f22891a01cfd2..680dcb246a3b803223a561ec6589eec011e2a207 100644
--- a/MaterialLib/SolidModels/Ehlers.cpp
+++ b/MaterialLib/SolidModels/Ehlers.cpp
@@ -473,14 +473,12 @@ boost::optional<std::tuple<typename SolidEhlers<DisplacementDim>::KelvinVector,
                                DisplacementDim>::MaterialStateVariables>,
                            typename SolidEhlers<DisplacementDim>::KelvinMatrix>>
 SolidEhlers<DisplacementDim>::integrateStress(
-    double const t,
-    ProcessLib::SpatialPosition const& x,
-    double const dt,
-    KelvinVector const& eps_prev,
-    KelvinVector const& eps,
+    double const t, ProcessLib::SpatialPosition const& x, double const dt,
+    KelvinVector const& eps_prev, KelvinVector const& eps,
     KelvinVector const& sigma_prev,
     typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
-        material_state_variables, double const /*T*/)
+        material_state_variables,
+    double const /*T*/) const
 {
     assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
                &material_state_variables) != nullptr);
diff --git a/MaterialLib/SolidModels/Ehlers.h b/MaterialLib/SolidModels/Ehlers.h
index 0be72d6005d89439e6a8db9adf1901415b1fe352..2958c1c94ed4150a3fb145091f9863c839d2538f 100644
--- a/MaterialLib/SolidModels/Ehlers.h
+++ b/MaterialLib/SolidModels/Ehlers.h
@@ -267,7 +267,7 @@ public:
 public:
     std::unique_ptr<
         typename MechanicsBase<DisplacementDim>::MaterialStateVariables>
-    createMaterialStateVariables() override
+    createMaterialStateVariables() const override
     {
         return std::make_unique<StateVariables<DisplacementDim>>();
     }
@@ -309,19 +309,17 @@ public:
         return (eps - eps_p.D - eps_p.V / 3 * identity2).dot(sigma) / 2;
     }
 
-    boost::optional<
-        std::tuple<KelvinVector, std::unique_ptr<typename MechanicsBase<
-                                     DisplacementDim>::MaterialStateVariables>,
-                   KelvinMatrix>>
+    boost::optional<std::tuple<KelvinVector,
+                               std::unique_ptr<typename MechanicsBase<
+                                   DisplacementDim>::MaterialStateVariables>,
+                               KelvinMatrix>>
     integrateStress(
-        double const t,
-        ProcessLib::SpatialPosition const& x,
-        double const dt,
-        KelvinVector const& eps_prev,
-        KelvinVector const& eps,
+        double const t, ProcessLib::SpatialPosition const& x, double const dt,
+        KelvinVector const& eps_prev, KelvinVector const& eps,
         KelvinVector const& sigma_prev,
         typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
-            material_state_variables, double const T) override;
+            material_state_variables,
+        double const T) const override;
 
     std::vector<typename MechanicsBase<DisplacementDim>::InternalVariable>
     getInternalVariables() const override;
diff --git a/MaterialLib/SolidModels/LinearElasticIsotropic.cpp b/MaterialLib/SolidModels/LinearElasticIsotropic.cpp
index 14d2c7bdcab8d92b5467f9131bebb43652e27799..e28133f1c425211b6c55bc8a9494df9d72286e1c 100644
--- a/MaterialLib/SolidModels/LinearElasticIsotropic.cpp
+++ b/MaterialLib/SolidModels/LinearElasticIsotropic.cpp
@@ -20,14 +20,12 @@ boost::optional<
                    DisplacementDim>::MaterialStateVariables>,
                typename LinearElasticIsotropic<DisplacementDim>::KelvinMatrix>>
 LinearElasticIsotropic<DisplacementDim>::integrateStress(
-    double const t,
-    ProcessLib::SpatialPosition const& x,
-    double const /*dt*/,
-    KelvinVector const& eps_prev,
-    KelvinVector const& eps,
+    double const t, ProcessLib::SpatialPosition const& x, double const /*dt*/,
+    KelvinVector const& eps_prev, KelvinVector const& eps,
     KelvinVector const& sigma_prev,
     typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
-        material_state_variables, double const T)
+        material_state_variables,
+    double const T) const
 {
     KelvinMatrix C = getElasticTensor(t, x, T);
 
diff --git a/MaterialLib/SolidModels/LinearElasticIsotropic.h b/MaterialLib/SolidModels/LinearElasticIsotropic.h
index 46bf4554f57416a38088cf7a693018f2d07112f3..2bd88582428d34a70668280d0816076ccb72ea17 100644
--- a/MaterialLib/SolidModels/LinearElasticIsotropic.h
+++ b/MaterialLib/SolidModels/LinearElasticIsotropic.h
@@ -75,7 +75,7 @@ public:
 
     std::unique_ptr<
         typename MechanicsBase<DisplacementDim>::MaterialStateVariables>
-    createMaterialStateVariables() override
+    createMaterialStateVariables() const override
     {
         return std::unique_ptr<
             typename MechanicsBase<DisplacementDim>::MaterialStateVariables>{
@@ -113,14 +113,12 @@ public:
                                    DisplacementDim>::MaterialStateVariables>,
                                KelvinMatrix>>
     integrateStress(
-        double const t,
-        ProcessLib::SpatialPosition const& x,
-        double const /*dt*/,
-        KelvinVector const& eps_prev,
-        KelvinVector const& eps,
-        KelvinVector const& sigma_prev,
+        double const t, ProcessLib::SpatialPosition const& x,
+        double const /*dt*/, KelvinVector const& eps_prev,
+        KelvinVector const& eps, KelvinVector const& sigma_prev,
         typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
-            material_state_variables, double const T) override;
+            material_state_variables,
+        double const T) const override;
 
     KelvinMatrix getElasticTensor(double const t,
                                   ProcessLib::SpatialPosition const& x,
diff --git a/MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.h b/MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.h
index 8edc0986ec478dd6ead0050245d25cf608a60c7e..003c68e70ffafcbb3a7413cf94d727e95ad5e825 100644
--- a/MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.h
+++ b/MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.h
@@ -41,7 +41,7 @@ public:
 
     std::unique_ptr<
         typename MechanicsBase<DisplacementDim>::MaterialStateVariables>
-    createMaterialStateVariables() override
+    createMaterialStateVariables() const override
     {
         return LinearElasticIsotropic<
             DisplacementDim>::createMaterialStateVariables();
@@ -52,14 +52,12 @@ public:
                                    DisplacementDim>::MaterialStateVariables>,
                                KelvinMatrix>>
     integrateStress(
-        double const t,
-        ProcessLib::SpatialPosition const& x,
-        double const dt,
-        KelvinVector const& eps_prev,
-        KelvinVector const& eps,
+        double const t, ProcessLib::SpatialPosition const& x, double const dt,
+        KelvinVector const& eps_prev, KelvinVector const& eps,
         KelvinVector const& sigma_prev,
         typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
-            material_state_variables, double const T) override
+            material_state_variables,
+        double const T) const override
     {
         return LinearElasticIsotropic<DisplacementDim>::integrateStress(
             t, x, dt, eps_prev, eps, sigma_prev, material_state_variables, T);
diff --git a/MaterialLib/SolidModels/Lubby2.cpp b/MaterialLib/SolidModels/Lubby2.cpp
index ac9160c34e8b749ca81b65c8c79266798ea7fde6..ab9024fd027e381e0fbe8d41183ba4d82ee88b2e 100644
--- a/MaterialLib/SolidModels/Lubby2.cpp
+++ b/MaterialLib/SolidModels/Lubby2.cpp
@@ -73,14 +73,12 @@ boost::optional<std::tuple<typename Lubby2<DisplacementDim>::KelvinVector,
                                DisplacementDim>::MaterialStateVariables>,
                            typename Lubby2<DisplacementDim>::KelvinMatrix>>
 Lubby2<DisplacementDim>::integrateStress(
-    double const t,
-    ProcessLib::SpatialPosition const& x,
-    double const dt,
-    KelvinVector const& /*eps_prev*/,
-    KelvinVector const& eps,
+    double const t, ProcessLib::SpatialPosition const& x, double const dt,
+    KelvinVector const& /*eps_prev*/, KelvinVector const& eps,
     KelvinVector const& /*sigma_prev*/,
     typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
-        material_state_variables, double const /*T*/)
+        material_state_variables,
+    double const /*T*/) const
 {
     using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
 
@@ -204,7 +202,7 @@ void Lubby2<DisplacementDim>::calculateResidualBurgers(
     KelvinVector& strain_Max_curr,
     const KelvinVector& strain_Max_t,
     ResidualVector& res,
-    detail::LocalLubby2Properties<DisplacementDim> const& properties)
+    detail::LocalLubby2Properties<DisplacementDim> const& properties) const
 {
     // calculate stress residual
     res.template segment<KelvinVectorSize>(0).noalias() =
@@ -231,7 +229,7 @@ void Lubby2<DisplacementDim>::calculateJacobianBurgers(
     double s_eff,
     const KelvinVector& sig_i,
     const KelvinVector& eps_K_i,
-    detail::LocalLubby2Properties<DisplacementDim> const& properties)
+    detail::LocalLubby2Properties<DisplacementDim> const& properties) const
 {
     Jac.setZero();
 
diff --git a/MaterialLib/SolidModels/Lubby2.h b/MaterialLib/SolidModels/Lubby2.h
index f8b6c5ed21ec395e18eb5c93cd694a26863fcdd2..e24b37d72b62c3dbbb18cc626104995f30109f6a 100644
--- a/MaterialLib/SolidModels/Lubby2.h
+++ b/MaterialLib/SolidModels/Lubby2.h
@@ -157,7 +157,7 @@ public:
 
     std::unique_ptr<
         typename MechanicsBase<DisplacementDim>::MaterialStateVariables>
-    createMaterialStateVariables() override
+    createMaterialStateVariables() const override
     {
         return std::unique_ptr<
             typename MechanicsBase<DisplacementDim>::MaterialStateVariables>{
@@ -224,14 +224,12 @@ public:
                                    DisplacementDim>::MaterialStateVariables>,
                                KelvinMatrix>>
     integrateStress(
-        double const t,
-        ProcessLib::SpatialPosition const& x,
-        double const dt,
-        KelvinVector const& eps_prev,
-        KelvinVector const& eps,
+        double const t, ProcessLib::SpatialPosition const& x, double const dt,
+        KelvinVector const& eps_prev, KelvinVector const& eps,
         KelvinVector const& sigma_prev,
         typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
-            material_state_variables, double const T) override;
+            material_state_variables,
+        double const T) const override;
 
 private:
     /// Calculates the 18x1 residual vector.
@@ -244,7 +242,7 @@ private:
         KelvinVector& strain_Max_curr,
         const KelvinVector& strain_Max_t,
         ResidualVector& res,
-        detail::LocalLubby2Properties<DisplacementDim> const& properties);
+        detail::LocalLubby2Properties<DisplacementDim> const& properties) const;
 
     /// Calculates the 18x18 Jacobian.
     void calculateJacobianBurgers(
@@ -255,7 +253,7 @@ private:
         double s_eff,
         const KelvinVector& sig_i,
         const KelvinVector& eps_K_i,
-        detail::LocalLubby2Properties<DisplacementDim> const& properties);
+        detail::LocalLubby2Properties<DisplacementDim> const& properties) const;
 
 private:
     NumLib::NewtonRaphsonSolverParameters const _nonlinear_solver_parameters;
diff --git a/MaterialLib/SolidModels/MechanicsBase.h b/MaterialLib/SolidModels/MechanicsBase.h
index a839aa3114088e442ddca7ca7bbff52361d7d10f..26f76b6622f1597aa2fc0c059525b252bf46cb19 100644
--- a/MaterialLib/SolidModels/MechanicsBase.h
+++ b/MaterialLib/SolidModels/MechanicsBase.h
@@ -66,7 +66,7 @@ struct MechanicsBase
     /// Polymorphic creator for MaterialStateVariables objects specific for a
     /// material model.
     virtual std::unique_ptr<MaterialStateVariables>
-    createMaterialStateVariables() = 0;
+    createMaterialStateVariables() const = 0;
 
     using KelvinVector =
         MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
@@ -77,9 +77,8 @@ struct MechanicsBase
     /// constitutive relation compute function.
     /// Returns nothing in case of errors in the computation if Newton
     /// iterations did not converge, for example.
-    boost::optional<std::tuple<KelvinVector,
-                               std::unique_ptr<MaterialStateVariables>,
-                               KelvinMatrix>>
+    boost::optional<std::tuple<
+        KelvinVector, std::unique_ptr<MaterialStateVariables>, KelvinMatrix>>
     integrateStress(double const t,
                     ProcessLib::SpatialPosition const& x,
                     double const dt,
@@ -87,7 +86,7 @@ struct MechanicsBase
                     Eigen::Matrix<double, Eigen::Dynamic, 1> const& eps,
                     Eigen::Matrix<double, Eigen::Dynamic, 1> const& sigma_prev,
                     MaterialStateVariables const& material_state_variables,
-                    double const T)
+                    double const T) const
     {
         // TODO Avoid copies of data:
         // Using MatrixBase<Derived> not possible because template functions
@@ -109,9 +108,8 @@ struct MechanicsBase
     /// wrapper function.
     /// Returns nothing in case of errors in the computation if Newton
     /// iterations did not converge, for example.
-    virtual boost::optional<std::tuple<KelvinVector,
-                                       std::unique_ptr<MaterialStateVariables>,
-                                       KelvinMatrix>>
+    virtual boost::optional<std::tuple<
+        KelvinVector, std::unique_ptr<MaterialStateVariables>, KelvinMatrix>>
     integrateStress(double const t,
                     ProcessLib::SpatialPosition const& x,
                     double const dt,
@@ -119,7 +117,7 @@ struct MechanicsBase
                     KelvinVector const& eps,
                     KelvinVector const& sigma_prev,
                     MaterialStateVariables const& material_state_variables,
-                    double const T) = 0;
+                    double const T) const = 0;
 
     /// Helper type for providing access to internal variables.
     struct InternalVariable
diff --git a/MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h b/MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h
new file mode 100644
index 0000000000000000000000000000000000000000..a038e1133eab288b96ec995c027ff2e25e153bb1
--- /dev/null
+++ b/MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h
@@ -0,0 +1,65 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2018, 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 "MechanicsBase.h"
+
+namespace MaterialLib
+{
+namespace Solids
+{
+/// Choose solid material model for given element id out of a set of models,
+/// possibly using the material ids.
+///
+/// Only two possibilities yield a valid result and result in OGS_FATAL call
+/// otherwise.
+/// 1. There is only one constitutive relation in the set of
+/// constitutive_relations, then the material and element ids are ignored and
+/// the only constitutive relation (under id 0) is returned.
+/// 2. For each material id chosen for the given element id there exists a
+/// constitutive relation in the set.
+template <int DisplacementDim>
+MechanicsBase<DisplacementDim>& selectSolidConstitutiveRelation(
+    std::map<int, std::unique_ptr<MechanicsBase<DisplacementDim>>> const&
+        constitutive_relations,
+    MeshLib::PropertyVector<int> const* const material_ids,
+    std::size_t const element_id)
+{
+    int material_id;
+    if (constitutive_relations.size() == 1 || material_ids == nullptr)
+    {
+        material_id = 0;
+    }
+    else
+    {
+        material_id = (*material_ids)[element_id];
+    }
+
+    auto const constitutive_relation = constitutive_relations.find(material_id);
+    if (constitutive_relation == end(constitutive_relations))
+    {
+        OGS_FATAL(
+            "No constitutive relation found for material id %d and element %d. "
+            "There are %d constitutive relations available.",
+            material_id, element_id, constitutive_relations.size());
+    }
+
+    if (constitutive_relation->second == nullptr)
+    {
+        OGS_FATAL(
+            "The constitutive relation found for material id %d and element %d "
+            "is a nullptr, which is impossible.",
+            material_id, element_id);
+    }
+
+    return *constitutive_relation->second;
+}
+}  // namespace Solids
+}  // namespace MaterialLib
diff --git a/MeshLib/Mesh.cpp b/MeshLib/Mesh.cpp
index d3958e8f91f1251fd73ac97efb1a2b14fdd8abab..57a148fa85e95bdfcbfb55a52500742e36ee210b 100644
--- a/MeshLib/Mesh.cpp
+++ b/MeshLib/Mesh.cpp
@@ -339,6 +339,14 @@ void scaleMeshPropertyVector(MeshLib::Mesh & mesh,
         v *= factor;
 }
 
+PropertyVector<int> const* materialIDs(Mesh const& mesh)
+{
+    auto const& properties = mesh.getProperties();
+    return properties.existsPropertyVector<int>("MaterialIDs")
+               ? properties.getPropertyVector<int>("MaterialIDs")
+               : nullptr;
+}
+
 std::unique_ptr<MeshLib::Mesh> createMeshFromElementSelection(
     std::string mesh_name, std::vector<MeshLib::Element*> const& elements)
 {
diff --git a/MeshLib/Mesh.h b/MeshLib/Mesh.h
index e62ba19d36881f546fe13a1bacd8226d8e877e94..de0b088e02992cd461647a473fe129b1b2950c60 100644
--- a/MeshLib/Mesh.h
+++ b/MeshLib/Mesh.h
@@ -286,6 +286,13 @@ PropertyVector<T>* getOrCreateMeshProperty(Mesh& mesh,
     return result;
 }
 
+/// Returns the material ids property vector defined on the mesh.
+///
+/// The material ids are always an \c int property named "MaterialIDs".
+/// If the property does not exists (or is of different type), a nullptr is
+/// returned.
+PropertyVector<int> const* materialIDs(Mesh const& mesh);
+
 /// Creates a new mesh from a vector of elements.
 ///
 /// \note The elements are owned by the returned mesh object as well as the
diff --git a/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp b/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp
index a3d8e4a949a17128861411e88f4d28bc09c0164c..ad862c4a53be0f5d39c445ccd96703446433245c 100644
--- a/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp
+++ b/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp
@@ -101,9 +101,8 @@ std::unique_ptr<Process> createHydroMechanicsProcess(
             variable_p->getNumberOfComponents());
     }
 
-    // Constitutive relation.
-    auto material =
-        MaterialLib::Solids::createConstitutiveRelation<DisplacementDim>(
+    auto solid_constitutive_relations =
+        MaterialLib::Solids::createConstitutiveRelations<DisplacementDim>(
             parameters, config);
 
     // Intrinsic permeability
@@ -187,15 +186,11 @@ std::unique_ptr<Process> createHydroMechanicsProcess(
             "reference_temperature", std::numeric_limits<double>::quiet_NaN());
 
     HydroMechanicsProcessData<DisplacementDim> process_data{
-        std::move(material),
-        intrinsic_permeability,
-        specific_storage,
-        fluid_viscosity,
-        fluid_density,
-        biot_coefficient,
-        porosity,
-        solid_density,
-        specific_body_force,
+        materialIDs(mesh),      std::move(solid_constitutive_relations),
+        intrinsic_permeability, specific_storage,
+        fluid_viscosity,        fluid_density,
+        biot_coefficient,       porosity,
+        solid_density,          specific_body_force,
         reference_temperature};
 
     SecondaryVariableCollection secondary_variables;
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
index 9e6d1641f28e693e8404b5a4aa4be097fc824145..ae05e4465c717fbee315070d85e7cf1e59a3bd2e 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
@@ -13,6 +13,7 @@
 
 #include "HydroMechanicsFEM.h"
 
+#include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
 #include "MathLib/KelvinVector.h"
 #include "NumLib/Function/Interpolation.h"
 #include "ProcessLib/CoupledSolutionsForStaggeredScheme.h"
@@ -53,10 +54,14 @@ HydroMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
                           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++)
     {
-        // displacement (subscript u)
-        _ip_data.emplace_back(*_process_data.material);
+        _ip_data.emplace_back(solid_material);
         auto& ip_data = _ip_data[ip];
         auto const& sm_u = shape_matrices_u[ip];
         _ip_data[ip].integration_weight =
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
index 26dd8d5a5c761268ecf285072cfe1f48a617c5bb..f48950a68cf664d725a0123f859a316870afdfbd 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
@@ -36,7 +36,8 @@ template <typename BMatricesType, typename ShapeMatrixTypeDisplacement,
 struct IntegrationPointData final
 {
     explicit IntegrationPointData(
-        MaterialLib::Solids::MechanicsBase<DisplacementDim>& solid_material)
+        MaterialLib::Solids::MechanicsBase<DisplacementDim> const&
+            solid_material)
         : solid_material(solid_material),
           material_state_variables(
               solid_material.createMaterialStateVariables())
@@ -55,7 +56,7 @@ struct IntegrationPointData final
     typename ShapeMatricesTypePressure::NodalRowVectorType N_p;
     typename ShapeMatricesTypePressure::GlobalDimNodalMatrixType dNdx_p;
 
-    MaterialLib::Solids::MechanicsBase<DisplacementDim>& solid_material;
+    MaterialLib::Solids::MechanicsBase<DisplacementDim> const& solid_material;
     std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
         DisplacementDim>::MaterialStateVariables>
         material_state_variables;
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsProcessData.h b/ProcessLib/HydroMechanics/HydroMechanicsProcessData.h
index fea5be6470ac90e93cc580107081320ab2fafff8..88022050f9d96d4a4c1cf4b95a005f34393f329f 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsProcessData.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsProcessData.h
@@ -32,8 +32,11 @@ template <int DisplacementDim>
 struct HydroMechanicsProcessData
 {
     HydroMechanicsProcessData(
-        std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>&&
-            material_,
+        MeshLib::PropertyVector<int> const* const material_ids_,
+        std::map<int,
+                 std::unique_ptr<
+                     MaterialLib::Solids::MechanicsBase<DisplacementDim>>>&&
+            solid_materials_,
         Parameter<double> const& intrinsic_permeability_,
         Parameter<double> const& specific_storage_,
         Parameter<double> const& fluid_viscosity_,
@@ -44,7 +47,8 @@ struct HydroMechanicsProcessData
         Eigen::Matrix<double, DisplacementDim, 1>
             specific_body_force_,
         double const reference_temperature_)
-        : material{std::move(material_)},
+        : material_ids(material_ids_),
+          solid_materials{std::move(solid_materials_)},
           intrinsic_permeability(intrinsic_permeability_),
           specific_storage(specific_storage_),
           fluid_viscosity(fluid_viscosity_),
@@ -57,21 +61,7 @@ struct HydroMechanicsProcessData
     {
     }
 
-    HydroMechanicsProcessData(HydroMechanicsProcessData&& other)
-        : material{std::move(other.material)},
-          intrinsic_permeability(other.intrinsic_permeability),
-          specific_storage(other.specific_storage),
-          fluid_viscosity(other.fluid_viscosity),
-          fluid_density(other.fluid_density),
-          biot_coefficient(other.biot_coefficient),
-          porosity(other.porosity),
-          solid_density(other.solid_density),
-          specific_body_force(other.specific_body_force),
-          dt(other.dt),
-          t(other.t),
-          reference_temperature(other.reference_temperature)
-    {
-    }
+    HydroMechanicsProcessData(HydroMechanicsProcessData&& other) = default;
 
     //! Copies are forbidden.
     HydroMechanicsProcessData(HydroMechanicsProcessData const&) = delete;
@@ -82,10 +72,14 @@ struct HydroMechanicsProcessData
     //! Assignments are not needed.
     void operator=(HydroMechanicsProcessData&&) = 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::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>
-        material;
+    std::map<
+        int,
+        std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>>
+        solid_materials;
     /// Permeability of the solid. A scalar quantity, Parameter<double>.
     Parameter<double> const& intrinsic_permeability;
     /// Volumetric average specific storage of the solid and fluid phases.
diff --git a/ProcessLib/LIE/HydroMechanics/CreateHydroMechanicsProcess.cpp b/ProcessLib/LIE/HydroMechanics/CreateHydroMechanicsProcess.cpp
index bac9e7073808a4c1fe941ecb828ecefd75b3610c..67c8e20bb186f2863f9e48265f36dcfb07d43f89 100644
--- a/ProcessLib/LIE/HydroMechanics/CreateHydroMechanicsProcess.cpp
+++ b/ProcessLib/LIE/HydroMechanics/CreateHydroMechanicsProcess.cpp
@@ -123,11 +123,9 @@ std::unique_ptr<Process> createHydroMechanicsProcess(
     else
         process_variables.push_back(std::move(p_u_process_variables));
 
-
-    // Constitutive relation.
-    auto material =
-        MaterialLib::Solids::createConstitutiveRelation<GlobalDim>(
-            parameters, config);
+    auto solid_constitutive_relations =
+        MaterialLib::Solids::createConstitutiveRelations<GlobalDim>(parameters,
+                                                                    config);
 
     // Intrinsic permeability
     auto& intrinsic_permeability = findParameter<double>(
@@ -309,7 +307,8 @@ std::unique_ptr<Process> createHydroMechanicsProcess(
             "reference_temperature", std::numeric_limits<double>::quiet_NaN());
 
     HydroMechanicsProcessData<GlobalDim> process_data{
-        std::move(material),
+        materialIDs(mesh),
+        std::move(solid_constitutive_relations),
         intrinsic_permeability,
         specific_storage,
         fluid_viscosity,
diff --git a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcessData.h b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcessData.h
index f19c93b24d4e7c788dde3931a3e9eb4d400e7379..7d86a1e93c2b425734054d6cf4d4fa09329b7055 100644
--- a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcessData.h
+++ b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcessData.h
@@ -36,8 +36,11 @@ template <unsigned GlobalDim>
 struct HydroMechanicsProcessData
 {
     HydroMechanicsProcessData(
-        std::unique_ptr<MaterialLib::Solids::MechanicsBase<GlobalDim>>&&
-            material_,
+        MeshLib::PropertyVector<int> const* const material_ids_,
+        std::map<
+            int,
+            std::unique_ptr<MaterialLib::Solids::MechanicsBase<GlobalDim>>>&&
+            solid_materials_,
         Parameter<double> const& intrinsic_permeability_,
         Parameter<double> const& specific_storage_,
         Parameter<double> const& fluid_viscosity_,
@@ -54,7 +57,8 @@ struct HydroMechanicsProcessData
         Parameter<double> const& initial_fracture_effective_stress_,
         bool const deactivate_matrix_in_flow_,
         double const reference_temperature_)
-        : material{std::move(material_)},
+        : material_ids(material_ids_),
+          solid_materials{std::move(solid_materials_)},
           intrinsic_permeability(intrinsic_permeability_),
           specific_storage(specific_storage_),
           fluid_viscosity(fluid_viscosity_),
@@ -72,29 +76,7 @@ struct HydroMechanicsProcessData
     {
     }
 
-    HydroMechanicsProcessData(HydroMechanicsProcessData&& other)
-        : material{std::move(other.material)},
-          intrinsic_permeability(other.intrinsic_permeability),
-          specific_storage(other.specific_storage),
-          fluid_viscosity(other.fluid_viscosity),
-          fluid_density(other.fluid_density),
-          biot_coefficient(other.biot_coefficient),
-          porosity(other.porosity),
-          solid_density(other.solid_density),
-          specific_body_force(other.specific_body_force),
-          fracture_model{std::move(other.fracture_model)},
-          fracture_property{std::move(other.fracture_property)},
-          initial_effective_stress(other.initial_effective_stress),
-          initial_fracture_effective_stress(
-              other.initial_fracture_effective_stress),
-          deactivate_matrix_in_flow(other.deactivate_matrix_in_flow),
-          p_element_status(std::move(other.p_element_status)),
-          p0(other.p0),
-          dt(other.dt),
-          t(other.t),
-          reference_temperature(other.reference_temperature)
-    {
-    }
+    HydroMechanicsProcessData(HydroMechanicsProcessData&& other) = default;
 
     //! Copies are forbidden.
     HydroMechanicsProcessData(HydroMechanicsProcessData const&) = delete;
@@ -105,7 +87,10 @@ struct HydroMechanicsProcessData
     //! Assignments are not needed.
     void operator=(HydroMechanicsProcessData&&) = delete;
 
-    std::unique_ptr<MaterialLib::Solids::MechanicsBase<GlobalDim>> material;
+    MeshLib::PropertyVector<int> const* const material_ids;
+    std::map<int,
+             std::unique_ptr<MaterialLib::Solids::MechanicsBase<GlobalDim>>>
+        solid_materials;
     Parameter<double> const& intrinsic_permeability;
     Parameter<double> const& specific_storage;
     Parameter<double> const& fluid_viscosity;
diff --git a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix-impl.h b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix-impl.h
index 1ff35e41d21e86456ad5e1d9bfdc3e4c1c59da80..5d6bcd9e6e4398bb2a69f035dbfd6bc79ef46caa 100644
--- a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix-impl.h
+++ b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix-impl.h
@@ -11,8 +11,9 @@
 
 #include "HydroMechanicsLocalAssemblerMatrix.h"
 
-#include "MathLib/KelvinVector.h"
 #include "MaterialLib/PhysicalConstant.h"
+#include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
+#include "MathLib/KelvinVector.h"
 #include "MeshLib/ElementStatus.h"
 #include "NumLib/Fem/CoordinatesMapping/NaturalNodeCoordinates.h"
 #include "ProcessLib/Deformation/LinearBMatrix.h"
@@ -61,13 +62,16 @@ HydroMechanicsLocalAssemblerMatrix<ShapeFunctionDisplacement,
                           IntegrationMethod, GlobalDim>(e, is_axially_symmetric,
                                                         integration_method);
 
+    auto& solid_material = MaterialLib::Solids::selectSolidConstitutiveRelation(
+        _process_data.solid_materials, _process_data.material_ids, e.getID());
+
     SpatialPosition x_position;
     x_position.setElementID(e.getID());
     for (unsigned ip = 0; ip < n_integration_points; ip++)
     {
         x_position.setIntegrationPoint(ip);
 
-        _ip_data.emplace_back(*_process_data.material);
+        _ip_data.emplace_back(solid_material);
         auto& ip_data = _ip_data[ip];
         auto const& sm_u = shape_matrices_u[ip];
         auto const& sm_p = shape_matrices_p[ip];
diff --git a/ProcessLib/LIE/SmallDeformation/CreateSmallDeformationProcess.cpp b/ProcessLib/LIE/SmallDeformation/CreateSmallDeformationProcess.cpp
index edcac74b40af1743435018ad89a835b2e3b6c011..30cf7a9bf4e02765f056f373e68093fe5b951d71 100644
--- a/ProcessLib/LIE/SmallDeformation/CreateSmallDeformationProcess.cpp
+++ b/ProcessLib/LIE/SmallDeformation/CreateSmallDeformationProcess.cpp
@@ -98,9 +98,8 @@ std::unique_ptr<Process> createSmallDeformationProcess(
         process_variables;
     process_variables.push_back(std::move(per_process_variables));
 
-    // Constitutive relation.
-    auto material =
-        MaterialLib::Solids::createConstitutiveRelation<DisplacementDim>(
+    auto solid_constitutive_relations =
+        MaterialLib::Solids::createConstitutiveRelations<DisplacementDim>(
             parameters, config);
 
     // Fracture constitutive relation.
@@ -174,8 +173,9 @@ std::unique_ptr<Process> createSmallDeformationProcess(
             "reference_temperature", std::numeric_limits<double>::quiet_NaN());
 
     SmallDeformationProcessData<DisplacementDim> process_data(
-        std::move(material), std::move(fracture_model),
-        std::move(vec_fracture_property), reference_temperature);
+        materialIDs(mesh), std::move(solid_constitutive_relations),
+        std::move(fracture_model), std::move(vec_fracture_property),
+        reference_temperature);
 
     SecondaryVariableCollection secondary_variables;
 
diff --git a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrix-impl.h b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrix-impl.h
index 8229c464b94b22a21433d56ddf3d3f8876bd2166..1b701ab91d31e1906966f4259d2ebf1113f0f4ea 100644
--- a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrix-impl.h
+++ b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrix-impl.h
@@ -17,6 +17,7 @@
 #include <Eigen/Eigen>
 
 #include "MaterialLib/PhysicalConstant.h"
+#include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
@@ -62,9 +63,12 @@ SmallDeformationLocalAssemblerMatrix<ShapeFunction, IntegrationMethod,
                           DisplacementDim>(e, is_axially_symmetric,
                                            _integration_method);
 
+    auto& 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(*_process_data._material);
+        _ip_data.emplace_back(solid_material);
         auto& ip_data = _ip_data[ip];
         auto const& sm = shape_matrices[ip];
         ip_data.N = sm.N;
diff --git a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture-impl.h b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture-impl.h
index 50c6b3c58abb1f8339ef5f9f655930a1f7976fee..341d7d2ab9e6b289a4f7b4f8da321eb38c3f5f8c 100644
--- a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture-impl.h
+++ b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture-impl.h
@@ -81,9 +81,12 @@ SmallDeformationLocalAssemblerMatrixNearFracture<ShapeFunction,
     _ip_data.reserve(n_integration_points);
     _secondary_data.N.resize(n_integration_points);
 
+    auto& 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(*_process_data._material);
+        _ip_data.emplace_back(solid_material);
         auto& ip_data = _ip_data[ip];
         auto const& sm = shape_matrices[ip];
         ip_data.N = sm.N;
diff --git a/ProcessLib/LIE/SmallDeformation/SmallDeformationProcessData.h b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcessData.h
index cba9b7f93e2f689deab79ba6cfd5031b65fc7cdb..9d6dc28fbf706d77be8d199b8f390234adb3ef2c 100644
--- a/ProcessLib/LIE/SmallDeformation/SmallDeformationProcessData.h
+++ b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcessData.h
@@ -33,27 +33,25 @@ template <int DisplacementDim>
 struct SmallDeformationProcessData
 {
     SmallDeformationProcessData(
-        std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>&&
-            material,
+        MeshLib::PropertyVector<int> const* const material_ids_,
+        std::map<int,
+                 std::unique_ptr<
+                     MaterialLib::Solids::MechanicsBase<DisplacementDim>>>&&
+            solid_materials_,
         std::unique_ptr<
             MaterialLib::Fracture::FractureModelBase<DisplacementDim>>&&
             fracture_model,
         std::vector<std::unique_ptr<FractureProperty>>&& vec_fracture_prop,
         double const reference_temperature)
-        : _material{std::move(material)},
+        : material_ids(material_ids_),
+          solid_materials{std::move(solid_materials_)},
           _fracture_model{std::move(fracture_model)},
           _vec_fracture_property(std::move(vec_fracture_prop)),
           _reference_temperature(reference_temperature)
     {
     }
 
-    SmallDeformationProcessData(SmallDeformationProcessData&& other)
-        : _material{std::move(other._material)},
-          _fracture_model{std::move(other._fracture_model)},
-          _vec_fracture_property(std::move(other._vec_fracture_property)),
-          _reference_temperature(other._reference_temperature)
-    {
-    }
+    SmallDeformationProcessData(SmallDeformationProcessData&& other) = default;
 
     //! Copies are forbidden.
     SmallDeformationProcessData(SmallDeformationProcessData const&) = delete;
@@ -64,8 +62,14 @@ struct SmallDeformationProcessData
     //! Assignments are not needed.
     void operator=(SmallDeformationProcessData&&) = delete;
 
-    std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>
-        _material;
+    MeshLib::PropertyVector<int> const* const material_ids;
+
+    /// The constitutive relation for the mechanical part.
+    std::map<
+        int,
+        std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>>
+        solid_materials;
+
     std::unique_ptr<MaterialLib::Fracture::FractureModelBase<DisplacementDim>>
         _fracture_model;
     std::vector<std::unique_ptr<FractureProperty>> _vec_fracture_property;
diff --git a/ProcessLib/PhaseField/PhaseFieldFEM.h b/ProcessLib/PhaseField/PhaseFieldFEM.h
index e4310b752473b3dd6f3fcd147a853351d432f3af..d8baae9231332189f29424d86ceb330afc52385d 100644
--- a/ProcessLib/PhaseField/PhaseFieldFEM.h
+++ b/ProcessLib/PhaseField/PhaseFieldFEM.h
@@ -32,7 +32,8 @@ template <typename BMatricesType, typename ShapeMatrixType, int DisplacementDim>
 struct IntegrationPointData final
 {
     explicit IntegrationPointData(
-        MaterialLib::Solids::MechanicsBase<DisplacementDim>& solid_material)
+        MaterialLib::Solids::MechanicsBase<DisplacementDim> const&
+            solid_material)
         : solid_material(solid_material),
           material_state_variables(
               solid_material.createMaterialStateVariables())
@@ -48,7 +49,7 @@ struct IntegrationPointData final
         sigma;
     double strain_energy_tensile, elastic_energy;
 
-    MaterialLib::Solids::MechanicsBase<DisplacementDim>& solid_material;
+    MaterialLib::Solids::MechanicsBase<DisplacementDim> const& solid_material;
     std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
         DisplacementDim>::MaterialStateVariables>
         material_state_variables;
@@ -72,7 +73,8 @@ struct IntegrationPointData final
                                     DisplacementVectorType const& /*u*/,
                                     double const degradation)
     {
-        static_cast<MaterialLib::Solids::PhaseFieldExtension<DisplacementDim>&>(
+        static_cast<
+            MaterialLib::Solids::PhaseFieldExtension<DisplacementDim> const&>(
             solid_material)
             .calculateDegradedStress(t, x_position, eps, strain_energy_tensile,
                                      sigma_tensile, sigma_compressive,
diff --git a/ProcessLib/RichardsFlow/CreateRichardsFlowMaterialProperties.cpp b/ProcessLib/RichardsFlow/CreateRichardsFlowMaterialProperties.cpp
index e6099f0c85f68a86e82429955104ecae6e87292e..26b22ca0d965747c9bbf24cb317d88f9604c7ee3 100644
--- a/ProcessLib/RichardsFlow/CreateRichardsFlowMaterialProperties.cpp
+++ b/ProcessLib/RichardsFlow/CreateRichardsFlowMaterialProperties.cpp
@@ -34,8 +34,7 @@ namespace RichardsFlow
 std::unique_ptr<RichardsFlowMaterialProperties>
 createRichardsFlowMaterialProperties(
     BaseLib::ConfigTree const& config,
-    boost::optional<MeshLib::PropertyVector<int> const&>
-        material_ids,
+    MeshLib::PropertyVector<int> const* const material_ids,
     std::vector<std::unique_ptr<ParameterBase>> const& parameters)
 {
     DBUG("Reading material properties of Richards flow process.");
diff --git a/ProcessLib/RichardsFlow/CreateRichardsFlowMaterialProperties.h b/ProcessLib/RichardsFlow/CreateRichardsFlowMaterialProperties.h
index f45efa6942870884405506a12a8685ac5c11ab06..193ed6b9cd91347af34ba72418a14ba547567e15 100644
--- a/ProcessLib/RichardsFlow/CreateRichardsFlowMaterialProperties.h
+++ b/ProcessLib/RichardsFlow/CreateRichardsFlowMaterialProperties.h
@@ -22,8 +22,7 @@ namespace RichardsFlow
 std::unique_ptr<RichardsFlowMaterialProperties>
 createRichardsFlowMaterialProperties(
     BaseLib::ConfigTree const& config,
-    boost::optional<MeshLib::PropertyVector<int> const&>
-        material_ids,
+    MeshLib::PropertyVector<int> const* material_ids,
     std::vector<std::unique_ptr<ParameterBase>> const& parameters);
 
 }  // end namespace
diff --git a/ProcessLib/RichardsFlow/CreateRichardsFlowProcess.cpp b/ProcessLib/RichardsFlow/CreateRichardsFlowProcess.cpp
index 007175036af7a11ce18ee47718e7723ee5f8d925..f0651c3825b8d45236815889cd977936a64e0279 100644
--- a/ProcessLib/RichardsFlow/CreateRichardsFlowProcess.cpp
+++ b/ProcessLib/RichardsFlow/CreateRichardsFlowProcess.cpp
@@ -80,12 +80,11 @@ std::unique_ptr<Process> createRichardsFlowProcess(
     //! \ogs_file_param{prj__processes__process__RICHARDS_FLOW__material_property}
     auto const& mat_config = config.getConfigSubtree("material_property");
 
-    boost::optional<MeshLib::PropertyVector<int> const&> material_ids;
-    if (mesh.getProperties().existsPropertyVector<int>("MaterialIDs"))
+    auto const material_ids =
+        mesh.getProperties().getPropertyVector<int>("MaterialIDs");
+    if (material_ids != nullptr)
     {
         INFO("The Richards flow is in heterogeneous porous media.");
-        material_ids =
-            *mesh.getProperties().getPropertyVector<int>("MaterialIDs");
     }
     else
     {
diff --git a/ProcessLib/RichardsFlow/RichardsFlowMaterialProperties.cpp b/ProcessLib/RichardsFlow/RichardsFlowMaterialProperties.cpp
index c1efa6e2b6d87fb5da59205bb36efa9e299c28bb..ebbf4186543d3e4642b1b38bb53974871c055654 100644
--- a/ProcessLib/RichardsFlow/RichardsFlowMaterialProperties.cpp
+++ b/ProcessLib/RichardsFlow/RichardsFlowMaterialProperties.cpp
@@ -29,7 +29,7 @@ namespace ProcessLib
 namespace RichardsFlow
 {
 RichardsFlowMaterialProperties::RichardsFlowMaterialProperties(
-    boost::optional<MeshLib::PropertyVector<int> const&> const material_ids,
+    MeshLib::PropertyVector<int> const* const material_ids,
     std::unique_ptr<MaterialLib::Fluid::FluidProperties>&& fluid_properties,
     std::vector<std::unique_ptr<MaterialLib::PorousMedium::Permeability>>&&
         intrinsic_permeability_models,
diff --git a/ProcessLib/RichardsFlow/RichardsFlowMaterialProperties.h b/ProcessLib/RichardsFlow/RichardsFlowMaterialProperties.h
index 3224975029d53bff873040374c6fd487806a543e..5711e3fcc07556f150c3c6904d73d52e3aadfedc 100644
--- a/ProcessLib/RichardsFlow/RichardsFlowMaterialProperties.h
+++ b/ProcessLib/RichardsFlow/RichardsFlowMaterialProperties.h
@@ -49,7 +49,7 @@ public:
     using ArrayType = MaterialLib::Fluid::FluidProperty::ArrayType;
 
     RichardsFlowMaterialProperties(
-        boost::optional<MeshLib::PropertyVector<int> const&> const material_ids,
+        MeshLib::PropertyVector<int> const* const material_ids,
         std::unique_ptr<MaterialLib::Fluid::FluidProperties>&& fluid_properties,
         std::vector<std::unique_ptr<MaterialLib::PorousMedium::Permeability>>&&
             intrinsic_permeability_models,
@@ -108,7 +108,7 @@ private:
     /**
     *  Material IDs must be given as mesh element properties.
     */
-    boost::optional<MeshLib::PropertyVector<int> const&> const _material_ids;
+    MeshLib::PropertyVector<int> const* const _material_ids;
 
     const std::unique_ptr<MaterialLib::Fluid::FluidProperties>
         _fluid_properties;
diff --git a/ProcessLib/RichardsMechanics/CreateRichardsMechanicsProcess.cpp b/ProcessLib/RichardsMechanics/CreateRichardsMechanicsProcess.cpp
index d901fb3f40a97c73bad8efd10a0ed9caff65e009..7234773c662a7cf2d2fd12961b59f70fc032d59a 100644
--- a/ProcessLib/RichardsMechanics/CreateRichardsMechanicsProcess.cpp
+++ b/ProcessLib/RichardsMechanics/CreateRichardsMechanicsProcess.cpp
@@ -101,9 +101,8 @@ std::unique_ptr<Process> createRichardsMechanicsProcess(
             variable_p->getNumberOfComponents());
     }
 
-    // Constitutive relation.
-    auto solid_material =
-        MaterialLib::Solids::createConstitutiveRelation<DisplacementDim>(
+    auto solid_constitutive_relations =
+        MaterialLib::Solids::createConstitutiveRelations<DisplacementDim>(
             parameters, config);
 
     // Intrinsic permeability
@@ -174,14 +173,15 @@ std::unique_ptr<Process> createRichardsMechanicsProcess(
         //! \ogs_file_param{prj__processes__process__RICHARDS_MECHANICS__material_property}
         config.getConfigSubtree("material_property");
 
-    boost::optional<MeshLib::PropertyVector<int> const&> material_ids;
-    if (mesh.getProperties().existsPropertyVector<int>("MaterialIDs"))
+    auto const material_ids =
+        mesh.getProperties().existsPropertyVector<int>("MaterialIDs")
+            ? mesh.getProperties().getPropertyVector<int>("MaterialIDs")
+            : nullptr;
+    if (material_ids != nullptr)
     {
         INFO(
             "MaterialIDs vector found; the Richards flow is in heterogeneous "
             "porous media.");
-        material_ids =
-            *mesh.getProperties().getPropertyVector<int>("MaterialIDs");
     }
     else
     {
@@ -194,8 +194,9 @@ std::unique_ptr<Process> createRichardsMechanicsProcess(
             flow_material_config, material_ids, parameters);
 
     RichardsMechanicsProcessData<DisplacementDim> process_data{
+        material_ids,
         std::move(flow_material),
-        std::move(solid_material),
+        std::move(solid_constitutive_relations),
         intrinsic_permeability,
         fluid_bulk_modulus,
         biot_coefficient,
diff --git a/ProcessLib/RichardsMechanics/IntegrationPointData.h b/ProcessLib/RichardsMechanics/IntegrationPointData.h
index cc03455f71169d0e915a049cbb58cb84ba0c84ea..2cc15863b219e9da6825661988977557eb4f2f71 100644
--- a/ProcessLib/RichardsMechanics/IntegrationPointData.h
+++ b/ProcessLib/RichardsMechanics/IntegrationPointData.h
@@ -24,7 +24,8 @@ template <typename BMatricesType, typename ShapeMatrixTypeDisplacement,
 struct IntegrationPointData final
 {
     explicit IntegrationPointData(
-        MaterialLib::Solids::MechanicsBase<DisplacementDim>& solid_material)
+        MaterialLib::Solids::MechanicsBase<DisplacementDim> const&
+            solid_material)
         : solid_material(solid_material),
           material_state_variables(
               solid_material.createMaterialStateVariables())
@@ -55,7 +56,7 @@ struct IntegrationPointData final
 
     double saturation;
 
-    MaterialLib::Solids::MechanicsBase<DisplacementDim>& solid_material;
+    MaterialLib::Solids::MechanicsBase<DisplacementDim> const& solid_material;
     std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
         DisplacementDim>::MaterialStateVariables>
         material_state_variables;
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
index 5d6f227360b9d98cca813a76856fe7e92175ee9f..c7d2fb5f04e0667c5a2b4ca9b7bea8c60b2c5eae 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
@@ -13,6 +13,7 @@
 
 #include "RichardsMechanicsFEM.h"
 
+#include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
 #include "MathLib/KelvinVector.h"
 #include "NumLib/Function/Interpolation.h"
 #include "ProcessLib/CoupledSolutionsForStaggeredScheme.h"
@@ -54,10 +55,14 @@ RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
                           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++)
     {
-        // displacement (subscript u)
-        _ip_data.emplace_back(*_process_data.solid_material);
+        _ip_data.emplace_back(solid_material);
         auto& ip_data = _ip_data[ip];
         auto const& sm_u = shape_matrices_u[ip];
         _ip_data[ip].integration_weight =
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h b/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h
index df5a6edf8d02bc2ef5e3bca2b89d08ca5ea9a624..f737f7ce28ada2b4bbc29dcfc7e26832ddf99a16 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h
@@ -34,11 +34,14 @@ template <int DisplacementDim>
 struct RichardsMechanicsProcessData
 {
     RichardsMechanicsProcessData(
+        MeshLib::PropertyVector<int> const* const material_ids_,
         std::unique_ptr<
             ProcessLib::RichardsFlow::RichardsFlowMaterialProperties>&&
             flow_material_,
-        std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>&&
-            solid_material_,
+        std::map<int,
+                 std::unique_ptr<
+                     MaterialLib::Solids::MechanicsBase<DisplacementDim>>>&&
+            solid_materials_,
         Parameter<double> const& intrinsic_permeability_,
         Parameter<double> const& fluid_bulk_modulus_,
         Parameter<double> const& biot_coefficient_,
@@ -47,8 +50,9 @@ struct RichardsMechanicsProcessData
         Parameter<double> const& temperature_,
         Eigen::Matrix<double, DisplacementDim, 1>
             specific_body_force_)
-        : flow_material{std::move(flow_material_)},
-          solid_material{std::move(solid_material_)},
+        : material_ids(material_ids_),
+          flow_material{std::move(flow_material_)},
+          solid_materials{std::move(solid_materials_)},
           intrinsic_permeability(intrinsic_permeability_),
           fluid_bulk_modulus(fluid_bulk_modulus_),
           biot_coefficient(biot_coefficient_),
@@ -69,12 +73,16 @@ struct RichardsMechanicsProcessData
     void operator=(RichardsMechanicsProcessData const&) = delete;
     void operator=(RichardsMechanicsProcessData&&) = delete;
 
+    MeshLib::PropertyVector<int> const* const material_ids;
+
     std::unique_ptr<ProcessLib::RichardsFlow::RichardsFlowMaterialProperties>
         flow_material;
 
     /// The constitutive relation for the mechanical part.
-    std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>
-        solid_material;
+    std::map<
+        int,
+        std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>>
+        solid_materials;
     /// Permeability of the solid. A scalar quantity, Parameter<double>.
     Parameter<double> const& intrinsic_permeability;
     /// Fluid's bulk modulus. A scalar quantity, Parameter<double>.
diff --git a/ProcessLib/SmallDeformation/CreateSmallDeformationProcess.cpp b/ProcessLib/SmallDeformation/CreateSmallDeformationProcess.cpp
index 15dee6ba84edde6ff6aef94d3d86c00007f46abd..ddb267a395797e29aad419c995b7958534728f35 100644
--- a/ProcessLib/SmallDeformation/CreateSmallDeformationProcess.cpp
+++ b/ProcessLib/SmallDeformation/CreateSmallDeformationProcess.cpp
@@ -64,9 +64,8 @@ createSmallDeformationProcess(
         process_variables;
     process_variables.push_back(std::move(per_process_variables));
 
-    // Constitutive relation.
-    auto material =
-        MaterialLib::Solids::createConstitutiveRelation<DisplacementDim>(
+    auto solid_constitutive_relations =
+        MaterialLib::Solids::createConstitutiveRelations<DisplacementDim>(
             parameters, config);
 
     // Solid density
@@ -100,8 +99,8 @@ createSmallDeformationProcess(
             "reference_temperature", std::numeric_limits<double>::quiet_NaN());
 
     SmallDeformationProcessData<DisplacementDim> process_data{
-        std::move(material), solid_density, specific_body_force,
-        reference_temperature};
+        materialIDs(mesh), std::move(solid_constitutive_relations),
+        solid_density, specific_body_force, reference_temperature};
 
     SecondaryVariableCollection secondary_variables;
 
diff --git a/ProcessLib/SmallDeformation/SmallDeformationFEM.h b/ProcessLib/SmallDeformation/SmallDeformationFEM.h
index d304a54cfbcd76de996d1e4d938f8c524ded7617..f15ad75ede39695bcd23a81cc616517503636d7e 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationFEM.h
+++ b/ProcessLib/SmallDeformation/SmallDeformationFEM.h
@@ -13,7 +13,7 @@
 #include <vector>
 
 #include "MaterialLib/PhysicalConstant.h"
-#include "MaterialLib/SolidModels/LinearElasticIsotropic.h"
+#include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "NumLib/Extrapolation/ExtrapolatableElement.h"
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
@@ -38,7 +38,8 @@ template <typename BMatricesType, typename ShapeMatricesType,
 struct IntegrationPointData final
 {
     explicit IntegrationPointData(
-        MaterialLib::Solids::MechanicsBase<DisplacementDim>& solid_material)
+        MaterialLib::Solids::MechanicsBase<DisplacementDim> const&
+            solid_material)
         : solid_material(solid_material),
           material_state_variables(
               solid_material.createMaterialStateVariables())
@@ -49,7 +50,7 @@ struct IntegrationPointData final
     typename BMatricesType::KelvinVectorType eps, eps_prev;
     double free_energy_density = 0;
 
-    MaterialLib::Solids::MechanicsBase<DisplacementDim>& solid_material;
+    MaterialLib::Solids::MechanicsBase<DisplacementDim> const& solid_material;
     std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
         DisplacementDim>::MaterialStateVariables>
         material_state_variables;
@@ -125,9 +126,15 @@ public:
                               IntegrationMethod, DisplacementDim>(
                 e, is_axially_symmetric, _integration_method);
 
+        auto& 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(*_process_data.material);
+            _ip_data.emplace_back(solid_material);
             auto& ip_data = _ip_data[ip];
             auto const& sm = shape_matrices[ip];
             _ip_data[ip].integration_weight =
diff --git a/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp b/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp
index 9b4e112c4de7d827e422ccb24875006d8d99d6cb..20cb7734353c852f7e48a3461512ebedf994f0f6 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp
+++ b/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp
@@ -112,9 +112,24 @@ void SmallDeformationProcess<DisplacementDim>::initializeConcreteProcess(
                          getExtrapolator(), _local_assemblers,
                          &LocalAssemblerInterface::getIntPtEpsilon));
 
+    //
     // enable output of internal variables defined by material models
-    auto const internal_variables =
-        _process_data.material->getInternalVariables();
+    //
+
+    // Collect the internal variables for all solid materials.
+    std::vector<typename MaterialLib::Solids::MechanicsBase<
+        DisplacementDim>::InternalVariable>
+        internal_variables;
+    for (auto const& material_id__solid_material :
+         _process_data.solid_materials)
+    {
+        auto const variables =
+            material_id__solid_material.second->getInternalVariables();
+        copy(begin(variables), end(variables),
+             back_inserter(internal_variables));
+    }
+
+    // Register the internal variables.
     for (auto const& internal_variable : internal_variables)
     {
         auto const& name = internal_variable.name;
diff --git a/ProcessLib/SmallDeformation/SmallDeformationProcessData.h b/ProcessLib/SmallDeformation/SmallDeformationProcessData.h
index 29a24c8e4262840ff26bdc95ab805f0d01c92487..751e6754abe94002940b9ad1db17ea95b4271a33 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationProcessData.h
+++ b/ProcessLib/SmallDeformation/SmallDeformationProcessData.h
@@ -32,28 +32,24 @@ template <int DisplacementDim>
 struct SmallDeformationProcessData
 {
     SmallDeformationProcessData(
-        std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>&&
-            material,
+        MeshLib::PropertyVector<int> const* const material_ids_,
+        std::map<int,
+                 std::unique_ptr<
+                     MaterialLib::Solids::MechanicsBase<DisplacementDim>>>&&
+            solid_materials_,
         Parameter<double> const& solid_density_,
         Eigen::Matrix<double, DisplacementDim, 1>
             specific_body_force_,
         double const reference_temperature_)
-        : material{std::move(material)},
+        : material_ids(material_ids_),
+          solid_materials{std::move(solid_materials_)},
           solid_density(solid_density_),
           specific_body_force(std::move(specific_body_force_)),
           reference_temperature(reference_temperature_)
     {
     }
 
-    SmallDeformationProcessData(SmallDeformationProcessData&& other)
-        : material{std::move(other.material)},
-          solid_density(other.solid_density),
-          specific_body_force(other.specific_body_force),
-          dt{other.dt},
-          t{other.t},
-          reference_temperature(other.reference_temperature)
-    {
-    }
+    SmallDeformationProcessData(SmallDeformationProcessData&& other) = default;
 
     //! Copies are forbidden.
     SmallDeformationProcessData(SmallDeformationProcessData const&) = delete;
@@ -64,8 +60,12 @@ struct SmallDeformationProcessData
     //! Assignments are not needed.
     void operator=(SmallDeformationProcessData&&) = delete;
 
-    std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>
-        material;
+    MeshLib::PropertyVector<int> const* const material_ids;
+
+    std::map<
+        int,
+        std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>>
+        solid_materials;
     /// Solid's density. A scalar quantity, Parameter<double>.
     Parameter<double> const& solid_density;
     /// Specific body forces applied to the solid.
diff --git a/ProcessLib/SmallDeformation/Tests.cmake b/ProcessLib/SmallDeformation/Tests.cmake
index 1787f24f1e42cd063841776e3d7a8858c2ab471a..12a969336368ab9d800e2ba66d71e0b2e0b4cf50 100644
--- a/ProcessLib/SmallDeformation/Tests.cmake
+++ b/ProcessLib/SmallDeformation/Tests.cmake
@@ -287,9 +287,9 @@ AddTest(
     ref_axisymmetric_sphere_pcs_0_ts_100_t_1.000000.vtu axisymmetric_sphere_pcs_0_ts_100_t_1.000000.vtu sigma sigma 1e-15 0
 )
 
-# Two materials in gravity field
+# Two materials in gravity field variation of density
 AddTest(
-    NAME SmallDeformation_SDL_two_material_gravity
+    NAME SmallDeformation_SDL_two_material_gravity_density
     PATH Mechanics/Linear
     EXECUTABLE_ARGS two_material_gravity.prj
     TESTER vtkdiff
@@ -300,6 +300,18 @@ AddTest(
     expected_two_material_gravity.vtu two_material_gravity_pcs_0_ts_1_t_1.000000.vtu epsilon epsilon 5e-14 0
 )
 
+# Two materials in gravity field variation of modulus of elasticity
+AddTest(
+    NAME SmallDeformation_SDL_two_material_gravity_Emodulus
+    PATH Mechanics/Linear
+    EXECUTABLE_ARGS two_material_gravity_Emodulus.prj
+    TESTER vtkdiff
+    REQUIREMENTS NOT OGS_USE_MPI
+    DIFF_DATA
+    expected_two_material_gravity_Emodulus.vtu two_material_gravity_Emodulus_pcs_0_ts_1_t_1.000000.vtu displacement displacement 5e-14 0
+    expected_two_material_gravity_Emodulus.vtu two_material_gravity_Emodulus_pcs_0_ts_1_t_1.000000.vtu sigma sigma 1e-13 0
+    expected_two_material_gravity_Emodulus.vtu two_material_gravity_Emodulus_pcs_0_ts_1_t_1.000000.vtu epsilon epsilon 5e-14 0
+)
 AddTest(
     NAME PythonBCSmallDeformationPiston
     PATH Mechanics/Linear/PythonPiston
diff --git a/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM.h b/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM.h
index 116e093ce817dd1ec946155f7f2ebe645ee1e258..d7741e0e8abcde82671d2b55392df10fe9cf0ec7 100644
--- a/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM.h
+++ b/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM.h
@@ -32,7 +32,8 @@ template <typename BMatricesType, typename ShapeMatrixType, int DisplacementDim>
 struct IntegrationPointData final
 {
     explicit IntegrationPointData(
-        MaterialLib::Solids::MechanicsBase<DisplacementDim>& solid_material)
+        MaterialLib::Solids::MechanicsBase<DisplacementDim> const&
+            solid_material)
         : solid_material(solid_material),
           material_state_variables(
               solid_material.createMaterialStateVariables())
@@ -50,7 +51,7 @@ struct IntegrationPointData final
     double strain_energy_tensile, elastic_energy;
     typename ShapeMatrixType::GlobalDimVectorType heatflux;
 
-    MaterialLib::Solids::MechanicsBase<DisplacementDim>& solid_material;
+    MaterialLib::Solids::MechanicsBase<DisplacementDim> const& solid_material;
     std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
         DisplacementDim>::MaterialStateVariables>
         material_state_variables;
@@ -78,7 +79,8 @@ struct IntegrationPointData final
     {
         eps_m.noalias() = eps - alpha * delta_T * Invariants::identity2;
 
-        static_cast<MaterialLib::Solids::PhaseFieldExtension<DisplacementDim>&>(
+        static_cast<
+            MaterialLib::Solids::PhaseFieldExtension<DisplacementDim> const&>(
             solid_material)
             .calculateDegradedStress(
                 t, x_position, eps_m, strain_energy_tensile, sigma_tensile,
diff --git a/ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.cpp b/ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.cpp
index 304ffa96647e9882b00bb010f9816d8e818d750e..db77bc7d96f80b30e532395e67294804c5c4c8d1 100644
--- a/ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.cpp
+++ b/ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.cpp
@@ -100,9 +100,8 @@ std::unique_ptr<Process> createThermoMechanicsProcess(
             variable_T->getNumberOfComponents());
     }
 
-    // Constitutive relation.
-    auto material =
-        MaterialLib::Solids::createConstitutiveRelation<DisplacementDim>(
+    auto solid_constitutive_relations =
+        MaterialLib::Solids::createConstitutiveRelations<DisplacementDim>(
             parameters, config);
 
     // Reference solid density
@@ -153,11 +152,9 @@ std::unique_ptr<Process> createThermoMechanicsProcess(
     }
 
     ThermoMechanicsProcessData<DisplacementDim> process_data{
-        std::move(material),
-        reference_solid_density,
-        linear_thermal_expansion_coefficient,
-        specific_heat_capacity,
-        thermal_conductivity,
+        materialIDs(mesh),       std::move(solid_constitutive_relations),
+        reference_solid_density, linear_thermal_expansion_coefficient,
+        specific_heat_capacity,  thermal_conductivity,
         specific_body_force};
 
     SecondaryVariableCollection secondary_variables;
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h b/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h
index 0cae1175f97cb6abdd68cc1bc9128f9941e0b963..76ef5ce9e782375dc9804918e68b979c52ba30a5 100644
--- a/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h
@@ -12,7 +12,7 @@
 #include <memory>
 #include <vector>
 
-#include "MaterialLib/SolidModels/MechanicsBase.h"
+#include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
 #include "MathLib/KelvinVector.h"
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "NumLib/Extrapolation/ExtrapolatableElement.h"
@@ -35,7 +35,8 @@ template <typename BMatricesType, typename ShapeMatricesType,
 struct IntegrationPointData final
 {
     explicit IntegrationPointData(
-        MaterialLib::Solids::MechanicsBase<DisplacementDim>& solid_material)
+        MaterialLib::Solids::MechanicsBase<DisplacementDim> const&
+            solid_material)
         : solid_material(solid_material),
           material_state_variables(
               solid_material.createMaterialStateVariables())
@@ -46,7 +47,7 @@ struct IntegrationPointData final
     typename BMatricesType::KelvinVectorType eps;
     typename BMatricesType::KelvinVectorType eps_m, eps_m_prev;
 
-    MaterialLib::Solids::MechanicsBase<DisplacementDim>& solid_material;
+    MaterialLib::Solids::MechanicsBase<DisplacementDim> const& solid_material;
     std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
         DisplacementDim>::MaterialStateVariables>
         material_state_variables;
@@ -123,10 +124,15 @@ public:
                               IntegrationMethod, DisplacementDim>(
                 e, is_axially_symmetric, _integration_method);
 
+        auto& solid_material =
+            MaterialLib::Solids::selectSolidConstitutiveRelation(
+                _process_data.solid_materials,
+                _process_data.material_ids,
+                e.getID());
+
         for (unsigned ip = 0; ip < n_integration_points; ip++)
         {
-            // displacement (subscript u)
-            _ip_data.emplace_back(*_process_data.material);
+            _ip_data.emplace_back(solid_material);
             auto& ip_data = _ip_data[ip];
             ip_data.integration_weight =
                 _integration_method.getWeightedPoint(ip).getWeight() *
@@ -327,14 +333,15 @@ public:
             //
             KuT.noalias() +=
                 B.transpose() * C * alpha * Invariants::identity2 * N * w;
-            if (_process_data.material->getConstitutiveModel() ==
+            if (_ip_data[ip].solid_material.getConstitutiveModel() ==
                 MaterialLib::Solids::ConstitutiveModel::CreepBGRa)
             {
                 auto const s = Invariants::deviatoric_projection * sigma;
                 double const norm_s = Invariants::FrobeniusNorm(s);
                 const double creep_coefficient =
-                    _process_data.material->getTemperatureRelatedCoefficient(
-                        t, dt, x_position, T_ip, norm_s);
+                    _ip_data[ip]
+                        .solid_material.getTemperatureRelatedCoefficient(
+                            t, dt, x_position, T_ip, norm_s);
                 KuT.noalias() += creep_coefficient * B.transpose() * s * N * w;
             }
 
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsProcessData.h b/ProcessLib/ThermoMechanics/ThermoMechanicsProcessData.h
index f369be3c35cc8a30bb90f0a3d129b662c1e562e0..63b205285cd4a4a5714ba4990f58e1420bafa773 100644
--- a/ProcessLib/ThermoMechanics/ThermoMechanicsProcessData.h
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsProcessData.h
@@ -32,14 +32,18 @@ template <int DisplacementDim>
 struct ThermoMechanicsProcessData
 {
     ThermoMechanicsProcessData(
-        std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>&&
-            material_,
+        MeshLib::PropertyVector<int> const* const material_ids_,
+        std::map<int,
+                 std::unique_ptr<
+                     MaterialLib::Solids::MechanicsBase<DisplacementDim>>>&&
+            solid_materials_,
         Parameter<double> const& reference_solid_density_,
         Parameter<double> const& linear_thermal_expansion_coefficient_,
         Parameter<double> const& specific_heat_capacity_,
         Parameter<double> const& thermal_conductivity_,
         Eigen::Matrix<double, DisplacementDim, 1> const& specific_body_force_)
-        : material{std::move(material_)},
+        : material_ids(material_ids_),
+          solid_materials{std::move(solid_materials_)},
           reference_solid_density(reference_solid_density_),
           linear_thermal_expansion_coefficient(
               linear_thermal_expansion_coefficient_),
@@ -49,18 +53,7 @@ struct ThermoMechanicsProcessData
     {
     }
 
-    ThermoMechanicsProcessData(ThermoMechanicsProcessData&& other)
-        : material{std::move(other.material)},
-          reference_solid_density(other.reference_solid_density),
-          linear_thermal_expansion_coefficient(
-              other.linear_thermal_expansion_coefficient),
-          specific_heat_capacity(other.specific_heat_capacity),
-          thermal_conductivity(other.thermal_conductivity),
-          specific_body_force(other.specific_body_force),
-          dt(other.dt),
-          t(other.t)
-    {
-    }
+    ThermoMechanicsProcessData(ThermoMechanicsProcessData&& other) = default;
 
     //! Copies are forbidden.
     ThermoMechanicsProcessData(ThermoMechanicsProcessData const&) = delete;
@@ -71,8 +64,14 @@ struct ThermoMechanicsProcessData
     //! Assignments are not needed.
     void operator=(ThermoMechanicsProcessData&&) = delete;
 
-    std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>
-        material;
+    MeshLib::PropertyVector<int> const* const material_ids;
+
+    /// The constitutive relation for the mechanical part.
+    std::map<
+        int,
+        std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>>
+        solid_materials;
+
     Parameter<double> const& reference_solid_density;
     Parameter<double> const& linear_thermal_expansion_coefficient;
     Parameter<double> const& specific_heat_capacity;
diff --git a/Tests/Data/Mechanics/Linear/expected_two_material_gravity_Emodulus.vtu b/Tests/Data/Mechanics/Linear/expected_two_material_gravity_Emodulus.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..e2b393d0bb576a47fcbf7697738ba3c1ba49bfc1
--- /dev/null
+++ b/Tests/Data/Mechanics/Linear/expected_two_material_gravity_Emodulus.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:2c53b1a459bca2030d2423acf52828beb4f245e31adca471a1819d01aa246c2d
+size 611781
diff --git a/Tests/Data/Mechanics/Linear/two_material_gravity_Emodulus.prj b/Tests/Data/Mechanics/Linear/two_material_gravity_Emodulus.prj
new file mode 100644
index 0000000000000000000000000000000000000000..0a41d8b26cc7f4a4ba2ffd30d6464e1ec5b2bc40
--- /dev/null
+++ b/Tests/Data/Mechanics/Linear/two_material_gravity_Emodulus.prj
@@ -0,0 +1,161 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+    <mesh>two_material_gravity_Emodulus.vtu</mesh>
+    <geometry>cube_1x1x1.gml</geometry>
+    <processes>
+        <process>
+            <name>SD</name>
+            <type>SMALL_DEFORMATION</type>
+            <integration_order>2</integration_order>
+            <constitutive_relation id="0">
+                <type>LinearElasticIsotropic</type>
+                <youngs_modulus>E0</youngs_modulus>
+                <poissons_ratio>nu</poissons_ratio>
+            </constitutive_relation>
+            <constitutive_relation id="1">
+                <type>LinearElasticIsotropic</type>
+                <youngs_modulus>E1</youngs_modulus>
+                <poissons_ratio>nu</poissons_ratio>
+            </constitutive_relation>
+            <solid_density>rho_sr</solid_density>
+            <specific_body_force>0 0 -1</specific_body_force>
+            <process_variables>
+                <process_variable>displacement</process_variable>
+            </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>
+        </process>
+    </processes>
+    <time_loop>
+        <processes>
+            <process ref="SD">
+                <nonlinear_solver>basic_newton</nonlinear_solver>
+                <convergence_criterion>
+                    <type>DeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <abstol>1e-13</abstol>
+                </convergence_criterion>
+                <time_discretization>
+                    <type>BackwardEuler</type>
+                </time_discretization>
+                <output>
+                    <variables>
+                        <variable>displacement</variable>
+                        <variable>sigma</variable>
+                        <variable>epsilon</variable>
+                    </variables>
+                </output>
+                <time_stepping>
+                    <type>FixedTimeStepping</type>
+                    <t_initial>0</t_initial>
+                    <t_end>1</t_end>
+                    <timesteps>
+                        <pair>
+                            <repeat>1</repeat>
+                            <delta_t>1</delta_t>
+                        </pair>
+                    </timesteps>
+                </time_stepping>
+            </process>
+        </processes>
+        <output>
+            <type>VTK</type>
+            <prefix>two_material_gravity_Emodulus</prefix>
+            <timesteps>
+                <pair>
+                    <repeat>1</repeat>
+                    <each_steps>10000000</each_steps>
+                </pair>
+            </timesteps>
+        </output>
+    </time_loop>
+    <parameters>
+        <parameter>
+            <name>E0</name>
+            <type>Constant</type>
+            <value>1000</value>
+        </parameter>
+        <parameter>
+            <name>E1</name>
+            <type>Constant</type>
+            <value>10000</value>
+        </parameter>
+        <parameter>
+            <name>nu</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+        <parameter>
+            <name>rho_sr</name>
+            <type>Constant</type>
+            <value>1</value>
+        </parameter>
+        <parameter>
+            <name>displacement0</name>
+            <type>Constant</type>
+            <values>0 0 0</values>
+        </parameter>
+        <parameter>
+            <name>dirichlet0</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+    </parameters>
+    <process_variables>
+        <process_variable>
+            <name>displacement</name>
+            <components>3</components>
+            <order>1</order>
+            <initial_condition>displacement0</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+                    <geometry>bottom</geometry>
+                    <type>Dirichlet</type>
+                    <component>0</component>
+                    <parameter>dirichlet0</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+                    <geometry>bottom</geometry>
+                    <type>Dirichlet</type>
+                    <component>1</component>
+                    <parameter>dirichlet0</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+                    <geometry>bottom</geometry>
+                    <type>Dirichlet</type>
+                    <component>2</component>
+                    <parameter>dirichlet0</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+    </process_variables>
+    <nonlinear_solvers>
+        <nonlinear_solver>
+            <name>basic_newton</name>
+            <type>Newton</type>
+            <max_iter>4</max_iter>
+            <linear_solver>general_linear_solver</linear_solver>
+        </nonlinear_solver>
+    </nonlinear_solvers>
+    <linear_solvers>
+        <linear_solver>
+            <name>general_linear_solver</name>
+            <lis>-i cg -p jacobi -tol 1e-16 -maxiter 10000</lis>
+            <eigen>
+                <solver_type>CG</solver_type>
+                <precon_type>DIAGONAL</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/Mechanics/Linear/two_material_gravity_Emodulus.vtu b/Tests/Data/Mechanics/Linear/two_material_gravity_Emodulus.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..aebeee38c560de47543d802214683f22b06f31a8
--- /dev/null
+++ b/Tests/Data/Mechanics/Linear/two_material_gravity_Emodulus.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:74078c2e63b1516f25fad9ec16bfcde0af6d43a1ae020ac40209ffc84ecfeb1b
+size 19892
diff --git a/Tests/Data/RichardsMechanics/mechanics_linear.prj b/Tests/Data/RichardsMechanics/mechanics_linear.prj
index 6d0e8ae8c35a2d9c296bea6de3e29dafe4da30c5..fe83ced3a89dcbd670552d3bb3881ff82f12227f 100644
--- a/Tests/Data/RichardsMechanics/mechanics_linear.prj
+++ b/Tests/Data/RichardsMechanics/mechanics_linear.prj
@@ -21,7 +21,7 @@
                 <relative_epsilons>1e-8 1e-8 1e-8 1e-8 1e-8 1e-8 1e-8 1e-8 1e-8 1e-8 1e-8 1e-8 1e-8 1e-8 1e-8 1e-8 1e-8 1e-8 1e-8 1e-8</relative_epsilons>
             </jacobian_assembler>
             -->
-            <constitutive_relation>
+            <constitutive_relation id="0">
                 <type>LinearElasticIsotropic</type>
                 <youngs_modulus>E</youngs_modulus>
                 <poissons_ratio>nu</poissons_ratio>
diff --git a/Tests/Data/ThermoMechanics/CreepBGRa/CreepAfterExcavation/CreepAfterExcavation.prj b/Tests/Data/ThermoMechanics/CreepBGRa/CreepAfterExcavation/CreepAfterExcavation.prj
index 1e00e2049f10e33f1844e0f77b470c43f53c486f..89f9157ff63c2baebe323d28e243cd7c8b677c96 100644
--- a/Tests/Data/ThermoMechanics/CreepBGRa/CreepAfterExcavation/CreepAfterExcavation.prj
+++ b/Tests/Data/ThermoMechanics/CreepBGRa/CreepAfterExcavation/CreepAfterExcavation.prj
@@ -11,10 +11,15 @@
             <name>ThermoMechanics</name>
             <type>THERMO_MECHANICS</type>
             <integration_order>2</integration_order>
-            <constitutive_relation>
+            <constitutive_relation id="0">
+                <type>LinearElasticIsotropic</type>
+                <youngs_modulus>E0</youngs_modulus>
+                <poissons_ratio>nu0</poissons_ratio>
+            </constitutive_relation>
+            <constitutive_relation id="1">
                 <type>CreepBGRa</type>
-                <youngs_modulus>E</youngs_modulus>
-                <poissons_ratio>nu</poissons_ratio>
+                <youngs_modulus>E1</youngs_modulus>
+                <poissons_ratio>nu1</poissons_ratio>
                 <a>A</a>
                 <n>n</n>
                 <sigma0>sigma_f</sigma0>
@@ -99,8 +104,8 @@
     <parameters>
         <parameter>
             <name>A</name>
-            <type>MeshElement</type>
-            <field_name>A</field_name>
+            <type>Constant</type>
+            <value>2.0833333333333333e-6</value>
         </parameter>
         <parameter>
             <name>n</name>
@@ -119,14 +124,24 @@
         </parameter>
 
         <parameter>
-            <name>E</name>
-            <type>MeshElement</type>
-            <field_name>E</field_name>
+            <name>E1</name>
+            <type>Constant</type>
+            <value>7.65e9</value>
         </parameter>
         <parameter>
-            <name>nu</name>
-            <type>MeshElement</type>
-            <field_name>nu</field_name>
+            <name>nu1</name>
+            <type>Constant</type>
+            <value>0.27</value>
+        </parameter>
+        <parameter>
+            <name>nu0</name>
+            <type>Constant</type>
+            <value>0.3</value>
+        </parameter>
+        <parameter>
+            <name>E0</name>
+            <type>Constant</type>
+            <value>7e9</value>
         </parameter>
         <parameter>
             <name>rho_sr</name>
diff --git a/Tests/Data/ThermoMechanics/CreepBGRa/CreepAfterExcavation/mesh.vtu b/Tests/Data/ThermoMechanics/CreepBGRa/CreepAfterExcavation/mesh.vtu
index 184d69d8a9345843f5c0f85343670878ad960150..93871b5c8c4cb695527933d0e18d0381d6c08d39 100644
--- a/Tests/Data/ThermoMechanics/CreepBGRa/CreepAfterExcavation/mesh.vtu
+++ b/Tests/Data/ThermoMechanics/CreepBGRa/CreepAfterExcavation/mesh.vtu
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:63351c97e8fea8910c6c59faec743b084b1950557086c58eaee19e66b533d8ad
-size 255945
+oid sha256:267930c68c85c52515682459994b19ea3fe69d233964be8639d80d39f1d311d6
+size 60171