diff --git a/MaterialLib/MPL/Component.cpp b/MaterialLib/MPL/Component.cpp
index 83fa318835160dbc5e04cd86e9c993919adc8ec8..96c44fcfe762d715aa041236656874b43a0ff140 100644
--- a/MaterialLib/MPL/Component.cpp
+++ b/MaterialLib/MPL/Component.cpp
@@ -31,9 +31,21 @@ Component::Component(std::string const& component_name,
 
 Property const& Component::property(PropertyType const& p) const
 {
+    Property const* const property = properties_[p].get();
+    if (property == nullptr)
+    {
+        OGS_FATAL("Trying to access undefined property '{:s}' of {:s}",
+                  property_enum_to_string[p], description());
+    }
+
     return *properties_[p];
 }
 
+Property const& Component::operator[](PropertyType const& p) const
+{
+    return property(p);
+}
+
 bool Component::hasProperty(PropertyType const& p) const
 {
     return properties_[p] != nullptr;
diff --git a/MaterialLib/MPL/Component.h b/MaterialLib/MPL/Component.h
index df3eadf0f02a48080fd9260faf167dbec2ea1931..7a995388ac49d2bf9131b8058809c278db0c74d4 100644
--- a/MaterialLib/MPL/Component.h
+++ b/MaterialLib/MPL/Component.h
@@ -34,6 +34,9 @@ public:
 
     /// A get-function for retrieving a certain property.
     Property const& property(PropertyType const& /*p*/) const;
+
+    Property const& operator[](PropertyType const& p) const;
+
     bool hasProperty(PropertyType const& p) const;
 
     template <typename T>
diff --git a/MaterialLib/MPL/Medium.cpp b/MaterialLib/MPL/Medium.cpp
index 2e99adaab1d427b26da8dabd9b48e3adbdac2e6a..0deb92d67afe8bce67e2afb8771d66dea5b92399 100644
--- a/MaterialLib/MPL/Medium.cpp
+++ b/MaterialLib/MPL/Medium.cpp
@@ -61,6 +61,11 @@ Property const& Medium::property(PropertyType const& p) const
     return *properties_[p];
 }
 
+Property const& Medium::operator[](PropertyType const& p) const
+{
+    return property(p);
+}
+
 bool Medium::hasProperty(PropertyType const& p) const
 {
     return properties_[p] != nullptr;
diff --git a/MaterialLib/MPL/Medium.h b/MaterialLib/MPL/Medium.h
index bfc0d290b8506b4d1e2e87a3d07ff65170d94bf6..aa79781c8d04a07bc156fed451bd3c6e1df13e01 100644
--- a/MaterialLib/MPL/Medium.h
+++ b/MaterialLib/MPL/Medium.h
@@ -45,6 +45,9 @@ public:
     /// A get-function for a property. The argument refers to the name of the
     /// property.
     Property const& property(PropertyType const& p) const;
+
+    Property const& operator[](PropertyType const& p) const;
+
     bool hasProperty(PropertyType const& p) const;
 
     /// A simple get-function for retrieving the number of phases the medium
diff --git a/MaterialLib/MPL/Phase.cpp b/MaterialLib/MPL/Phase.cpp
index 68921e4770c9dfd2162101eadb407c8379b3d822..f65c660d80e5ad80bd19875b4a39a867a177e23f 100644
--- a/MaterialLib/MPL/Phase.cpp
+++ b/MaterialLib/MPL/Phase.cpp
@@ -60,6 +60,11 @@ Property const& Phase::property(PropertyType const& p) const
     return *properties_[p];
 }
 
+Property const& Phase::operator[](PropertyType const& p) const
+{
+    return property(p);
+}
+
 bool Phase::hasProperty(PropertyType const& p) const
 {
     return properties_[p] != nullptr;
diff --git a/MaterialLib/MPL/Phase.h b/MaterialLib/MPL/Phase.h
index b4fbe03590a79ff03802a5d2c6555738dccd8167..141c3eb4a9d363cbdc6fb1f392ba307fda4cf3a8 100644
--- a/MaterialLib/MPL/Phase.h
+++ b/MaterialLib/MPL/Phase.h
@@ -46,6 +46,9 @@ public:
     /// A get-function for a property. The argument refers to the name of the
     /// property.
     Property const& property(PropertyType const& p) const;
+
+    Property const& operator[](PropertyType const& p) const;
+
     bool hasProperty(PropertyType const& p) const;
 
     /// A get-function for retrieving the number of components in this phase.
diff --git a/ProcessLib/ComponentTransport/ComponentTransportFEM.h b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
index f98ee53784fd65035274b14912782947c4c1ffdf..76fc6bdd6a5f9af990e356fdff3bdbdab540c131 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportFEM.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
@@ -206,7 +206,7 @@ public:
                                       GlobalDim>(element, is_axially_symmetric,
                                                  _integration_method);
         auto const& medium =
-            _process_data.media_map->getMedium(_element.getID());
+            *_process_data.media_map->getMedium(_element.getID());
         for (unsigned ip = 0; ip < n_integration_points; ip++)
         {
             _ip_data.emplace_back(
@@ -218,7 +218,7 @@ public:
             pos.setIntegrationPoint(ip);
 
             _ip_data[ip].porosity =
-                medium->property(MaterialPropertyLib::PropertyType::porosity)
+                medium[MaterialPropertyLib::PropertyType::porosity]
                     .template initialValue<double>(
                         pos, std::numeric_limits<double>::quiet_NaN() /*t*/);
 
@@ -474,16 +474,13 @@ public:
                 MaterialPropertyLib::Variable::phase_pressure)] = p_int_pt;
 
             // update according to a particular porosity model
-            porosity =
-                medium.property(MaterialPropertyLib::PropertyType::porosity)
-                    .template value<double>(vars, pos, t, dt);
+            porosity = medium[MaterialPropertyLib::PropertyType::porosity]
+                           .template value<double>(vars, pos, t, dt);
             vars[static_cast<int>(MaterialPropertyLib::Variable::porosity)] =
                 porosity;
 
             auto const& retardation_factor =
-                component
-                    .property(
-                        MaterialPropertyLib::PropertyType::retardation_factor)
+                component[MaterialPropertyLib::PropertyType::retardation_factor]
                     .template value<double>(vars, pos, t, dt);
 
             auto const& solute_dispersivity_transverse = medium.template value<
@@ -499,29 +496,25 @@ public:
             // TODO (renchao): concentration of which component as the argument
             // for calculation of fluid density
             auto const density =
-                phase.property(MaterialPropertyLib::PropertyType::density)
+                phase[MaterialPropertyLib::PropertyType::density]
                     .template value<double>(vars, pos, t, dt);
 
             auto const decay_rate =
-                component
-                    .property(MaterialPropertyLib::PropertyType::decay_rate)
+                component[MaterialPropertyLib::PropertyType::decay_rate]
                     .template value<double>(vars, pos, t, dt);
 
             auto const& pore_diffusion_coefficient =
                 MaterialPropertyLib::formEigenTensor<GlobalDim>(
-                    component
-                        .property(
-                            MaterialPropertyLib::PropertyType::pore_diffusion)
+                    component[MaterialPropertyLib::PropertyType::pore_diffusion]
                         .value(vars, pos, t, dt));
 
             auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
-                medium.property(MaterialPropertyLib::PropertyType::permeability)
-                    .value(vars, pos, t, dt));
+                medium[MaterialPropertyLib::PropertyType::permeability].value(
+                    vars, pos, t, dt));
 
             // Use the viscosity model to compute the viscosity
-            auto const mu =
-                phase.property(MaterialPropertyLib::PropertyType::viscosity)
-                    .template value<double>(vars, pos, t, dt);
+            auto const mu = phase[MaterialPropertyLib::PropertyType::viscosity]
+                                .template value<double>(vars, pos, t, dt);
 
             GlobalDimMatrixType const K_over_mu = K / mu;
             GlobalDimVectorType const velocity =
@@ -531,13 +524,13 @@ public:
                     : GlobalDimVectorType(-K_over_mu * dNdx * p_nodal_values);
 
             const double drho_dp =
-                phase.property(MaterialPropertyLib::PropertyType::density)
+                phase[MaterialPropertyLib::PropertyType::density]
                     .template dValue<double>(
                         vars, MaterialPropertyLib::Variable::phase_pressure,
                         pos, t, dt);
 
             const double drho_dC =
-                phase.property(MaterialPropertyLib::PropertyType::density)
+                phase[MaterialPropertyLib::PropertyType::density]
                     .template dValue<double>(
                         vars, MaterialPropertyLib::Variable::concentration, pos,
                         t, dt);
@@ -683,9 +676,7 @@ public:
                 porosity =
                     _process_data.chemically_induced_porosity_change
                         ? porosity_prev
-                        : medium
-                              .property(
-                                  MaterialPropertyLib::PropertyType::porosity)
+                        : medium[MaterialPropertyLib::PropertyType::porosity]
                               .template value<double>(vars, vars_prev, pos, t,
                                                       dt);
 
@@ -697,27 +688,26 @@ public:
             // TODO: Concentration of which component as one of arguments for
             // calculation of fluid density
             auto const density =
-                phase.property(MaterialPropertyLib::PropertyType::density)
+                phase[MaterialPropertyLib::PropertyType::density]
                     .template value<double>(vars, pos, t, dt);
 
             auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
-                medium.property(MaterialPropertyLib::PropertyType::permeability)
-                    .value(vars, pos, t, dt));
+                medium[MaterialPropertyLib::PropertyType::permeability].value(
+                    vars, pos, t, dt));
 
             // Use the viscosity model to compute the viscosity
-            auto const mu =
-                phase.property(MaterialPropertyLib::PropertyType::viscosity)
-                    .template value<double>(vars, pos, t, dt);
+            auto const mu = phase[MaterialPropertyLib::PropertyType::viscosity]
+                                .template value<double>(vars, pos, t, dt);
 
             GlobalDimMatrixType const K_over_mu = K / mu;
 
             const double drho_dp =
-                phase.property(MaterialPropertyLib::PropertyType::density)
+                phase[MaterialPropertyLib::PropertyType::density]
                     .template dValue<double>(
                         vars, MaterialPropertyLib::Variable::phase_pressure,
                         pos, t, dt);
             const double drho_dC =
-                phase.property(MaterialPropertyLib::PropertyType::density)
+                phase[MaterialPropertyLib::PropertyType::density]
                     .template dValue<double>(
                         vars, MaterialPropertyLib::Variable::concentration, pos,
                         t, dt);
@@ -816,9 +806,7 @@ public:
                 porosity =
                     _process_data.chemically_induced_porosity_change
                         ? porosity_prev
-                        : medium
-                              .property(
-                                  MaterialPropertyLib::PropertyType::porosity)
+                        : medium[MaterialPropertyLib::PropertyType::porosity]
                               .template value<double>(vars, vars_prev, pos, t,
                                                       dt);
 
@@ -827,9 +815,7 @@ public:
             }
 
             auto const& retardation_factor =
-                component
-                    .property(
-                        MaterialPropertyLib::PropertyType::retardation_factor)
+                component[MaterialPropertyLib::PropertyType::retardation_factor]
                     .template value<double>(vars, pos, t, dt);
 
             auto const& solute_dispersivity_transverse = medium.template value<
@@ -842,27 +828,23 @@ public:
 
             // Use the fluid density model to compute the density
             auto const density =
-                phase.property(MaterialPropertyLib::PropertyType::density)
+                phase[MaterialPropertyLib::PropertyType::density]
                     .template value<double>(vars, pos, t, dt);
             auto const decay_rate =
-                component
-                    .property(MaterialPropertyLib::PropertyType::decay_rate)
+                component[MaterialPropertyLib::PropertyType::decay_rate]
                     .template value<double>(vars, pos, t, dt);
 
             auto const& pore_diffusion_coefficient =
                 MaterialPropertyLib::formEigenTensor<GlobalDim>(
-                    component
-                        .property(
-                            MaterialPropertyLib::PropertyType::pore_diffusion)
+                    component[MaterialPropertyLib::PropertyType::pore_diffusion]
                         .value(vars, pos, t, dt));
 
             auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
-                medium.property(MaterialPropertyLib::PropertyType::permeability)
-                    .value(vars, pos, t, dt));
+                medium[MaterialPropertyLib::PropertyType::permeability].value(
+                    vars, pos, t, dt));
             // Use the viscosity model to compute the viscosity
-            auto const mu =
-                phase.property(MaterialPropertyLib::PropertyType::viscosity)
-                    .template value<double>(vars, pos, t, dt);
+            auto const mu = phase[MaterialPropertyLib::PropertyType::viscosity]
+                                .template value<double>(vars, pos, t, dt);
 
             GlobalDimMatrixType const K_over_mu = K / mu;
             GlobalDimVectorType const velocity =
@@ -893,7 +875,7 @@ public:
             if (_process_data.non_advective_form)
             {
                 const double drho_dC =
-                    phase.property(MaterialPropertyLib::PropertyType::density)
+                    phase[MaterialPropertyLib::PropertyType::density]
                         .template dValue<double>(
                             vars, MaterialPropertyLib::Variable::concentration,
                             pos, t, dt);
@@ -911,7 +893,7 @@ public:
 
                 NumLib::shapeFunctionInterpolate(local_pdot, N, dot_p_int_pt);
                 const double drho_dp =
-                    phase.property(MaterialPropertyLib::PropertyType::density)
+                    phase[MaterialPropertyLib::PropertyType::density]
                         .template dValue<double>(
                             vars, MaterialPropertyLib::Variable::phase_pressure,
                             pos, t, dt);
@@ -1023,18 +1005,17 @@ public:
             // models. Need extension of secondary variables interface.
             double const dt = std::numeric_limits<double>::quiet_NaN();
             auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
-                medium.property(MaterialPropertyLib::PropertyType::permeability)
-                    .value(vars, pos, t, dt));
-            auto const mu =
-                phase.property(MaterialPropertyLib::PropertyType::viscosity)
-                    .template value<double>(vars, pos, t, dt);
+                medium[MaterialPropertyLib::PropertyType::permeability].value(
+                    vars, pos, t, dt));
+            auto const mu = phase[MaterialPropertyLib::PropertyType::viscosity]
+                                .template value<double>(vars, pos, t, dt);
             GlobalDimMatrixType const K_over_mu = K / mu;
 
             cache_mat.col(ip).noalias() = -K_over_mu * dNdx * p_nodal_values;
             if (_process_data.has_gravity)
             {
                 auto const rho_w =
-                    phase.property(MaterialPropertyLib::PropertyType::density)
+                    phase[MaterialPropertyLib::PropertyType::density]
                         .template value<double>(vars, pos, t, dt);
                 auto const b = _process_data.specific_body_force;
                 // here it is assumed that the vector b is directed 'downwards'
@@ -1096,18 +1077,16 @@ public:
         // Need extension of secondary variables interface.
         double const dt = std::numeric_limits<double>::quiet_NaN();
         auto const K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
-            medium.property(MaterialPropertyLib::PropertyType::permeability)
-                .value(vars, pos, t, dt));
+            medium[MaterialPropertyLib::PropertyType::permeability].value(
+                vars, pos, t, dt));
 
-        auto const mu =
-            phase.property(MaterialPropertyLib::PropertyType::viscosity)
-                .template value<double>(vars, pos, t, dt);
+        auto const mu = phase[MaterialPropertyLib::PropertyType::viscosity]
+                            .template value<double>(vars, pos, t, dt);
         GlobalDimMatrixType const K_over_mu = K / mu;
 
         GlobalDimVectorType q = -K_over_mu * shape_matrices.dNdx * local_p;
-        auto const rho_w =
-            phase.property(MaterialPropertyLib::PropertyType::density)
-                .template value<double>(vars, pos, t, dt);
+        auto const rho_w = phase[MaterialPropertyLib::PropertyType::density]
+                               .template value<double>(vars, pos, t, dt);
         if (_process_data.has_gravity)
         {
             auto const b = _process_data.specific_body_force;