diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
index 5102fd192d5409888b3738016e3abccf9ce5f335..3bc2cebedc2d2f21da60573cd775a0810cc59af8 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 21ac46ffd983eb7906ba6ef5e8053eaeaaa85f4d..1bb61385358e7eb347fa127010fc2143613934eb 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 916e952f4cd560f5c55c64c0b0c7f90503d31dee..3da2948367554e7a3318bd3caf0ceed827e88c53 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;
 };