diff --git a/ProcessLib/ThermoMechanics/LocalAssemblerInterface.h b/ProcessLib/ThermoMechanics/LocalAssemblerInterface.h
index 23fcaa4b9d0aba263bb812ed897d08f4a752c6be..22a541c8cdcaa2a872f0518ee330e55a85971905 100644
--- a/ProcessLib/ThermoMechanics/LocalAssemblerInterface.h
+++ b/ProcessLib/ThermoMechanics/LocalAssemblerInterface.h
@@ -28,6 +28,10 @@ struct ThermoMechanicsLocalAssemblerInterface
virtual std::vector<double> getSigma() const = 0;
+ virtual std::vector<double> getEpsilon() const = 0;
+
+ virtual std::vector<double> getEpsilonMechanical() 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
index ffdaf39ba3ad6b2ee1e506de08cd051b19c63327..82080cbc3767d999ca8e66caf8bbfe0eae1cbf0b 100644
--- a/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h
@@ -186,6 +186,14 @@ public:
{
return setSigma(values);
}
+ else if (name == "epsilon_ip")
+ {
+ return setEpsilon(values);
+ }
+ else if (name == "epsilon_m_ip")
+ {
+ return setEpsilonMechanical(values);
+ }
return 0;
}
@@ -415,7 +423,6 @@ private:
unsigned const n_integration_points =
_integration_method.getNumberOfPoints();
- std::vector<double> ip_sigma_values;
auto sigma_values =
Eigen::Map<Eigen::Matrix<double, kelvin_vector_size, Eigen::Dynamic,
Eigen::ColMajor> const>(
@@ -484,6 +491,52 @@ private:
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();
+
+ 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*/,
@@ -511,6 +564,52 @@ private:
return cache;
}
+ std::size_t setEpsilonMechanical(double const* values)
+ {
+ auto const kelvin_vector_size =
+ MathLib::KelvinVector::KelvinVectorDimensions<
+ DisplacementDim>::value;
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ auto epsilon_m_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_m =
+ MathLib::KelvinVector::symmetricTensorToKelvinVector(
+ epsilon_m_values.col(ip));
+ }
+
+ return n_integration_points;
+ }
+
+ std::vector<double> getEpsilonMechanical() 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_m_values;
+ auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
+ double, Eigen::Dynamic, kelvin_vector_size, Eigen::RowMajor>>(
+ ip_epsilon_m_values, n_integration_points, kelvin_vector_size);
+
+ for (unsigned ip = 0; ip < n_integration_points; ++ip)
+ {
+ auto const& eps_m = _ip_data[ip].eps_m;
+ cache_mat.row(ip) =
+ MathLib::KelvinVector::kelvinVectorToSymmetricTensor(eps_m);
+ }
+
+ return ip_epsilon_m_values;
+ }
+
ThermoMechanicsProcessData<DisplacementDim>& _process_data;
std::vector<
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.cpp b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.cpp
index 042783e58534d37ec9e7ac97c9f6a85c363e83c2..c819806691ca022883bd20a5c9140e839e4bb153 100644
--- a/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.cpp
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.cpp
@@ -39,7 +39,8 @@ ThermoMechanicsProcess<DisplacementDim>::ThermoMechanicsProcess(
_process_data(std::move(process_data))
{
_integration_point_writer.emplace_back(
- std::make_unique<SigmaIntegrationPointWriter>(
+ std::make_unique<KelvinVectorIntegrationPointWriter>(
+ "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,46 @@ ThermoMechanicsProcess<DisplacementDim>::ThermoMechanicsProcess(
result[i] = local_asm.getSigma();
}
+ return result;
+ }));
+
+ _integration_point_writer.emplace_back(
+ std::make_unique<KelvinVectorIntegrationPointWriter>(
+ "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;
+ }));
+
+ _integration_point_writer.emplace_back(
+ std::make_unique<KelvinVectorIntegrationPointWriter>(
+ "epsilon_m_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.getEpsilonMechanical();
+ }
+
return result;
}));
}
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.h b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.h
index b124b068a473b877697f9be83a7bc40b6e0be095..869e7e7e41cc8ed2c880e481448d067bdded7d13 100644
--- a/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.h
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.h
@@ -20,14 +20,16 @@ namespace ThermoMechanics
{
struct ThermoMechanicsLocalAssemblerInterface;
-struct SigmaIntegrationPointWriter final : public IntegrationPointWriter
+struct KelvinVectorIntegrationPointWriter final : public IntegrationPointWriter
{
- explicit SigmaIntegrationPointWriter(
+ explicit KelvinVectorIntegrationPointWriter(
+ 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;