[TM] support epsilon in restart

ufz · Jan 17, 2019 · 164f3f3 · 164f3f3
1 parent 8972d34
commit 164f3f3
Show file tree

Hide file tree

Showing 4 changed files with 79 additions and 2 deletions.
diff --git a/ProcessLib/ThermoMechanics/LocalAssemblerInterface.h b/ProcessLib/ThermoMechanics/LocalAssemblerInterface.h
@@ -28,6 +28,8 @@ struct ThermoMechanicsLocalAssemblerInterface
 
     virtual std::vector<double> getSigma() const = 0;
 
+    virtual std::vector<double> getEpsilon() const = 0;
+
     virtual std::vector<double> const& getIntPtSigma(
         const double /*t*/,
         GlobalVector const& /*current_solution*/,

diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h b/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h
@@ -186,6 +186,10 @@ class ThermoMechanicsLocalAssembler
         {
             return setSigma(values);
         }
+        else if (name == "epsilon_ip")
+        {
+            return setEpsilon(values);
+        }
 
         return 0;
     }
@@ -484,6 +488,53 @@ class ThermoMechanicsLocalAssembler
         return cache;
     }
 
+    std::size_t setEpsilon(double const* values)
+    {
+        auto const kelvin_vector_size =
+            MathLib::KelvinVector::KelvinVectorDimensions<
+                DisplacementDim>::value;
+        unsigned const n_integration_points =
+            _integration_method.getNumberOfPoints();
+
+        std::vector<double> ip_epsilon_values;
+        auto epsilon_values =
+            Eigen::Map<Eigen::Matrix<double, kelvin_vector_size, Eigen::Dynamic,
+                                     Eigen::ColMajor> const>(
+                values, kelvin_vector_size, n_integration_points);
+
+        for (unsigned ip = 0; ip < n_integration_points; ++ip)
+        {
+            _ip_data[ip].eps =
+                MathLib::KelvinVector::symmetricTensorToKelvinVector(
+                    epsilon_values.col(ip));
+        }
+
+        return n_integration_points;
+    }
+
+    std::vector<double> getEpsilon() const override
+    {
+        auto const kelvin_vector_size =
+            MathLib::KelvinVector::KelvinVectorDimensions<
+                DisplacementDim>::value;
+        unsigned const n_integration_points =
+            _integration_method.getNumberOfPoints();
+
+        std::vector<double> ip_epsilon_values;
+        auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
+            double, Eigen::Dynamic, kelvin_vector_size, Eigen::RowMajor>>(
+            ip_epsilon_values, n_integration_points, kelvin_vector_size);
+
+        for (unsigned ip = 0; ip < n_integration_points; ++ip)
+        {
+            auto const& eps = _ip_data[ip].eps;
+            cache_mat.row(ip) =
+                MathLib::KelvinVector::kelvinVectorToSymmetricTensor(eps);
+        }
+
+        return ip_epsilon_values;
+    }
+
     virtual std::vector<double> const& getIntPtEpsilon(
         const double /*t*/,
         GlobalVector const& /*current_solution*/,

diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.cpp b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.cpp
@@ -40,6 +40,7 @@ ThermoMechanicsProcess<DisplacementDim>::ThermoMechanicsProcess(
 {
     _integration_point_writer.emplace_back(
         std::make_unique<SigmaIntegrationPointWriter>(
+            "sigma_ip",
             static_cast<int>(mesh.getDimension() == 2 ? 4 : 6) /*n components*/,
             2 /*integration order*/, [this]() {
                 // Result containing integration point data for each local
@@ -54,6 +55,26 @@ ThermoMechanicsProcess<DisplacementDim>::ThermoMechanicsProcess(
                     result[i] = local_asm.getSigma();
                 }
 
+                return result;
+            }));
+
+    _integration_point_writer.emplace_back(
+        std::make_unique<SigmaIntegrationPointWriter>(
+            "epsilon_ip",
+            static_cast<int>(mesh.getDimension() == 2 ? 4 : 6) /*n components*/,
+            2 /*integration order*/, [this]() {
+                // Result containing integration point data for each local
+                // assembler.
+                std::vector<std::vector<double>> result;
+                result.resize(_local_assemblers.size());
+
+                for (std::size_t i = 0; i < _local_assemblers.size(); ++i)
+                {
+                    auto const& local_asm = *_local_assemblers[i];
+
+                    result[i] = local_asm.getEpsilon();
+                }
+
                 return result;
             }));
 }

diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.h b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.h
@@ -23,11 +23,13 @@ struct ThermoMechanicsLocalAssemblerInterface;
 struct SigmaIntegrationPointWriter final : public IntegrationPointWriter
 {
     explicit SigmaIntegrationPointWriter(
+        std::string const& name,
         int const n_components,
         int const integration_order,
         std::function<std::vector<std::vector<double>>()>
             callback)
-        : _n_components(n_components),
+        : _name(name),
+          _n_components(n_components),
           _integration_order(integration_order),
           _callback(callback)
     {
@@ -41,7 +43,7 @@ struct SigmaIntegrationPointWriter final : public IntegrationPointWriter
         // TODO (naumov) remove ip suffix. Probably needs modification of the
         // mesh properties, s.t. there is no "overlapping" with cell/point data.
         // See getOrCreateMeshProperty.
-        return "sigma_ip";
+        return _name;
     }
 
     std::vector<std::vector<double>> values() const override
@@ -50,6 +52,7 @@ struct SigmaIntegrationPointWriter final : public IntegrationPointWriter
     }
 
 private:
+    std::string const _name;
     int const _n_components;
     int const _integration_order;
     std::function<std::vector<std::vector<double>>()> _callback;