From 65af28fd7babca003b3796d83bed3ff7d8499d58 Mon Sep 17 00:00:00 2001
From: Emilie Sauvage <emilie.sauvage@uclouvain.be>
Date: Tue, 8 Jan 2013 17:12:28 +0000
Subject: [PATCH] Changed intersection so that intersection(M,M) = M
---
Geo/STensor3.cpp | 42 ++++++++++++++++++++++++++++++++++++++++++
1 file changed, 42 insertions(+)
diff --git a/Geo/STensor3.cpp b/Geo/STensor3.cpp
index 47a6da59ca..7096523b71 100644
--- a/Geo/STensor3.cpp
+++ b/Geo/STensor3.cpp
@@ -22,6 +22,8 @@ void STensor3::print (const char *s) const
(*this)(2,0), (*this)(2,1), (*this)(2,2));
}
+//// ORIGINAL CODE FROM GMSH
+/*
SMetric3 intersection (const SMetric3 &m1, const SMetric3 &m2)
{
SMetric3 im1 = m1.invert();
@@ -39,6 +41,46 @@ SMetric3 intersection (const SMetric3 &m1, const SMetric3 &m2)
return iv;
}
+*/
+
+//// MODIFIED INTERSECTION
+SMetric3 intersection (const SMetric3 &m1, const SMetric3 &m2)
+{
+ SMetric3 im1 = m1.invert();
+ fullMatrix<double> V(3,3);
+ fullVector<double> S(3);
+ im1 *= m2;
+ im1.eig(V,S,true);
+ SVector3 v0(V(0,0),V(1,0),V(2,0));
+ SVector3 v1(V(0,1),V(1,1),V(2,1));
+ SVector3 v2(V(0,2),V(1,2),V(2,2));
+ double l0 = std::max(dot(v0,m1,v0),dot(v0,m2,v0));
+ double l1 = std::max(dot(v1,m1,v1),dot(v1,m2,v1));
+ double l2 = std::max(dot(v2,m1,v2),dot(v2,m2,v2));
+
+
+ // Correction as explained in the PhD thesis of Frederic Alauzet p.16
+ double max_eig = std::max(l0, std::max(l1,l2));
+ double min_eig = std::min(l0, std::min(l1,l2));
+ double range_eig = (max_eig - min_eig);
+
+ if( range_eig < 1.0e-1) // Change this value if you think it's too big ...
+ {
+ if(max_eig <= 1.0)
+ {
+ return m1;
+ }
+
+ else
+ {
+ return m2;
+ }
+ }
+
+ SMetric3 iv(l0,l1,l2,v0,v1,v2);
+
+ return iv;
+}
// preserve orientation of m1 !!!
SMetric3 intersection_conserveM1 (const SMetric3 &m1, const SMetric3 &m2)
--
GitLab