From 1be4290d51c0f0bb6ab54b30efd0e1e73a921e24 Mon Sep 17 00:00:00 2001
From: Christoph Lehmann <christoph.lehmann@ufz.de>
Date: Tue, 26 Jul 2022 14:34:04 +0200
Subject: [PATCH] [PL] Moved porosity initialization

---
 .../LocalAssemblerInterface.h                 | 38 +++++++++++++++++++
 .../ThermoRichardsMechanicsFEM-impl.h         | 26 -------------
 2 files changed, 38 insertions(+), 26 deletions(-)

diff --git a/ProcessLib/ThermoRichardsMechanics/LocalAssemblerInterface.h b/ProcessLib/ThermoRichardsMechanics/LocalAssemblerInterface.h
index 95fc8994c39..c2971c7081c 100644
--- a/ProcessLib/ThermoRichardsMechanics/LocalAssemblerInterface.h
+++ b/ProcessLib/ThermoRichardsMechanics/LocalAssemblerInterface.h
@@ -50,6 +50,44 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
         {
             material_states_.emplace_back(solid_material_);
         }
+
+        ParameterLib::SpatialPosition x_position;
+        x_position.setElementID(e.getID());
+
+        auto const& medium = process_data_.media_map->getMedium(e.getID());
+
+        for (unsigned ip = 0; ip < n_integration_points; ip++)
+        {
+            namespace MPL = MaterialPropertyLib;
+
+            x_position.setIntegrationPoint(ip);
+
+            auto& current_state = current_states_[ip];
+
+            // Initial porosity. Could be read from integration point data or
+            // mesh.
+            current_state.poro_data.phi =
+                medium->property(MPL::porosity)
+                    .template initialValue<double>(
+                        x_position,
+                        std::numeric_limits<
+                            double>::quiet_NaN() /* t independent */);
+
+            if (medium->hasProperty(MPL::PropertyType::transport_porosity))
+            {
+                current_state.transport_poro_data.phi =
+                    medium->property(MPL::transport_porosity)
+                        .template initialValue<double>(
+                            x_position,
+                            std::numeric_limits<
+                                double>::quiet_NaN() /* t independent */);
+            }
+            else
+            {
+                current_state.transport_poro_data.phi =
+                    current_state.poro_data.phi;
+            }
+        }
     }
 
     std::size_t setIPDataInitialConditions(std::string const& name,
diff --git a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h
index c2bcedfe026..19b54a300c2 100644
--- a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h
+++ b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h
@@ -62,8 +62,6 @@ ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunction,
                                   DisplacementDim>(e, is_axially_symmetric,
                                                    integration_method);
 
-    auto const& medium = this->process_data_.media_map->getMedium(e.getID());
-
     ParameterLib::SpatialPosition x_position;
     x_position.setElementID(e.getID());
     for (unsigned ip = 0; ip < n_integration_points; ip++)
@@ -93,30 +91,6 @@ ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunction,
         ip_data.N_p = shape_matrices[ip].N;
         ip_data.dNdx_p = shape_matrices[ip].dNdx;
 
-        auto& current_state = this->current_states_[ip];
-
-        // Initial porosity. Could be read from integration point data or mesh.
-        current_state.poro_data.phi =
-            medium->property(MPL::porosity)
-                .template initialValue<double>(
-                    x_position,
-                    std::numeric_limits<
-                        double>::quiet_NaN() /* t independent */);
-
-        if (medium->hasProperty(MPL::PropertyType::transport_porosity))
-        {
-            current_state.transport_poro_data.phi =
-                medium->property(MPL::transport_porosity)
-                    .template initialValue<double>(
-                        x_position,
-                        std::numeric_limits<
-                            double>::quiet_NaN() /* t independent */);
-        }
-        else
-        {
-            current_state.transport_poro_data.phi = current_state.poro_data.phi;
-        }
-
         secondary_data_.N_u[ip] = shape_matrices_u[ip].N;
     }
 }
-- 
GitLab