From c462da9cd6b651108287da36a74b0ee1f6eb9a70 Mon Sep 17 00:00:00 2001
From: Norbert Grunwald <n.g.boettcher@gmail.com>
Date: Fri, 12 Jan 2018 15:05:32 +0100
Subject: [PATCH] changed deviator trace assertion

Used diagonal of tensor and fixed norm
---
 MathLib/KelvinVector-impl.h | 15 +++++++++++----
 MathLib/KelvinVector.h      |  5 +++++
 2 files changed, 16 insertions(+), 4 deletions(-)

diff --git a/MathLib/KelvinVector-impl.h b/MathLib/KelvinVector-impl.h
index 8916fba496c..aabed08116a 100644
--- a/MathLib/KelvinVector-impl.h
+++ b/MathLib/KelvinVector-impl.h
@@ -16,7 +16,7 @@ double Invariants<KelvinVectorSize>::equivalentStress(
     Eigen::Matrix<double, KelvinVectorSize, 1> const& deviatoric_v)
 {
     assert(std::abs(trace(deviatoric_v)) <=
-           1e-16 * std::abs(deviatoric_v.squaredNorm()));
+           5e-14 * diagonal(deviatoric_v).norm());
     return std::sqrt(3 * J2(deviatoric_v));
 }
 
@@ -25,7 +25,7 @@ double Invariants<KelvinVectorSize>::J2(
     Eigen::Matrix<double, KelvinVectorSize, 1> const& deviatoric_v)
 {
     assert(std::abs(trace(deviatoric_v)) <=
-           1e-16 * std::abs(deviatoric_v.squaredNorm()));
+           5e-14 * diagonal(deviatoric_v).norm());
     return 0.5 * deviatoric_v.transpose() * deviatoric_v;
 }
 
@@ -36,16 +36,23 @@ double Invariants<KelvinVectorSize>::J3(
     Eigen::Matrix<double, KelvinVectorSize, 1> const& deviatoric_v)
 {
     assert(std::abs(trace(deviatoric_v)) <=
-           1e-16 * std::abs(deviatoric_v.squaredNorm()));
+           5e-14 * diagonal(deviatoric_v).norm());
     return determinant(deviatoric_v);
 }
 
+template <int KelvinVectorSize>
+Eigen::Vector3d Invariants<KelvinVectorSize>::diagonal(
+    Eigen::Matrix<double, KelvinVectorSize, 1> const& v)
+{
+    return v.template topLeftCorner<3, 1>();
+}
+
 /// Trace of the corresponding tensor.
 template <int KelvinVectorSize>
 double Invariants<KelvinVectorSize>::trace(
     Eigen::Matrix<double, KelvinVectorSize, 1> const& v)
 {
-    return v.template topLeftCorner<3, 1>().sum();
+    return diagonal(v).sum();
 }
 
 //
diff --git a/MathLib/KelvinVector.h b/MathLib/KelvinVector.h
index bf77b34534b..0b2034610a2 100644
--- a/MathLib/KelvinVector.h
+++ b/MathLib/KelvinVector.h
@@ -99,6 +99,11 @@ struct Invariants final
 
     /// Trace of the corresponding tensor.
     static double trace(Eigen::Matrix<double, KelvinVectorSize, 1> const& v);
+
+    /// Diagonal of the corresponding tensor which is always of length 3 in 2D
+    /// and 3D cases.
+    static Eigen::Vector3d diagonal(
+        Eigen::Matrix<double, KelvinVectorSize, 1> const& v);
 };
 
 //
-- 
GitLab