L'introduction d'une relaxation dans XYtoUV permet de resoudre pas mal de...

L'introduction d'une relaxation dans XYtoUV permet de resoudre pas mal de problemes (cf. bench/bugs/xy2uv

Mesh/2D_Mesh.cpp

-// $Id: 2D_Mesh.cpp,v 1.39 2002-02-01 14:34:05 remacle Exp $
+// $Id: 2D_Mesh.cpp,v 1.40 2002-02-12 20:11:34 geuzaine Exp $
 
 /*
    Maillage Delaunay d'une surface (Point insertion Technique)
@@ -145,7 +145,7 @@ void Calcule_Z (void *data, void *dum){
   if (v->Frozen || v->Num < 0)
     return;
 
-  XYtoUV (THESUPPORT, &v->Pos.X, &v->Pos.Y, &U, &V, &Z);
+  XYtoUV (THESUPPORT, &v->Pos.X, &v->Pos.Y, &U, &V, &Z, 1.0);
   v->Pos.Z = Z;
 }
 

Mesh/SecondOrder.cpp

-// $Id: SecondOrder.cpp,v 1.7 2001-10-29 08:52:20 geuzaine Exp $
+// $Id: SecondOrder.cpp,v 1.8 2002-02-12 20:11:34 geuzaine Exp $
 
 #include "Gmsh.h"
 #include "Geo.h"
@@ -45,8 +45,8 @@ Vertex *middleface (Vertex * v1, Vertex * v2){
   if (THES->Typ == MSH_SURF_PLAN)
     return NULL;
 
-  XYZtoUV ( THES , v1->Pos.X , v1->Pos.Y , v1->Pos.Z, &U1 , &V1 );    
-  XYZtoUV ( THES , v2->Pos.X , v2->Pos.Y , v2->Pos.Z, &U2 , &V2 );    
+  XYZtoUV ( THES , v1->Pos.X , v1->Pos.Y , v1->Pos.Z, &U1 , &V1 , 1.0);    
+  XYZtoUV ( THES , v2->Pos.X , v2->Pos.Y , v2->Pos.Z, &U2 , &V2 , 1.0);    
 
   U = 0.5 *(U1+U2);
   V = 0.5 *(V1+V2);

Mesh/Utils.cpp

-// $Id: Utils.cpp,v 1.8 2001-12-16 05:16:37 remacle Exp $
+// $Id: Utils.cpp,v 1.9 2002-02-12 20:11:34 geuzaine Exp $
 
 #include "Gmsh.h"
 #include "Numeric.h"
@@ -178,7 +178,7 @@ end:
 
 
 #define  Precision 1.e-10
-#define  MaxIter 20
+#define  MaxIter 50
 
 void find_bestuv (Surface * s, double X, double Y,
                   double *U, double *V, double *Z, int N){
@@ -253,11 +253,24 @@ void invert_singular_matrix(double **M, int n, double **I){
   free_dvector(W,1,n);
 }
 
-void XYZtoUV (Surface *s, double X, double Y, double Z, double *U, double *V) {
+void XYZtoUV (Surface *s, double X, double Y, double Z, double *U, double *V,
+	      double relax) {
   double Unew,Vnew,err;
   int    iter;
   Vertex D_u,D_v,P;
   double **mat, **jac ;
+  double umin, umax, vmin, vmax;
+
+  if (s->Typ == MSH_SURF_NURBS){
+    umin = s->ku[0];
+    umax = s->ku[s->OrderU + s->Nu];
+    vmin = s->kv[0];
+    vmax = s->kv[s->OrderV + s->Nv];
+  }
+  else{
+    umin = vmin = 0.0;
+    umax = vmax = 1.0;
+  }
 
   mat = dmatrix(1,3,1,3);
   jac = dmatrix(1,3,1,3);
@@ -282,8 +295,10 @@ void XYZtoUV (Surface *s, double X, double Y, double Z, double *U, double *V) {
     mat[3][3] = 0.; 
     invert_singular_matrix(mat,3,jac);
 
-    Unew = *U + jac[1][1] * (X-P.Pos.X) + jac[2][1] * (Y-P.Pos.Y) + jac[3][1] * (Z-P.Pos.Z) ;
-    Vnew = *V + jac[1][2] * (X-P.Pos.X) + jac[2][2] * (Y-P.Pos.Y) + jac[3][2] * (Z-P.Pos.Z) ;
+    Unew = *U + relax * (jac[1][1] * (X-P.Pos.X) + jac[2][1] * (Y-P.Pos.Y) + 
+			 jac[3][1] * (Z-P.Pos.Z)) ;
+    Vnew = *V + relax * (jac[1][2] * (X-P.Pos.X) + jac[2][2] * (Y-P.Pos.Y) +
+			 jac[3][2] * (Z-P.Pos.Z)) ;
 
     err = DSQR(Unew - *U) + DSQR(Vnew - *V) ;
 
@@ -292,18 +307,23 @@ void XYZtoUV (Surface *s, double X, double Y, double Z, double *U, double *V) {
     *V = Vnew;
   }
 
-  if(iter > 10){
-    if(iter == MaxIter) Msg(WARNING, "Could not converge in XYZtoUV");
-    else Msg(WARNING, "Many (%d) iterations in XYZtoUV", iter);
-  }
-
   free_dmatrix(mat,1,3,1,3);
   free_dmatrix(jac,1,3,1,3);
 
+  if(iter == MaxIter ||
+     (Unew > umax || Vnew > vmax || Unew < umin || Vnew < vmin)){
+    if(relax < 1.e-12)
+      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);
+    }
+  }
+
 }
 
 void XYtoUV (Surface * s, double *X, double *Y,
-             double *U, double *V, double *Z){
+             double *U, double *V, double *Z, double relax){
 
   double det, Unew, Vnew, err, mat[2][2], jac[2][2];
   int iter;
@@ -345,8 +365,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 + 1.0 * (jac[0][0] * (*X - P.Pos.X) + jac[1][0] * (*Y - P.Pos.Y));
-    Vnew = *V + 1.0 * (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);
     
@@ -357,33 +377,30 @@ void XYtoUV (Surface * s, double *X, double *Y,
   
   *Z = P.Pos.Z;
 
-  if(iter > 10){
-    if(iter == MaxIter) Msg(WARNING, "Could not converge in XYtoUV");
-    else Msg(WARNING, "Many (%d) iterations in XYtoUV...", iter);
-  }
-
   int thresh = Unew > umax || Vnew > vmax || Unew < umin || Vnew < vmin;
 
-  if (thresh){
-    //Msg(WARNING, "(U,V) thresholded in XYtoUV (surface mesh may be wrong)");
-    if(Unew > umax) *U = umax;
-    if(Vnew > vmax) *V = vmax;
-    if(Unew < umin) *U = umin;
-    if(Vnew < vmin) *V = vmin;
+  if(iter == MaxIter || thresh){
+    if(relax < 1.e-6){
+      Msg(GERROR, "Could not converge: surface mesh will probably be wrong");
+      if(Unew > umax) *U = umax;
+      if(Vnew > vmax) *V = vmax;
+      if(Unew < umin) *U = umin;
+      if(Vnew < vmin) *V = vmin;
+      find_bestuv (s, *X, *Y, U, V, Z, 30);
+      P = InterpolateSurface (s, *U, *V, 0, 0);
+      *X = P.Pos.X;
+      *Y = P.Pos.Y;
+      *Z = P.Pos.Z;
+    }
+    else{
+      Msg(INFO, "Relaxation factor = %g", 0.75*relax);
+      XYtoUV (s, X, Y, U, V, Z, 0.75*relax);
+    }
   }
 
-#if 1
-  if (iter == MaxIter || thresh){
-    Msg(WARNING, "Entering rescue mode in XYtoUV: surface mesh may be wrong");
-    find_bestuv (s, *X, *Y, U, V, Z, 30);
-    P = InterpolateSurface (s, *U, *V, 0, 0);
-    *X = P.Pos.X;
-    *Y = P.Pos.Y;
-    *Z = P.Pos.Z;
-  }
-#endif
 }
 
+
 int Oriente (List_T * cu, double n[3]){
   int N, i, a, b, c;
   double cosa, sina, sum, v[3], w[3], u[3];

Mesh/Utils.h

@@ -7,8 +7,9 @@ void MeanPlane(List_T *point, Surface *s);
 void find_bestuv (Surface * s, double X, double Y,
                   double *U, double *V, double *Z, int N);
 void XYtoUV (Surface * s, double *X, double *Y,
-             double *U, double *V, double *Z);
-void XYZtoUV (Surface *s, double X, double Y, double Z, double *U, double *V);
+             double *U, double *V, double *Z, double relax);
+void XYZtoUV (Surface *s, double X, double Y, double Z, 
+	      double *U, double *V, double relax);
 int Oriente (List_T * cu, double n[3]);
 double angle_3p (Vertex * V, Vertex * P1, Vertex * P2);
 double angle_plan (Vertex * V, Vertex * P1, Vertex * P2, double n[3]);