From 1036163a0bf56835abe5d841a1a31f4886e76fe4 Mon Sep 17 00:00:00 2001
From: Florian Zill <florian.zill@ufz.de>
Date: Tue, 15 Oct 2019 14:25:06 +0200
Subject: [PATCH] [HM] Output Principal Stress and Vectors
Averaged per Element
Stored as MeshLib::Property Vectors
-> saved by default
---
.../HydroMechanics/HydroMechanicsFEM-impl.h | 29 +++++++++++++++++++
.../HydroMechanics/HydroMechanicsProcess.cpp | 20 +++++++++++++
.../HydroMechanicsProcessData.h | 3 ++
3 files changed, 52 insertions(+)
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
index 5102fd192d5..3bc2cebedc2 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
@@ -764,6 +764,35 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
DisplacementDim>(_element, _is_axially_symmetric, p,
*_process_data.pressure_interpolated);
+
+ auto sigma_eff_sum = MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
+ Eigen::Matrix<double, 3, 3>::Zero());
+
+ for (auto const& ip_data : _ip_data)
+ {
+ sigma_eff_sum += ip_data.sigma_eff;
+ }
+
+ Eigen::Matrix<double, 3, 3, 0, 3, 3> sigma_avg =
+ MathLib::KelvinVector::kelvinVectorToTensor(sigma_eff_sum) /
+ _integration_method.getNumberOfPoints();
+
+ Eigen::SelfAdjointEigenSolver<Eigen::Matrix<double, 3, 3>> e_s(sigma_avg);
+
+ int const elem_id = _element.getID();
+
+ Eigen::Map<Eigen::Vector3d>(
+ &(*_process_data.principal_stress_values)[elem_id * 3], 3) =
+ e_s.eigenvalues();
+
+ auto eigen_vectors = e_s.eigenvectors();
+
+ for (auto i = 0; i < 3; i++)
+ {
+ Eigen::Map<Eigen::Vector3d>(
+ &(*_process_data.principal_stress_vector[i])[elem_id * 3], 3) =
+ eigen_vectors.col(i);
+ }
}
} // namespace HydroMechanics
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsProcess.cpp b/ProcessLib/HydroMechanics/HydroMechanicsProcess.cpp
index 21ac46ffd98..1bb61385358 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsProcess.cpp
+++ b/ProcessLib/HydroMechanics/HydroMechanicsProcess.cpp
@@ -287,6 +287,26 @@ void HydroMechanicsProcess<DisplacementDim>::initializeConcreteProcess(
const_cast<MeshLib::Mesh&>(mesh), "pressure_interpolated",
MeshLib::MeshItemType::Node, 1);
+ _process_data.principal_stress_vector[0] =
+ MeshLib::getOrCreateMeshProperty<double>(
+ const_cast<MeshLib::Mesh&>(mesh), "principal_stress_vector_1",
+ MeshLib::MeshItemType::Cell, 3);
+
+ _process_data.principal_stress_vector[1] =
+ MeshLib::getOrCreateMeshProperty<double>(
+ const_cast<MeshLib::Mesh&>(mesh), "principal_stress_vector_2",
+ MeshLib::MeshItemType::Cell, 3);
+
+ _process_data.principal_stress_vector[2] =
+ MeshLib::getOrCreateMeshProperty<double>(
+ const_cast<MeshLib::Mesh&>(mesh), "principal_stress_vector_3",
+ MeshLib::MeshItemType::Cell, 3);
+
+ _process_data.principal_stress_values =
+ MeshLib::getOrCreateMeshProperty<double>(
+ const_cast<MeshLib::Mesh&>(mesh), "principal_stress_values",
+ MeshLib::MeshItemType::Cell, 3);
+
// Set initial conditions for integration point data.
for (auto const& ip_writer : _integration_point_writer)
{
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsProcessData.h b/ProcessLib/HydroMechanics/HydroMechanicsProcessData.h
index 916e952f4cd..3da29483675 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsProcessData.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsProcessData.h
@@ -120,6 +120,9 @@ struct HydroMechanicsProcessData
}
MeshLib::PropertyVector<double>* pressure_interpolated = nullptr;
+ std::array<MeshLib::PropertyVector<double>*, 3> principal_stress_vector = {
+ nullptr, nullptr, nullptr};
+ MeshLib::PropertyVector<double>* principal_stress_values = nullptr;
EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
};
--
GitLab