diff --git a/Applications/ApplicationsLib/ProjectData.h b/Applications/ApplicationsLib/ProjectData.h
index 53704d198b746fd0eaf0326e96d804d744193dd6..d309380822c5aa15100b1f9db95655ba04546aa2 100644
--- a/Applications/ApplicationsLib/ProjectData.h
+++ b/Applications/ApplicationsLib/ProjectData.h
@@ -91,6 +91,12 @@ public:
 
     MeshLib::Mesh* getMesh(std::string const& mesh_name) const;
 
+    std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const&
+    getMedia() const
+    {
+        return _media;
+    }
+
 private:
     /// Parses the process variables configuration and creates new variables for
     /// each variable entry passing the corresponding subtree to the process
diff --git a/Applications/ApplicationsLib/Simulation.cpp b/Applications/ApplicationsLib/Simulation.cpp
index 7b5656150b8d3077af7d959c09ba81d3a4a88404..53f07cba841f6fa033d7719f04dad9b221d7a44c 100644
--- a/Applications/ApplicationsLib/Simulation.cpp
+++ b/Applications/ApplicationsLib/Simulation.cpp
@@ -101,7 +101,7 @@ void Simulation::initializeDataStructures(
     INFO("Initialize processes.");
     for (auto& p : project_data->getProcesses())
     {
-        p->initialize();
+        p->initialize(project_data->getMedia());
     }
 
     // Check intermediately that config parsing went fine.
diff --git a/ProcessLib/BoundaryConditionAndSourceTerm/BoundaryConditionCollection.cpp b/ProcessLib/BoundaryConditionAndSourceTerm/BoundaryConditionCollection.cpp
index 769f29d8f1ba59b548c4bd7208f294fa15bad77c..9288714d7e8e14e6979125de661c0fe42b67f2e3 100644
--- a/ProcessLib/BoundaryConditionAndSourceTerm/BoundaryConditionCollection.cpp
+++ b/ProcessLib/BoundaryConditionAndSourceTerm/BoundaryConditionCollection.cpp
@@ -26,16 +26,17 @@ void BoundaryConditionCollection::addBCsForProcessVariables(
     std::vector<std::reference_wrapper<ProcessVariable>> const&
         process_variables,
     NumLib::LocalToGlobalIndexMap const& dof_table,
-    unsigned const integration_order, Process const& process)
+    unsigned const integration_order, Process const& process,
+    std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media)
 {
     for (int variable_id = 0;
          variable_id < static_cast<int>(process_variables.size());
          ++variable_id)
     {
         ProcessVariable& pv = process_variables[variable_id];
-        auto bcs = pv.createBoundaryConditions(dof_table, variable_id,
-                                               integration_order, _parameters,
-                                               process, process_variables);
+        auto bcs = pv.createBoundaryConditions(
+            dof_table, variable_id, integration_order, _parameters, process,
+            process_variables, media);
 
         std::move(bcs.begin(), bcs.end(),
                   std::back_inserter(_boundary_conditions));
diff --git a/ProcessLib/BoundaryConditionAndSourceTerm/BoundaryConditionCollection.h b/ProcessLib/BoundaryConditionAndSourceTerm/BoundaryConditionCollection.h
index cdc9bbff0021ab743dcd8d5b63a48e2d15d0d321..16d4cb067e66c8f7d21ad790ddca8f2e6223bb91 100644
--- a/ProcessLib/BoundaryConditionAndSourceTerm/BoundaryConditionCollection.h
+++ b/ProcessLib/BoundaryConditionAndSourceTerm/BoundaryConditionCollection.h
@@ -11,6 +11,7 @@
 #pragma once
 
 #include "BoundaryCondition.h"
+#include "MaterialLib/MPL/Medium.h"
 #include "NumLib/IndexValueVector.h"
 #include "ProcessLib/ProcessVariable.h"
 
@@ -47,7 +48,9 @@ public:
         std::vector<std::reference_wrapper<ProcessVariable>> const&
             process_variables,
         NumLib::LocalToGlobalIndexMap const& dof_table,
-        unsigned const integration_order, Process const& process);
+        unsigned const integration_order, Process const& process,
+        std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const&
+            media);
 
     void addBoundaryCondition(std::unique_ptr<BoundaryCondition>&& bc)
     {
diff --git a/ProcessLib/BoundaryConditionAndSourceTerm/CreateBoundaryCondition.cpp b/ProcessLib/BoundaryConditionAndSourceTerm/CreateBoundaryCondition.cpp
index d9a8f5436b6352ebf7b5434813506521a4ff4133..d817590d8c62ac6e890cfe335635ea5f9b1a2b96 100644
--- a/ProcessLib/BoundaryConditionAndSourceTerm/CreateBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryConditionAndSourceTerm/CreateBoundaryCondition.cpp
@@ -36,7 +36,8 @@ std::unique_ptr<BoundaryCondition> createBoundaryCondition(
     const std::vector<std::unique_ptr<ParameterLib::ParameterBase>>& parameters,
     const Process& process,
     [[maybe_unused]] std::vector<std::reference_wrapper<ProcessVariable>> const&
-        all_process_variables_for_this_process)
+        all_process_variables_for_this_process,
+    std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& /*media*/)
 {
     // Surface mesh and bulk mesh must have equal axial symmetry flags!
     if (config.boundary_mesh.isAxiallySymmetric() !=
diff --git a/ProcessLib/BoundaryConditionAndSourceTerm/CreateBoundaryCondition.h b/ProcessLib/BoundaryConditionAndSourceTerm/CreateBoundaryCondition.h
index d3ba048f884b4b185ebfff5440b92db5f206cf3c..cd7fada96a3ba44c066296a4cf7cf93d923e8738 100644
--- a/ProcessLib/BoundaryConditionAndSourceTerm/CreateBoundaryCondition.h
+++ b/ProcessLib/BoundaryConditionAndSourceTerm/CreateBoundaryCondition.h
@@ -13,6 +13,7 @@
 #include <memory>
 #include <vector>
 
+#include "MaterialLib/MPL/CreateMaterialSpatialDistributionMap.h"
 namespace MeshLib
 {
 class Mesh;
@@ -41,6 +42,7 @@ std::unique_ptr<BoundaryCondition> createBoundaryCondition(
     const std::vector<std::unique_ptr<ParameterLib::ParameterBase>>& parameters,
     const Process& process,
     std::vector<std::reference_wrapper<ProcessVariable>> const&
-        all_process_variables_for_this_process);
+        all_process_variables_for_this_process,
+    std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media);
 
 }  // namespace ProcessLib
diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index e4b44bb6a9474ffe24a079bbd4bf995b7661c05e..7b56e5ee628ecab953a255fb5bf408a596b73eff 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -68,20 +68,22 @@ Process::Process(
 }
 
 void Process::initializeProcessBoundaryConditionsAndSourceTerms(
-    const NumLib::LocalToGlobalIndexMap& dof_table, const int process_id)
+    const NumLib::LocalToGlobalIndexMap& dof_table, const int process_id,
+    std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media)
 {
     auto const& per_process_variables = _process_variables[process_id];
     auto& per_process_BCs = _boundary_conditions[process_id];
 
     per_process_BCs.addBCsForProcessVariables(per_process_variables, dof_table,
-                                              _integration_order, *this);
+                                              _integration_order, *this, media);
 
     auto& per_process_sts = _source_term_collections[process_id];
     per_process_sts.addSourceTermsForProcessVariables(
         per_process_variables, dof_table, _integration_order);
 }
 
-void Process::initializeBoundaryConditions()
+void Process::initializeBoundaryConditions(
+    std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media)
 {
     // The number of processes is identical to the size of _process_variables,
     // the vector contains variables for different processes. See the
@@ -90,11 +92,12 @@ void Process::initializeBoundaryConditions()
     for (std::size_t pcs_id = 0; pcs_id < number_of_processes; pcs_id++)
     {
         initializeProcessBoundaryConditionsAndSourceTerms(
-            *_local_to_global_index_map, pcs_id);
+            *_local_to_global_index_map, pcs_id, media);
     }
 }
 
-void Process::initialize()
+void Process::initialize(
+    std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media)
 {
     DBUG("Initialize process.");
 
@@ -111,7 +114,7 @@ void Process::initialize()
                               _integration_order);
 
     DBUG("Initialize boundary conditions.");
-    initializeBoundaryConditions();
+    initializeBoundaryConditions(media);
 }
 
 void Process::setInitialConditions(
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index f02a2cc6a70972d4e5316e30b4c3faafaa98837e..8b5ccbff63fdca2038c8603cf113ca20bd501f72 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -14,6 +14,7 @@
 #include <tuple>
 
 #include "AbstractJacobianAssembler.h"
+#include "MaterialLib/MPL/Medium.h"
 #include "MathLib/LinAlg/GlobalMatrixVectorTypes.h"
 #include "MeshLib/Utils/IntegrationPointWriter.h"
 #include "NumLib/ODESolver/NonlinearSolver.h"
@@ -87,7 +88,9 @@ public:
 
     NumLib::IterationResult postIteration(GlobalVector const& x) final;
 
-    void initialize();
+    void initialize(
+        std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const&
+            media);
 
     void setInitialConditions(
         std::vector<GlobalVector*>& process_solutions,
@@ -205,7 +208,9 @@ protected:
      * initializeBoundaryConditions().
      */
     void initializeProcessBoundaryConditionsAndSourceTerms(
-        const NumLib::LocalToGlobalIndexMap& dof_table, const int process_id);
+        const NumLib::LocalToGlobalIndexMap& dof_table, const int process_id,
+        std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const&
+            media);
 
 private:
     /// Process specific initialization called by initialize().
@@ -216,7 +221,9 @@ private:
 
     /// Member function to initialize the boundary conditions for all coupled
     /// processes. It is called by initialize().
-    virtual void initializeBoundaryConditions();
+    virtual void initializeBoundaryConditions(
+        std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const&
+            media);
 
     virtual void setInitialConditionsConcreteProcess(
         std::vector<GlobalVector*>& /*x*/,
diff --git a/ProcessLib/ProcessVariable.cpp b/ProcessLib/ProcessVariable.cpp
index 82946e3eb37fcbf9e936cd60758c6d77f7315255..d1fa21db0fa6f7ada9a5481db527c05685baeb0c 100644
--- a/ProcessLib/ProcessVariable.cpp
+++ b/ProcessLib/ProcessVariable.cpp
@@ -225,7 +225,8 @@ ProcessVariable::createBoundaryConditions(
     std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
     Process const& process,
     std::vector<std::reference_wrapper<ProcessVariable>> const&
-        all_process_variables_for_this_process)
+        all_process_variables_for_this_process,
+    std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media)
 {
     std::vector<std::unique_ptr<BoundaryCondition>> bcs;
     bcs.reserve(_bc_configs.size());
@@ -235,7 +236,7 @@ ProcessVariable::createBoundaryConditions(
         auto bc = createBoundaryCondition(
             config, dof_table, _mesh, variable_id, integration_order,
             _shapefunction_order, parameters, process,
-            all_process_variables_for_this_process);
+            all_process_variables_for_this_process, media);
 #ifdef USE_PETSC
         if (bc == nullptr)
         {
diff --git a/ProcessLib/ProcessVariable.h b/ProcessLib/ProcessVariable.h
index fcdc12d6738255630287b5417aa5d17c22138f2d..7d743ca76abcd7cf0411174763cd1095835be904 100644
--- a/ProcessLib/ProcessVariable.h
+++ b/ProcessLib/ProcessVariable.h
@@ -16,6 +16,7 @@
 #include <vector>
 
 #include "BaseLib/ConfigTree-fwd.h"
+#include "MaterialLib/MPL/Medium.h"
 
 namespace MathLib
 {
@@ -98,7 +99,9 @@ public:
             parameters,
         Process const& process,
         std::vector<std::reference_wrapper<ProcessVariable>> const&
-            all_process_variables_for_this_process);
+            all_process_variables_for_this_process,
+        std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const&
+            media);
 
     std::vector<std::unique_ptr<SourceTerm>> createSourceTerms(
         const NumLib::LocalToGlobalIndexMap& dof_table, const int variable_id,