From 51777770b562f1fdca9930fb7202e6fda0193c17 Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <dmitri.naumov@ufz.de>
Date: Mon, 2 Dec 2019 16:14:44 +0100
Subject: [PATCH] [PL] HM Add MFront internal state variables output

---
 .../HydroMechanics/HydroMechanicsFEM-impl.h   |  20 +++
 ProcessLib/HydroMechanics/HydroMechanicsFEM.h |   9 +-
 .../HydroMechanics/HydroMechanicsProcess.cpp  | 114 ++++++++----------
 .../HydroMechanics/HydroMechanicsProcess.h    |   4 +-
 .../HydroMechanics/LocalAssemblerInterface.h  |  14 ++-
 5 files changed, 94 insertions(+), 67 deletions(-)

diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
index 00177dac076..8e9e177562d 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
@@ -796,5 +796,25 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
     }
 }
 
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+          typename IntegrationMethod, int DisplacementDim>
+unsigned HydroMechanicsLocalAssembler<
+    ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
+    DisplacementDim>::getNumberOfIntegrationPoints() const
+{
+    return _integration_method.getNumberOfPoints();
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+          typename IntegrationMethod, int DisplacementDim>
+typename MaterialLib::Solids::MechanicsBase<
+    DisplacementDim>::MaterialStateVariables const&
+HydroMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
+                             IntegrationMethod, DisplacementDim>::
+    getMaterialStateVariablesAt(unsigned integration_point) const
+{
+    return *_ip_data[integration_point].material_state_variables;
+}
+
 }  // namespace HydroMechanics
 }  // namespace ProcessLib
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
index c3cfe2fade5..217a7689ef4 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
@@ -106,7 +106,8 @@ struct SecondaryData
 
 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
           typename IntegrationMethod, int DisplacementDim>
-class HydroMechanicsLocalAssembler : public LocalAssemblerInterface
+class HydroMechanicsLocalAssembler
+    : public LocalAssemblerInterface<DisplacementDim>
 {
 public:
     using ShapeMatricesTypeDisplacement =
@@ -436,6 +437,12 @@ private:
         std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
         LocalCoupledSolutions const& local_coupled_solutions);
 
+    unsigned getNumberOfIntegrationPoints() const override;
+
+    typename MaterialLib::Solids::MechanicsBase<
+        DisplacementDim>::MaterialStateVariables const&
+    getMaterialStateVariablesAt(unsigned integration_point) const override;
+
 private:
     HydroMechanicsProcessData<DisplacementDim>& _process_data;
 
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsProcess.cpp b/ProcessLib/HydroMechanics/HydroMechanicsProcess.cpp
index a4557810dc9..96467d80c85 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsProcess.cpp
+++ b/ProcessLib/HydroMechanics/HydroMechanicsProcess.cpp
@@ -12,14 +12,13 @@
 
 #include <cassert>
 
+#include "HydroMechanicsFEM.h"
+#include "HydroMechanicsProcessData.h"
 #include "MeshLib/Elements/Utils.h"
 #include "NumLib/DOF/ComputeSparsityPattern.h"
+#include "ProcessLib/Deformation/SolidMaterialInternalToSecondaryVariables.h"
 #include "ProcessLib/HydroMechanics/CreateLocalAssemblers.h"
 #include "ProcessLib/Process.h"
-
-#include "HydroMechanicsFEM.h"
-#include "HydroMechanicsProcessData.h"
-
 namespace ProcessLib
 {
 namespace HydroMechanics
@@ -208,77 +207,68 @@ void HydroMechanicsProcess<DisplacementDim>::initializeConcreteProcess(
         _local_assemblers, mesh.isAxiallySymmetric(), integration_order,
         _process_data);
 
-    _secondary_variables.addSecondaryVariable(
-        "sigma_xx",
-        makeExtrapolator(1, getExtrapolator(), _local_assemblers,
-                         &LocalAssemblerInterface::getIntPtSigmaXX));
+    auto add_secondary_variable = [&](std::string const& name,
+                                      int const num_components,
+                                      auto get_ip_values_function) {
+        _secondary_variables.addSecondaryVariable(
+            name,
+            makeExtrapolator(num_components, getExtrapolator(),
+                             _local_assemblers,
+                             std::move(get_ip_values_function)));
+    };
+
+    add_secondary_variable(
+        "sigma_xx", 1, &LocalAssemblerIF::getIntPtSigmaXX);
 
-    _secondary_variables.addSecondaryVariable(
-        "sigma_yy",
-        makeExtrapolator(1, getExtrapolator(), _local_assemblers,
-                         &LocalAssemblerInterface::getIntPtSigmaYY));
+    add_secondary_variable(
+        "sigma_yy", 1, &LocalAssemblerIF::getIntPtSigmaYY);
 
-    _secondary_variables.addSecondaryVariable(
-        "sigma_zz",
-        makeExtrapolator(1, getExtrapolator(), _local_assemblers,
-                         &LocalAssemblerInterface::getIntPtSigmaZZ));
+    add_secondary_variable(
+        "sigma_zz", 1, &LocalAssemblerIF::getIntPtSigmaZZ);
 
-    _secondary_variables.addSecondaryVariable(
-        "sigma_xy",
-        makeExtrapolator(1, getExtrapolator(), _local_assemblers,
-                         &LocalAssemblerInterface::getIntPtSigmaXY));
+    add_secondary_variable(
+        "sigma_xy", 1, &LocalAssemblerIF::getIntPtSigmaXY);
 
     if (DisplacementDim == 3)
     {
-        _secondary_variables.addSecondaryVariable(
-            "sigma_xz",
-            makeExtrapolator(1, getExtrapolator(), _local_assemblers,
-                             &LocalAssemblerInterface::getIntPtSigmaXZ));
+        add_secondary_variable(
+            "sigma_xz", 1, &LocalAssemblerIF::getIntPtSigmaXZ);
 
-        _secondary_variables.addSecondaryVariable(
-            "sigma_yz",
-            makeExtrapolator(1, getExtrapolator(), _local_assemblers,
-                             &LocalAssemblerInterface::getIntPtSigmaYZ));
+        add_secondary_variable(
+            "sigma_yz", 1, &LocalAssemblerIF::getIntPtSigmaYZ);
     }
 
-    _secondary_variables.addSecondaryVariable(
-        "epsilon_xx",
-        makeExtrapolator(1, getExtrapolator(), _local_assemblers,
-                         &LocalAssemblerInterface::getIntPtEpsilonXX));
+    add_secondary_variable(
+        "epsilon_xx", 1, &LocalAssemblerIF::getIntPtEpsilonXX);
 
-    _secondary_variables.addSecondaryVariable(
-        "epsilon_yy",
-        makeExtrapolator(1, getExtrapolator(), _local_assemblers,
-                         &LocalAssemblerInterface::getIntPtEpsilonYY));
+    add_secondary_variable(
+        "epsilon_yy", 1, &LocalAssemblerIF::getIntPtEpsilonYY);
 
-    _secondary_variables.addSecondaryVariable(
-        "epsilon_zz",
-        makeExtrapolator(1, getExtrapolator(), _local_assemblers,
-                         &LocalAssemblerInterface::getIntPtEpsilonZZ));
+    add_secondary_variable(
+        "epsilon_zz", 1, &LocalAssemblerIF::getIntPtEpsilonZZ);
 
-    _secondary_variables.addSecondaryVariable(
-        "epsilon_xy",
-        makeExtrapolator(1, getExtrapolator(), _local_assemblers,
-                         &LocalAssemblerInterface::getIntPtEpsilonXY));
+    add_secondary_variable(
+        "epsilon_xy", 1, &LocalAssemblerIF::getIntPtEpsilonXY);
 
     if (DisplacementDim == 3)
     {
-        _secondary_variables.addSecondaryVariable(
-            "epsilon_xz",
-            makeExtrapolator(1, getExtrapolator(), _local_assemblers,
-                             &LocalAssemblerInterface::getIntPtEpsilonXZ));
+        add_secondary_variable(
+            "epsilon_xz", 1, &LocalAssemblerIF::getIntPtEpsilonXZ);
 
-        _secondary_variables.addSecondaryVariable(
-            "epsilon_yz",
-            makeExtrapolator(1, getExtrapolator(), _local_assemblers,
-                             &LocalAssemblerInterface::getIntPtEpsilonYZ));
+        add_secondary_variable(
+            "epsilon_yz", 1, &LocalAssemblerIF::getIntPtEpsilonYZ);
     }
 
-    _secondary_variables.addSecondaryVariable(
-        "velocity",
-        makeExtrapolator(mesh.getDimension(), getExtrapolator(),
-                         _local_assemblers,
-                         &LocalAssemblerInterface::getIntPtDarcyVelocity));
+    add_secondary_variable("velocity",
+                           DisplacementDim,
+                           &LocalAssemblerIF::getIntPtDarcyVelocity);
+
+    //
+    // enable output of internal variables defined by material models
+    //
+    ProcessLib::Deformation::solidMaterialInternalToSecondaryVariables<
+        LocalAssemblerIF>(_process_data.solid_materials,
+                                 add_secondary_variable);
 
     _process_data.pressure_interpolated =
         MeshLib::getOrCreateMeshProperty<double>(
@@ -357,7 +347,7 @@ void HydroMechanicsProcess<DisplacementDim>::initializeConcreteProcess(
 
     // Initialize local assemblers after all variables have been set.
     GlobalExecutor::executeMemberOnDereferenced(
-        &LocalAssemblerInterface::initialize, _local_assemblers,
+        &LocalAssemblerIF::initialize, _local_assemblers,
         *_local_to_global_index_map);
 }
 
@@ -483,7 +473,7 @@ void HydroMechanicsProcess<DisplacementDim>::preTimestepConcreteProcess(
         ProcessLib::ProcessVariable const& pv =
             getProcessVariables(process_id)[0];
         GlobalExecutor::executeSelectedMemberOnDereferenced(
-            &LocalAssemblerInterface::preTimestep, _local_assemblers,
+            &LocalAssemblerIF::preTimestep, _local_assemblers,
             pv.getActiveElementIDs(), *_local_to_global_index_map,
             *x[process_id], t, dt);
     }
@@ -497,7 +487,7 @@ void HydroMechanicsProcess<DisplacementDim>::postTimestepConcreteProcess(
     DBUG("PostTimestep HydroMechanicsProcess.");
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
     GlobalExecutor::executeSelectedMemberOnDereferenced(
-        &LocalAssemblerInterface::postTimestep, _local_assemblers,
+        &LocalAssemblerIF::postTimestep, _local_assemblers,
         pv.getActiveElementIDs(), getDOFTable(process_id), *x[process_id], t,
         dt);
 }
@@ -516,7 +506,7 @@ void HydroMechanicsProcess<DisplacementDim>::postNonLinearSolverConcreteProcess(
     // Calculate strain, stress or other internal variables of mechanics.
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
     GlobalExecutor::executeSelectedMemberOnDereferenced(
-        &LocalAssemblerInterface::postNonLinearSolver, _local_assemblers,
+        &LocalAssemblerIF::postNonLinearSolver, _local_assemblers,
         pv.getActiveElementIDs(), getDOFTable(process_id), x, t, dt,
         _use_monolithic_scheme);
 }
@@ -528,7 +518,7 @@ void HydroMechanicsProcess<DisplacementDim>::computeSecondaryVariableConcrete(
     DBUG("Compute the secondary variables for HydroMechanicsProcess.");
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
     GlobalExecutor::executeSelectedMemberOnDereferenced(
-        &LocalAssemblerInterface::computeSecondaryVariable, _local_assemblers,
+        &LocalAssemblerIF::computeSecondaryVariable, _local_assemblers,
         pv.getActiveElementIDs(), getDOFTable(process_id), t, x,
         _coupled_solutions);
 }
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsProcess.h b/ProcessLib/HydroMechanics/HydroMechanicsProcess.h
index 5768ae43870..d1e254334bc 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsProcess.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsProcess.h
@@ -62,6 +62,8 @@ public:
         const int process_id) const override;
 
 private:
+    using LocalAssemblerIF = LocalAssemblerInterface<DisplacementDim>;
+
     void constructDofTable() override;
 
     void initializeConcreteProcess(
@@ -103,7 +105,7 @@ private:
     std::unique_ptr<MeshLib::MeshSubset const> _mesh_subset_base_nodes;
     HydroMechanicsProcessData<DisplacementDim> _process_data;
 
-    std::vector<std::unique_ptr<LocalAssemblerInterface>> _local_assemblers;
+    std::vector<std::unique_ptr<LocalAssemblerIF>> _local_assemblers;
 
     std::unique_ptr<NumLib::LocalToGlobalIndexMap>
         _local_to_global_index_map_single_component;
diff --git a/ProcessLib/HydroMechanics/LocalAssemblerInterface.h b/ProcessLib/HydroMechanics/LocalAssemblerInterface.h
index c9f999626b5..7781266d0fb 100644
--- a/ProcessLib/HydroMechanics/LocalAssemblerInterface.h
+++ b/ProcessLib/HydroMechanics/LocalAssemblerInterface.h
@@ -10,6 +10,7 @@
 
 #pragma once
 
+#include "MaterialLib/SolidModels/MechanicsBase.h"
 #include "NumLib/Extrapolation/ExtrapolatableElement.h"
 #include "ProcessLib/LocalAssemblerInterface.h"
 
@@ -17,9 +18,9 @@ namespace ProcessLib
 {
 namespace HydroMechanics
 {
-struct LocalAssemblerInterface
-    : public ProcessLib::LocalAssemblerInterface,
-      public NumLib::ExtrapolatableElement
+template <int DisplacementDim>
+struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
+                                 public NumLib::ExtrapolatableElement
 {
     virtual std::size_t setIPDataInitialConditions(
         std::string const& name, double const* values,
@@ -106,6 +107,13 @@ struct LocalAssemblerInterface
         std::vector<GlobalVector*> const& x,
         std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
         std::vector<double>& cache) const = 0;
+
+    // TODO move to NumLib::ExtrapolatableElement
+    virtual unsigned getNumberOfIntegrationPoints() const = 0;
+
+    virtual typename MaterialLib::Solids::MechanicsBase<
+        DisplacementDim>::MaterialStateVariables const&
+    getMaterialStateVariablesAt(unsigned /*integration_point*/) const = 0;
 };
 
 }  // namespace HydroMechanics
-- 
GitLab