[PL/TH2M] Extract constitutive variables.

A class of variables similar to those stored in the
integration point data but needed only temporary for
the assembly process.

ProcessLib/TH2M/ConstitutiveVariables.h

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ */
+
+#pragma once
+
+#include "MathLib/KelvinVector.h"
+
+namespace ProcessLib::TH2M
+{
+/// Variables needed only for the assembly process. The values are not preserved
+/// throughout the iterations contrary to the variables in IntegrationPointData.
+template <int DisplacementDim>
+struct ConstitutiveVariables
+{
+    using KelvinMatrixType =
+        MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>;
+
+    KelvinMatrixType C;
+
+    EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
+};
+
+}  // namespace ProcessLib::TH2M

ProcessLib/TH2M/TH2MFEM-impl.h

@@ -93,8 +93,9 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
 
 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
           typename IntegrationMethod, int DisplacementDim>
-void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
-                        IntegrationMethod, DisplacementDim>::
+std::vector<ConstitutiveVariables<DisplacementDim>>
+TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
+                   IntegrationMethod, DisplacementDim>::
     updateConstitutiveVariables(Eigen::VectorXd const& local_x,
                                 Eigen::VectorXd const& local_x_dot,
                                 double const t, double const dt)
@@ -127,10 +128,15 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
 
     unsigned const n_integration_points =
         _integration_method.getNumberOfPoints();
+
+    std::vector<ConstitutiveVariables<DisplacementDim>>
+        ip_constitutive_variables(n_integration_points);
+
     for (unsigned ip = 0; ip < n_integration_points; ip++)
     {
         pos.setIntegrationPoint(ip);
         auto& ip_data = _ip_data[ip];
+        auto& ip_cv = ip_constitutive_variables[ip];
 
         auto const& Np = ip_data.N_p;
         auto const& NT = Np;
@@ -260,7 +266,8 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
         auto const rhoSR = rho_ref_SR;
 #endif  // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
 
-        ip_data.updateConstitutiveRelation(vars, t, pos, dt, T - T_dot * dt);
+        ip_cv.C = ip_data.updateConstitutiveRelation(vars, t, pos, dt,
+                                                     T - T_dot * dt);
 
         // constitutive model object as specified in process creation
         auto& ptm = *_process_data.phase_transition_model_;
@@ -330,6 +337,8 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
         ip_data.h_L = c.hL;
         ip_data.pWGR = c.pWGR;
     }
+
+    return ip_constitutive_variables;
 }
 
 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,

ProcessLib/TH2M/TH2MFEM.h

@@ -13,6 +13,7 @@
 #include <memory>
 #include <vector>
 
+#include "ConstitutiveVariables.h"
 #include "IntegrationPointData.h"
 #include "LocalAssemblerInterface.h"
 #include "MaterialLib/PhysicalConstant.h"
@@ -169,9 +170,10 @@ private:
         std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
         std::vector<double>& cache) const override;
 
-    void updateConstitutiveVariables(Eigen::VectorXd const& local_x,
-                                     Eigen::VectorXd const& local_x_dot,
-                                     double const t, double const dt);
+    std::vector<ConstitutiveVariables<DisplacementDim>>
+    updateConstitutiveVariables(Eigen::VectorXd const& local_x,
+                                Eigen::VectorXd const& local_x_dot,
+                                double const t, double const dt);
 
     std::size_t setSigma(double const* values)
     {