From a8dd8a8fa259468e76b61d45dce7d064e4274b6d Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Mon, 4 Mar 2013 17:11:23 +0000
Subject: [PATCH] removed a bug in cross fields computation (surfaces) hope it
 is the last one ;-)

---
 Mesh/BackgroundMesh.cpp               | 18 +++++++++++++++++-
 Mesh/meshGRegionDelaunayInsertion.cpp |  4 ++--
 Mesh/surfaceFiller.cpp                |  6 ++++--
 3 files changed, 23 insertions(+), 5 deletions(-)

diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index 3dd3273077..59b7afc3aa 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -763,6 +763,13 @@ void backgroundMesh::propagateCrossFieldByDistance(GFace *_gf)
 #endif
 }
 
+inline double myAngle (const SVector3 &a, const SVector3 &b, const SVector3 &d){
+  double cosTheta = dot(a,b);
+  double sinTheta = dot(crossprod(a,b),d);
+  return atan2 (sinTheta,cosTheta);  
+}
+
+
 void backgroundMesh::propagatecrossField(GFace *_gf)
 {
 
@@ -783,15 +790,24 @@ void backgroundMesh::propagatecrossField(GFace *_gf)
         reparamMeshEdgeOnFace(v[0],v[1],_gf,p1,p2);
 	Pair<SVector3, SVector3> der = _gf->firstDer((p1+p2)*.5);
 	SVector3 t1 = der.first();
+	SVector3 s2 = der.second();
+	SVector3 n = crossprod(t1,s2);
+	n.normalize();
 	SVector3 t2 (v[1]->x()-v[0]->x(),v[1]->y()-v[0]->y(),v[1]->z()-v[0]->z());
 	t1.normalize();
 	t2.normalize();
-	double _angle = angle (t1,t2);	
+	double _angle = myAngle (t1,t2,n);	
+	//	printf("GFACE %d %g %g %g %g\n",_gf->tag(),t1.x(),t1.y(),t1.z(),_angle*180/M_PI);
 	//	printf("angle = %12.5E\n",_angle);
 	//	angle_ = atan2 ( p1.y()-p2.y() , p1.x()-p2.x() );
         //double angle = atan2 ( v[0]->y()-v[1]->y() , v[0]->x()- v[1]->x() );
         //double angle = atan2 ( v0->y()-v1->y() , v0->x()- v1->x() );
         crossField2d::normalizeAngle (_angle);
+	//	SVector3 s2 = der.second();
+	//	s2.normalize();
+	//	SVector3 x = t1 * cos (_angle) + s2 * sin (_angle);
+	//	printf("%g %g %g vs %g %g %g\n",x.x(),x.y(),x.z(),t2.x(),t2.y(),t2.z());
+	//	printf("GFACE %d GEDGE %d %g %g %g %g\n",_gf->tag(),(*it)->tag(),t1.x(),t1.y(),t1.z(),_angle*180/M_PI);
         for (int i=0;i<2;i++){
           std::map<MVertex*,double>::iterator itc = _cosines4.find(v[i]);
           std::map<MVertex*,double>::iterator its = _sines4.find(v[i]);
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index 2b33eef562..b1298f6367 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -305,8 +305,8 @@ bool insertVertexB(std::list<faceXtet> &shell,
     ++it;
   }
   // OK, the cavity is star shaped
-  if (onePointIsTooClose)printf("One point is too close\n");
-  if (fabs(oldVolume - newVolume) > 1.e-10 * oldVolume)printf("Volume do not match %22.15E %22.15E %22.15E\n",oldVolume,newVolume,fabs(oldVolume-newVolume)/newVolume);
+  //  if (onePointIsTooClose)printf("One point is too close\n");
+  //  if (fabs(oldVolume - newVolume) > 1.e-10 * oldVolume)printf("Volume do not match %22.15E %22.15E %22.15E\n",oldVolume,newVolume,fabs(oldVolume-newVolume)/newVolume);
   if (fabs(oldVolume - newVolume) < 1.e-10 * oldVolume &&
       !onePointIsTooClose){
     connectTets_vector(new_cavity.begin(), new_cavity.end());
diff --git a/Mesh/surfaceFiller.cpp b/Mesh/surfaceFiller.cpp
index 490bdccdbb..94d58a6b62 100644
--- a/Mesh/surfaceFiller.cpp
+++ b/Mesh/surfaceFiller.cpp
@@ -210,12 +210,15 @@ bool compute4neighbors (GFace *gf,   // the surface
   SVector3 n = crossprod(s1,s2);
   n.normalize();
   
+  //  printf("%d %g %g %g\n",gf->tag(),s1.x(),s1.y(),s1.z());
+
   for (int DIR = 0 ; DIR < NUMDIR ; DIR ++){  
     double quadAngle  = backgroundMesh::current()->getAngle (midpoint[0],midpoint[1],0) + DIRS[DIR];
     
     // normalize vector t1 that is tangent to gf at midpoint
     SVector3 t1 = s1 * cos (quadAngle) + s2 * sin (quadAngle);
     t1.normalize();
+    //    printf("%d %g %g %g -- %g %g %g\n",gf->tag(),s1.x(),s1.y(),s1.z(),t1.x(),t1.y(),t1.z());
     
     // compute the second direction t2 and normalize (t1,t2,n) is the tangent frame
     SVector3 t2 = crossprod(t1,n);
@@ -230,8 +233,7 @@ bool compute4neighbors (GFace *gf,   // the surface
     double N = dot(s2,s2);
     double E = dot(s1,s2);
     double metric[2][2] = {{M,E},{E,N}};
-    
-    
+        
     // compute covariant coordinates of t1 and t2
     // t1 = a s1 + b s2 -->
     // t1 . s1 = a M + b E
-- 
GitLab