From b7ce1bc30136e5981cc0504f1bd5ba59a5fb837c Mon Sep 17 00:00:00 2001
From: Thomas Fischer <thomas.fischer@ufz.de>
Date: Thu, 23 Jul 2020 14:03:57 +0200
Subject: [PATCH] [PL/HeatConduction] Use MPL for thermal_conductivity.

---
 Applications/ApplicationsLib/ProjectData.cpp  |  2 +-
 .../CreateHeatConductionProcess.cpp           | 32 ++++++++----
 .../CreateHeatConductionProcess.h             |  8 ++-
 ProcessLib/HeatConduction/HeatConductionFEM.h | 52 +++++++++++++++++--
 .../HeatConductionProcessData.h               |  7 ++-
 5 files changed, 84 insertions(+), 17 deletions(-)

diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index b23622d7490..38684db5536 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/Applications/ApplicationsLib/ProjectData.cpp
@@ -633,7 +633,7 @@ void ProjectData::parseProcesses(
             process = ProcessLib::HeatConduction::createHeatConductionProcess(
                 name, *_mesh_vec[0], std::move(jacobian_assembler),
                 _process_variables, _parameters, integration_order,
-                process_config);
+                process_config, _media);
         }
         else
 #endif
diff --git a/ProcessLib/HeatConduction/CreateHeatConductionProcess.cpp b/ProcessLib/HeatConduction/CreateHeatConductionProcess.cpp
index f54ea4fbfb3..35822e7f650 100644
--- a/ProcessLib/HeatConduction/CreateHeatConductionProcess.cpp
+++ b/ProcessLib/HeatConduction/CreateHeatConductionProcess.cpp
@@ -12,6 +12,9 @@
 
 #include "HeatConductionProcess.h"
 #include "HeatConductionProcessData.h"
+#include "MaterialLib/MPL/CheckMaterialSpatialDistributionMap.h"
+#include "MaterialLib/MPL/CreateMaterialSpatialDistributionMap.h"
+#include "MaterialLib/MPL/MaterialSpatialDistributionMap.h"
 #include "ParameterLib/Utils.h"
 #include "ProcessLib/Output/CreateSecondaryVariables.h"
 #include "ProcessLib/Utils/ProcessUtils.h"
@@ -20,6 +23,18 @@ namespace ProcessLib
 {
 namespace HeatConduction
 {
+void checkMPLProperties(
+    MeshLib::Mesh const& mesh,
+    MaterialPropertyLib::MaterialSpatialDistributionMap const& media_map)
+{
+    std::array const required_medium_properties = {
+        MaterialPropertyLib::PropertyType::thermal_conductivity};
+    std::array<MaterialPropertyLib::PropertyType, 0> const empty{};
+
+    MaterialPropertyLib::checkMaterialSpatialDistributionMap(
+        mesh, media_map, required_medium_properties, empty, empty);
+}
+
 std::unique_ptr<Process> createHeatConductionProcess(
     std::string name,
     MeshLib::Mesh& mesh,
@@ -27,7 +42,8 @@ std::unique_ptr<Process> createHeatConductionProcess(
     std::vector<ProcessVariable> const& variables,
     std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
     unsigned const integration_order,
-    BaseLib::ConfigTree const& config)
+    BaseLib::ConfigTree const& config,
+    std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media)
 {
     //! \ogs_file_param{prj__processes__process__type}
     config.checkConfigParameter("type", "HEAT_CONDUCTION");
@@ -47,14 +63,12 @@ std::unique_ptr<Process> createHeatConductionProcess(
          "process_variable"});
     process_variables.push_back(std::move(per_process_variables));
 
-    // thermal conductivity parameter.
-    auto& thermal_conductivity = ParameterLib::findParameter<double>(
-        config,
-        //! \ogs_file_param_special{prj__processes__process__HEAT_CONDUCTION__thermal_conductivity}
-        "thermal_conductivity", parameters, 1, &mesh);
+    auto media_map =
+        MaterialPropertyLib::createMaterialSpatialDistributionMap(media, mesh);
 
-    DBUG("Use '{:s}' as thermal conductivity parameter.",
-         thermal_conductivity.name);
+    DBUG("Check the media properties of heat conduction process ...");
+    checkMPLProperties(mesh, *media_map);
+    DBUG("Media properties verified.");
 
     // heat capacity parameter.
     auto& heat_capacity = ParameterLib::findParameter<double>(
@@ -76,7 +90,7 @@ std::unique_ptr<Process> createHeatConductionProcess(
         //! \ogs_file_param{prj__processes__process__HEAT_CONDUCTION__mass_lumping}
         config.getConfigParameter<bool>("mass_lumping", false);
 
-    HeatConductionProcessData process_data{thermal_conductivity, heat_capacity,
+    HeatConductionProcessData process_data{std::move(media_map), heat_capacity,
                                            density, mass_lumping};
 
     SecondaryVariableCollection secondary_variables;
diff --git a/ProcessLib/HeatConduction/CreateHeatConductionProcess.h b/ProcessLib/HeatConduction/CreateHeatConductionProcess.h
index 95e9fcec8d4..0b47295043d 100644
--- a/ProcessLib/HeatConduction/CreateHeatConductionProcess.h
+++ b/ProcessLib/HeatConduction/CreateHeatConductionProcess.h
@@ -13,6 +13,11 @@
 #include <memory>
 #include "ProcessLib/Process.h"
 
+namespace MaterialPropertyLib
+{
+class Medium;
+}
+
 namespace ProcessLib
 {
 namespace HeatConduction
@@ -24,7 +29,8 @@ std::unique_ptr<Process> createHeatConductionProcess(
     std::vector<ProcessVariable> const& variables,
     std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
     unsigned const integration_order,
-    BaseLib::ConfigTree const& config);
+    BaseLib::ConfigTree const& config,
+    std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media);
 
 }  // namespace HeatConduction
 }  // namespace ProcessLib
diff --git a/ProcessLib/HeatConduction/HeatConductionFEM.h b/ProcessLib/HeatConduction/HeatConductionFEM.h
index df9d2cf4186..84b4147bd3b 100644
--- a/ProcessLib/HeatConduction/HeatConductionFEM.h
+++ b/ProcessLib/HeatConduction/HeatConductionFEM.h
@@ -13,6 +13,9 @@
 #include <vector>
 
 #include "HeatConductionProcessData.h"
+#include "MaterialLib/MPL/Medium.h"
+#include "MaterialLib/MPL/Utils/FormEigenTensor.h"
+#include "MaterialLib/MPL/VariableType.h"
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "NumLib/Extrapolation/ExtrapolatableElement.h"
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
@@ -90,7 +93,7 @@ public:
         (void)local_matrix_size;
     }
 
-    void assemble(double const t, double const /*dt*/,
+    void assemble(double const t, double const dt,
                   std::vector<double> const& local_x,
                   std::vector<double> const& /*local_xdot*/,
                   std::vector<double>& local_M_data,
@@ -113,12 +116,25 @@ public:
         ParameterLib::SpatialPosition pos;
         pos.setElementID(_element.getID());
 
+        auto const& medium =
+            *_process_data.media_map->getMedium(_element.getID());
+        MaterialPropertyLib::VariableArray vars;
+        vars[static_cast<int>(MaterialPropertyLib::Variable::temperature)] =
+            medium
+                .property(
+                    MaterialPropertyLib::PropertyType::reference_temperature)
+                .template value<double>(vars, pos, t, dt);
+
         for (unsigned ip = 0; ip < n_integration_points; ip++)
         {
             pos.setIntegrationPoint(ip);
             auto const& sm = _shape_matrices[ip];
             auto const& wp = _integration_method.getWeightedPoint(ip);
-            auto const k = _process_data.thermal_conductivity(t, pos)[0];
+            auto const k = MaterialPropertyLib::formEigenTensor<GlobalDim>(
+                medium
+                    .property(
+                        MaterialPropertyLib::PropertyType::thermal_conductivity)
+                    .value(vars, pos, t, dt));
             auto const heat_capacity = _process_data.heat_capacity(t, pos)[0];
             auto const density = _process_data.density(t, pos)[0];
 
@@ -170,6 +186,15 @@ public:
         ParameterLib::SpatialPosition pos;
         pos.setElementID(_element.getID());
 
+        auto const& medium =
+            *_process_data.media_map->getMedium(_element.getID());
+        MaterialPropertyLib::VariableArray vars;
+        vars[static_cast<int>(MaterialPropertyLib::Variable::temperature)] =
+            medium
+                .property(
+                    MaterialPropertyLib::PropertyType::reference_temperature)
+                .template value<double>(vars, pos, t, dt);
+
         for (unsigned ip = 0; ip < n_integration_points; ip++)
         {
             pos.setIntegrationPoint(ip);
@@ -178,7 +203,11 @@ public:
                 _integration_method.getWeightedPoint(ip).getWeight() * sm.detJ *
                 sm.integralMeasure;
 
-            auto const k = _process_data.thermal_conductivity(t, pos)[0];
+            auto const k = MaterialPropertyLib::formEigenTensor<GlobalDim>(
+                medium
+                    .property(
+                        MaterialPropertyLib::PropertyType::thermal_conductivity)
+                    .value(vars, pos, t, dt));
             auto const heat_capacity = _process_data.heat_capacity(t, pos)[0];
             auto const density = _process_data.density(t, pos)[0];
 
@@ -196,7 +225,7 @@ public:
     }
 
     void computeSecondaryVariableConcrete(
-        double const t, double const /*dt*/, std::vector<double> const& local_x,
+        double const t, double const dt, std::vector<double> const& local_x,
         std::vector<double> const& /*local_x_dot*/) override
     {
         auto const local_matrix_size = local_x.size();
@@ -212,11 +241,24 @@ public:
         const auto local_x_vec =
             MathLib::toVector<NodalVectorType>(local_x, local_matrix_size);
 
+        auto const& medium =
+            *_process_data.media_map->getMedium(_element.getID());
+        MaterialPropertyLib::VariableArray vars;
+        vars[static_cast<int>(MaterialPropertyLib::Variable::temperature)] =
+            medium
+                .property(
+                    MaterialPropertyLib::PropertyType::reference_temperature)
+                .template value<double>(vars, pos, t, dt);
+
         for (unsigned ip = 0; ip < n_integration_points; ip++)
         {
             pos.setIntegrationPoint(ip);
             auto const& sm = _shape_matrices[ip];
-            auto const k = _process_data.thermal_conductivity(t, pos)[0];
+            auto const k = MaterialPropertyLib::formEigenTensor<GlobalDim>(
+                medium
+                    .property(
+                        MaterialPropertyLib::PropertyType::thermal_conductivity)
+                    .value(vars, pos, t, dt));
             // heat flux only computed for output.
             GlobalDimVectorType const heat_flux = -k * sm.dNdx * local_x_vec;
 
diff --git a/ProcessLib/HeatConduction/HeatConductionProcessData.h b/ProcessLib/HeatConduction/HeatConductionProcessData.h
index 0801fdf0822..e201c923d7b 100644
--- a/ProcessLib/HeatConduction/HeatConductionProcessData.h
+++ b/ProcessLib/HeatConduction/HeatConductionProcessData.h
@@ -10,6 +10,9 @@
 
 #pragma once
 
+#include <memory>
+#include "MaterialLib/MPL/MaterialSpatialDistributionMap.h"
+
 namespace ParameterLib
 {
 template <typename T>
@@ -20,7 +23,9 @@ namespace ProcessLib::HeatConduction
 {
 struct HeatConductionProcessData
 {
-    ParameterLib::Parameter<double> const& thermal_conductivity;
+    std::unique_ptr<MaterialPropertyLib::MaterialSpatialDistributionMap>
+        media_map;
+
     ParameterLib::Parameter<double> const& heat_capacity;
     ParameterLib::Parameter<double> const& density;
 
-- 
GitLab