Improved xytouv() and xyztouv() by trying different initial guesses in the
newton algo (before relaxing). Should be more robust than what we did before.

Mesh/Utils.cpp

-// $Id: Utils.cpp,v 1.24 2004-04-13 18:49:58 geuzaine Exp $
+// $Id: Utils.cpp,v 1.25 2004-04-18 03:23:02 geuzaine Exp $
 //
 // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 //
@@ -244,7 +244,8 @@ end:
 
 
 #define  Precision 1.e-10
-#define  MaxIter 50
+#define  MaxIter 25
+#define  NumInitGuess 9
 
 void find_bestuv(Surface * s, double X, double Y,
                  double *U, double *V, double *Z, int N)
@@ -361,6 +362,7 @@ void XYZtoUV(Surface * s, double X, double Y, double Z, double *U, double *V,
   Vertex D_u, D_v, P;
   double mat[3][3], jac[3][3];
   double umin, umax, vmin, vmax;
+  double init[NumInitGuess] = {0.487, 0.6, 0.4, 0.7, 0.3, 0.8, 0.2, 0.9, 0.1};
   
   if(s->Typ == MSH_SURF_NURBS) {
     umin = s->ku[0];
@@ -373,7 +375,11 @@ void XYZtoUV(Surface * s, double X, double Y, double Z, double *U, double *V,
     umax = vmax = 1.0;
   }
 
-  *U = *V = 0.487;
+  for(int i = 0; i < NumInitGuess; i++){
+    for(int j = 0; j < NumInitGuess; j++){
+    
+      *U = init[i];
+      *V = init[j];
       err = 1.0;
       iter = 1;
       
@@ -396,8 +402,8 @@ void XYZtoUV(Surface * s, double X, double Y, double Z, double *U, double *V,
 	Unew = *U + relax *
 	  (jac[0][0] * (X - P.Pos.X) + jac[1][0] * (Y - P.Pos.Y) +
 	   jac[2][0] * (Z - P.Pos.Z));
-    Vnew =
-      *V + relax * (jac[0][1] * (X - P.Pos.X) + jac[1][1] * (Y - P.Pos.Y) +
+	Vnew = *V + relax * 
+	  (jac[0][1] * (X - P.Pos.X) + jac[1][1] * (Y - P.Pos.Y) +
 	   jac[2][1] * (Z - P.Pos.Z));
 	
 	err = DSQR(Unew - *U) + DSQR(Vnew - *V);
@@ -407,15 +413,21 @@ void XYZtoUV(Surface * s, double X, double Y, double Z, double *U, double *V,
 	*V = Vnew;
       }
 
-  if(iter == MaxIter ||
-     (Unew > umax || Vnew > vmax || Unew < umin || Vnew < vmin)) {
-    if(relax < 1.e-12)
+      if(iter < MaxIter && err <= Precision &&
+	 Unew <= umax && Vnew <= vmax && 
+	 Unew >= umin && Vnew >= vmin){
+	//printf("converged for i=%d j=%d (err=%g iter=%d)\n", i, j, err, iter);
+	return;	
+      }
+    }
+  }
+
+  if(relax < 1.e-6)
     Msg(GERROR, "Could not converge: surface mesh will be wrong");
   else {
     Msg(INFO, "Relaxation factor = %g", 0.75 * relax);
     XYZtoUV(s, X, Y, Z, U, V, 0.75 * relax);
   }
-  }
 
 }
 
@@ -426,6 +438,7 @@ void XYtoUV(Surface * s, double *X, double *Y,
   int iter;
   Vertex D_u, D_v, P;
   double umin, umax, vmin, vmax;
+  double init[NumInitGuess] = {0.487, 0.6, 0.4, 0.7, 0.3, 0.8, 0.2, 0.9, 0.1};
 
   if(s->Typ == MSH_SURF_NURBS) {
     umin = s->ku[0];
@@ -438,7 +451,11 @@ void XYtoUV(Surface * s, double *X, double *Y,
     umax = vmax = 1.0;
   }
 
-  *U = *V = 0.487;
+  for(int i = 0; i < NumInitGuess; i++){
+    for(int j = 0; j < NumInitGuess; j++){
+
+      *U = init[i];
+      *V = init[j];
       err = 1.0;
       iter = 1;
 
@@ -462,10 +479,8 @@ void XYtoUV(Surface * s, double *X, double *Y,
 	jac[1][0] = -mat[1][0] / det;
 	jac[1][1] = mat[0][0] / det;
 	
-    Unew =
-      *U + relax * (jac[0][0] * (*X - P.Pos.X) + jac[1][0] * (*Y - P.Pos.Y));
-    Vnew =
-      *V + relax * (jac[0][1] * (*X - P.Pos.X) + jac[1][1] * (*Y - P.Pos.Y));
+	Unew = *U + relax * (jac[0][0] * (*X - P.Pos.X) + jac[1][0] * (*Y - P.Pos.Y));
+	Vnew = *V + relax * (jac[0][1] * (*X - P.Pos.X) + jac[1][1] * (*Y - P.Pos.Y));
 	
 	err = DSQR(Unew - *U) + DSQR(Vnew - *V);
 	
@@ -476,9 +491,15 @@ void XYtoUV(Surface * s, double *X, double *Y,
 
       *Z = P.Pos.Z;
 
-  int thresh = Unew > umax || Vnew > vmax || Unew < umin || Vnew < vmin;
+      if(iter < MaxIter && err <= Precision &&
+	 Unew <= umax && Vnew <= vmax && 
+	 Unew >= umin && Vnew >= vmin){
+	//printf("converged for i=%d j=%d (err=%g iter=%d)\n", i, j, err, iter);
+	return;	
+      }
+    }
+  }
 
-  if(iter == MaxIter || thresh) {
   if(relax < 1.e-6) {
     Msg(GERROR, "Could not converge: surface mesh will probably be wrong");
     if(Unew > umax)
@@ -499,11 +520,9 @@ void XYtoUV(Surface * s, double *X, double *Y,
     Msg(INFO, "Relaxation factor = %g", 0.75 * relax);
     XYtoUV(s, X, Y, U, V, Z, 0.75 * relax);
   }
-  }
 
 }
 
-
 int Oriente(List_T * cu, double n[3])
 {
   int N, i, a, b, c;