diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 7c1d108e5a9e883160545c47f08678ca0df235e5..6ee192e4bad1d6ca1219ece0742013076568a234 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -1,4 +1,4 @@
-// $Id: GFace.cpp,v 1.24 2006-12-03 00:04:31 geuzaine Exp $
+// $Id: GFace.cpp,v 1.25 2006-12-03 00:35:56 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -399,37 +399,34 @@ double GFace::curvature (const SPoint2 &param) const
 }
 
 void GFace::XYZtoUV(const double X, const double Y, const double Z, 
-		    double &U, double &V,
-		    const double relax) const
+		    double &U, double &V, const double relax, 
+		    bool onSurface) const
 {
   const double Precision = 1.e-8;
   const int MaxIter = 25;
   const int NumInitGuess = 11;
 
-  double Unew = 0., Vnew = 0., err, err_xyz;
-  int iter;
+  double Unew = 0., Vnew = 0.;
   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,0,1};
-  
+  double init[NumInitGuess] = {0.487, 0.6, 0.4, 0.7, 0.3, 0.8, 0.2, 0.9, 0.1, 0, 1};
   Range<double> ru = parBounds(0);
   Range<double> rv = parBounds(1);
-  umin = ru.low();
-  umax = ru.high();
-  vmin = rv.low();
-  vmax = rv.high();
+  double umin = ru.low();
+  double umax = ru.high();
+  double vmin = rv.low();
+  double vmax = rv.high();
 
   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;
+      int err = 1.0, err_xyz = 0.;
+      int iter = 1;
       
       while(err > Precision && iter < MaxIter) {
-	GPoint P = point(U,V);
-	Pair<SVector3,SVector3> der = firstDer(SPoint2(U,V)); 
+	GPoint P = point(U, V);
+	Pair<SVector3, SVector3> der = firstDer(SPoint2(U, V)); 
 	mat[0][0] = der.left().x();
 	mat[0][1] = der.left().y();
 	mat[0][2] = der.left().z();
@@ -449,7 +446,8 @@ void GFace::XYZtoUV(const double X, const double Y, const double Z,
 	   jac[2][1] * (Z - P.z()));
 	
 	err = DSQR(Unew - U) + DSQR(Vnew - V);
-	err_xyz = DSQR(X - P.x()) + DSQR(Y - P.y()) + DSQR(Z - P.z());
+	if(onSurface)
+	  err_xyz = DSQR(X - P.x()) + DSQR(Y - P.y()) + DSQR(Z - P.z());
 	iter++;
 	U = Unew;
 	V = Vnew;
@@ -458,7 +456,7 @@ void GFace::XYZtoUV(const double X, const double Y, const double Z,
       if(iter < MaxIter && err <= Precision &&
 	 Unew <= umax && Vnew <= vmax && 
 	 Unew >= umin && Vnew >= vmin){
-	if(err_xyz > 1.e-5)
+	if(onSurface && err_xyz > Precision * CTX.lc)
 	  Msg(WARNING,"converged for i=%d j=%d (err=%g iter=%d), but err_xyz = %g", 
 	      i, j, err, iter, err_xyz);
 	return;	
@@ -466,21 +464,20 @@ void GFace::XYZtoUV(const double X, const double Y, const double Z,
     }
   }
   
-  if(relax < 1.e-6)
+  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(X, Y, Z, U, V, 0.75 * relax);
+    XYZtoUV(X, Y, Z, U, V, 0.75 * relax, onSurface);
   }  
 }
 
-
 SPoint2 GFace::parFromPoint(const SPoint3 &p) const
 {
-  double U,V;
+  double U, V;
   
-  XYZtoUV(p.x(),p.y(),p.z(),U,V,1.0);
+  XYZtoUV(p.x(), p.y(), p.z(), U, V, 1.0);
 
-  return SPoint2(U,V);
+  return SPoint2(U, V);
 }
-
diff --git a/Geo/GFace.h b/Geo/GFace.h
index 94803e19ae3106b3acefd99b19fa9a764bcf6a88..886e55e50d5a3b206b5ace865173bd9189062799 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -76,10 +76,10 @@ class GFace : public GEntity
   virtual int dim() const {return 2;}
   virtual void setVisibility(char val, bool recursive=false);
 
-  // compute XYZ from parametric UV
+  // compute the parameters UV from a point XYZ
   void XYZtoUV(const double X, const double Y, const double Z, 
-	       double &U, double &V,
-	       const double relax) const;
+	       double &U, double &V, const double relax,
+	       bool onSurface=true) const;
 
   // The bounding box
   virtual SBoundingBox3d bounds() const; 
@@ -94,7 +94,7 @@ class GFace : public GEntity
 
   // Return the parmater location on the face given a point in space
   // that is on the face.
-  virtual SPoint2            parFromPoint(const SPoint3 &) const;
+  virtual SPoint2 parFromPoint(const SPoint3 &) const;
 
   // True if the parameter value is interior to the face.
   virtual int containsParam(const SPoint2 &pt) const = 0;
diff --git a/Mesh/meshGFaceTransfinite.cpp b/Mesh/meshGFaceTransfinite.cpp
index 335b9e2a3a7cff6e3c1d1a777d5729e28812044e..833d9fbb5979e52b6392997cd7f5bbeb11947ed4 100644
--- a/Mesh/meshGFaceTransfinite.cpp
+++ b/Mesh/meshGFaceTransfinite.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFaceTransfinite.cpp,v 1.10 2006-12-02 22:48:25 geuzaine Exp $
+// $Id: meshGFaceTransfinite.cpp,v 1.11 2006-12-03 00:35:56 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -221,9 +221,8 @@ int MeshTransfiniteSurface( GFace *gf)
 	double zp = TRAN_TRI(m_vertices[iP1]->z(), m_vertices[iP2]->z(), 
 			     m_vertices[iP3]->z(), m_vertices[N1]->z(),
 			     m_vertices[N2]->z(), m_vertices[N3]->z(), u, v);
-	SPoint2 param = gf->parFromPoint(SPoint3(xp, yp, zp));
-	double Up = param.x();
-	double Vp = param.y();
+	double Up, Vp;
+	gf->XYZtoUV(xp, yp, zp, Up, Vp, 1.0, false);
 #endif
 	GPoint gp = gf->point(SPoint2(Up, Vp));
 	MFaceVertex *newv = new MFaceVertex(gp.x(), gp.y(), gp.z(), gf, Up, Vp);