diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index 3dd3273077a4277afb184031c15576da8b96d506..59b7afc3aa53edd4aaab7796855a7f8244db64e4 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 2b33eef5629c330d1c48b6e336af17e77ebf3129..b1298f63675cf1adaa652287d43df54cccb4b843 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 490bdccdbb33d6920b71c4f5c860cba68cc5a2ca..94d58a6b62eb1744eb935331a42a7c66c9133ae5 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