diff --git a/Common/Context.h b/Common/Context.h
index 2fad261388eba2088268596d32471a3e8f449609..a468e0c57749d82b4ac99cb8cc187147316ba38b 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -152,6 +152,7 @@ class Context_T {
     double tolerance;
     double snap[3];
     int occ_fix_small_edges, occ_fix_small_faces, occ_sew_faces;
+    int sphere_geodesic;
   } geom;
 
   // mesh options 
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index b20666a9cf95b7cc57bde2ec374783acf1806ff1..a852255dc9c99968f448085a5c96a113e266ad5b 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -871,6 +871,8 @@ StringXNumber GeometryOptions_Number[] = {
     "Snapping grid spacing along the Y-axis" },
   { F|O, "SnapZ" , opt_geometry_snap2 , 0.1 ,
     "Snapping grid spacing along the Z-axis" },
+  { F|O, "SphereGeodesic" , opt_geometry_sphere_geodesic , 0. ,
+    "Use projection to compute sphere geodesic" },
   { F|O, "Surfaces" , opt_geometry_surfaces , 0. , 
     "Display geometry surfaces?" },
   { F|O, "SurfaceNumbers" , opt_geometry_surfaces_num , 0. , 
diff --git a/Common/Options.cpp b/Common/Options.cpp
index 486f8479ba647da12312e712230df0e2266c6a49..4c0b3782dc6b6f79d1f118c45170df9123c611de 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -1,4 +1,4 @@
-// $Id: Options.cpp,v 1.399 2008-06-05 12:37:03 samtech Exp $
+// $Id: Options.cpp,v 1.400 2008-06-20 05:51:36 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -4035,6 +4035,13 @@ double opt_geometry_tolerance(OPT_ARGS_NUM)
   return CTX.geom.tolerance;
 }
 
+double opt_geometry_sphere_geodesic(OPT_ARGS_NUM)
+{
+  if(action & GMSH_SET)
+    CTX.geom.sphere_geodesic = val;
+  return CTX.geom.sphere_geodesic;
+}
+
 double opt_geometry_normals(OPT_ARGS_NUM)
 {
   if(action & GMSH_SET)
diff --git a/Common/Options.h b/Common/Options.h
index 5824ef267a8729f570aa491180d4858e1082a406..0beacbcfca7b44f28c569f1a69a683cd0d900c05 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -399,6 +399,7 @@ double opt_general_light53(OPT_ARGS_NUM);
 double opt_geometry_auto_coherence(OPT_ARGS_NUM);
 double opt_geometry_highlight_orphans(OPT_ARGS_NUM);
 double opt_geometry_tolerance(OPT_ARGS_NUM);
+double opt_geometry_sphere_geodesic(OPT_ARGS_NUM);
 double opt_geometry_normals(OPT_ARGS_NUM);
 double opt_geometry_tangents(OPT_ARGS_NUM);
 double opt_geometry_points(OPT_ARGS_NUM);
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index b3af0fe7cae076b7b494f309d9e39b6d479fe22a..ff8ac7b3e313fa97532e2cdcae7f3d1536d525d7 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -1,4 +1,4 @@
-// $Id: Geo.cpp,v 1.114 2008-06-19 15:58:41 geuzaine Exp $
+// $Id: Geo.cpp,v 1.115 2008-06-20 05:51:36 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -2971,10 +2971,12 @@ bool ProjectPointOnSurface(Surface *s, Vertex &p, double u[2])
   double UMAX = 1.;
   double VMIN = 0.;
   double VMAX = 1.;
-  while(1) {
+  double eps = 1.e-6;
+  int nn = 0;
+  while(++nn < 20) {
     newt(x, 2, &check, projectPS);
     vv = InterpolateSurface(s, x[1], x[2], 0, 0);
-    if(x[1] >= UMIN && x[1] <= UMAX && x[2] >= VMIN && x[2] <= VMAX)
+    if(x[1] > UMIN-eps && x[1] < UMAX+eps && x[2] > VMIN-eps && x[2] < VMAX+eps)
       break;
     x[1] = UMIN + (UMAX - UMIN) * ((rand() % 10000) / 10000.);
     x[2] = VMIN + (VMAX - VMIN) * ((rand() % 10000) / 10000.);
@@ -2992,11 +2994,9 @@ bool ProjectPointOnSurface(Surface *s, Vertex &p, double u[2])
 
 void Projette(Vertex *v, double mat[3][3])
 {
-  double X, Y, Z;
-
-  X = v->Pos.X * mat[0][0] + v->Pos.Y * mat[0][1] + v->Pos.Z * mat[0][2];
-  Y = v->Pos.X * mat[1][0] + v->Pos.Y * mat[1][1] + v->Pos.Z * mat[1][2];
-  Z = v->Pos.X * mat[2][0] + v->Pos.Y * mat[2][1] + v->Pos.Z * mat[2][2];
+  double X = v->Pos.X * mat[0][0] + v->Pos.Y * mat[0][1] + v->Pos.Z * mat[0][2];
+  double Y = v->Pos.X * mat[1][0] + v->Pos.Y * mat[1][1] + v->Pos.Z * mat[1][2];
+  double Z = v->Pos.X * mat[2][0] + v->Pos.Y * mat[2][1] + v->Pos.Z * mat[2][2];
   v->Pos.X = X;
   v->Pos.Y = Y;
   v->Pos.Z = Z;
diff --git a/Geo/Makefile b/Geo/Makefile
index 0572e42d00758a29c6cf13fea9b3be665157d283..cba91cb51a420682d984cfd47596ffaae2e8168c 100644
--- a/Geo/Makefile
+++ b/Geo/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.213 2008-06-07 17:20:46 geuzaine Exp $
+# $Id: Makefile,v 1.214 2008-06-20 05:51:36 geuzaine Exp $
 #
 # Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 #
@@ -126,7 +126,7 @@ gmshFace.o: gmshFace.cpp GModel.h GVertex.h GEntity.h Range.h SPoint3.h \
   gmshSurface.h ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
   ../Common/ListUtils.h ../Common/TreeUtils.h ../Common/avl.h \
   ../Common/ListUtils.h ExtrudeParams.h ../Common/SmoothData.h \
-  GeoInterpolation.h ../Common/Message.h
+  GeoInterpolation.h ../Common/Message.h ../Common/Context.h
 gmshRegion.o: gmshRegion.cpp GModel.h GVertex.h GEntity.h Range.h \
   SPoint3.h SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h SVector3.h \
   GFace.h GEdgeLoop.h Pair.h GRegion.h gmshRegion.h Geo.h \
diff --git a/Geo/gmshFace.cpp b/Geo/gmshFace.cpp
index 1f5c9b1430d9e65e83ada2913c3d0c31b735199a..85c9c9ada36fafe69c297ce39488b7ed1d72e0dc 100644
--- a/Geo/gmshFace.cpp
+++ b/Geo/gmshFace.cpp
@@ -1,4 +1,4 @@
-// $Id: gmshFace.cpp,v 1.60 2008-06-10 12:59:12 remacle Exp $
+// $Id: gmshFace.cpp,v 1.61 2008-06-20 05:51:37 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -26,6 +26,9 @@
 #include "GeoInterpolation.h"
 #include "Numeric.h"
 #include "Message.h"
+#include "Context.h"
+
+extern Context_T CTX;
 
 gmshFace::gmshFace(GModel *m, Surface *face)
   : GFace(m, face->Num), s(face)
@@ -85,7 +88,7 @@ gmshFace::gmshFace(GModel *m, Surface *face)
         Msg::Error("Unknown point %d", v->Num);
     }
   }
-  isSphere = iSRuledSurfaceASphere (s,center,radius);
+  isSphere = iSRuledSurfaceASphere(s, center, radius);
 }
 
 double gmshFace::getMetricEigenvalue(const SPoint2 &pt)
@@ -246,21 +249,23 @@ GEntity::GeomType gmshFace::geomType() const
 // by default we assume that straight lines are geodesics
 SPoint2 gmshFace::geodesic(const SPoint2 &pt1 , const SPoint2 &pt2 , double t)
 {
-  if (isSphere){
-    GPoint gp1 = point (pt1.x(),pt1.y());
-    GPoint gp2 = point (pt2.x(),pt2.y());    
+  if(isSphere && CTX.geom.sphere_geodesic){
+    // FIXME: this is buggy -- remove the CTX option once we do it in
+    // a robust manner
+    GPoint gp1 = point(pt1.x(), pt1.y());
+    GPoint gp2 = point(pt2.x(), pt2.y());
     SPoint2 guess = pt1 + (pt2 - pt1) * t;
-    GPoint gp = closestPoint(SPoint3(gp1.x()+t*(gp2.x()-gp1.x()),
-				     gp1.y()+t*(gp2.y()-gp1.y()),
-				     gp1.z()+t*(gp2.z()-gp1.z())),(double*)guess);
-    return SPoint2(gp.u(),gp.v());
+    GPoint gp = closestPoint(SPoint3(gp1.x() + t * (gp2.x() - gp1.x()),
+				     gp1.y() + t * (gp2.y() - gp1.y()),
+				     gp1.z() + t * (gp2.z() - gp1.z())),
+			     (double*)guess);
+    return SPoint2(gp.u(), gp.v());
   }
   else{
     return pt1 + (pt2 - pt1) * t;
   }
 }
 
-
 int gmshFace::containsPoint(const SPoint3 &pt) const
 { 
   if(geomType() == Plane){
diff --git a/Mesh/Makefile b/Mesh/Makefile
index 2277626b21c5c09854c1bf4ea83edbfba843dcd4..096bb0bf07cd291b29e0fa9946dd6b1dfbbbedbe 100644
--- a/Mesh/Makefile
+++ b/Mesh/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.225 2008-06-10 08:37:34 remacle Exp $
+# $Id: Makefile,v 1.226 2008-06-20 05:51:37 geuzaine Exp $
 #
 # Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 #
@@ -147,22 +147,23 @@ meshGFace.o: meshGFace.cpp meshGFace.h meshGFaceBDS.h \
   meshGFaceDelaunayInsertion.h ../Geo/MElement.h ../Common/GmshDefines.h \
   ../Geo/MVertex.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
   ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/MFace.h ../Geo/MVertex.h \
-  ../Geo/SVector3.h meshGFaceOptimize.h DivideAndConquer.h \
-  BackgroundMesh.h ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h \
-  ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h \
-  ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/GEdge.h ../Geo/GEntity.h \
-  ../Geo/GVertex.h ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint2.h \
-  ../Geo/GFace.h ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/GEdgeLoop.h \
-  ../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
-  ../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEdge.h ../Geo/GFace.h \
-  ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SBoundingBox3d.h \
-  ../Common/Context.h ../Common/Message.h ../Numeric/Numeric.h \
-  ../Numeric/NumericEmbedded.h BDS.h ../Post/PView.h qualityMeasures.h \
-  Field.h ../Geo/Geo.h ../Geo/gmshSurface.h ../Geo/Pair.h ../Geo/Range.h \
-  ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/SVector3.h \
-  ../Geo/SBoundingBox3d.h ../Common/ListUtils.h ../Common/TreeUtils.h \
-  ../Common/avl.h ../Common/ListUtils.h ../Geo/SPoint2.h \
-  ../Geo/ExtrudeParams.h ../Common/SmoothData.h ../Common/OS.h
+  ../Geo/SVector3.h meshGFaceQuadrilateralize.h meshGFaceOptimize.h \
+  DivideAndConquer.h BackgroundMesh.h ../Geo/GVertex.h ../Geo/GEntity.h \
+  ../Geo/Range.h ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h \
+  ../Geo/SPoint3.h ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/GEdge.h \
+  ../Geo/GEntity.h ../Geo/GVertex.h ../Geo/SVector3.h ../Geo/SPoint3.h \
+  ../Geo/SPoint2.h ../Geo/GFace.h ../Geo/GEntity.h ../Geo/GPoint.h \
+  ../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h \
+  ../Geo/Pair.h ../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEdge.h \
+  ../Geo/GFace.h ../Geo/GRegion.h ../Geo/GEntity.h \
+  ../Geo/SBoundingBox3d.h ../Common/Context.h ../Common/Message.h \
+  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h BDS.h ../Post/PView.h \
+  qualityMeasures.h Field.h ../Geo/Geo.h ../Geo/gmshSurface.h \
+  ../Geo/Pair.h ../Geo/Range.h ../Geo/SPoint2.h ../Geo/SPoint3.h \
+  ../Geo/SVector3.h ../Geo/SBoundingBox3d.h ../Common/ListUtils.h \
+  ../Common/TreeUtils.h ../Common/avl.h ../Common/ListUtils.h \
+  ../Geo/SPoint2.h ../Geo/ExtrudeParams.h ../Common/SmoothData.h \
+  ../Common/OS.h
 meshGFaceTransfinite.o: meshGFaceTransfinite.cpp meshGFace.h \
   ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
   ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
@@ -231,6 +232,19 @@ meshGFaceOptimize.o: meshGFaceOptimize.cpp meshGFaceOptimize.h \
   ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/SVector3.h \
   ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/SPoint2.h ../Geo/SVector3.h \
   ../Geo/Pair.h BackgroundMesh.h
+meshGFaceQuadrilateralize.o: meshGFaceQuadrilateralize.cpp \
+  meshGFaceQuadrilateralize.h ../Common/Message.h ../Numeric/Numeric.h \
+  ../Numeric/NumericEmbedded.h meshGFaceDelaunayInsertion.h \
+  ../Geo/MElement.h ../Common/GmshDefines.h ../Geo/MVertex.h \
+  ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
+  ../Geo/SPoint3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
+  meshGFaceOptimize.h meshGFaceBDS.h BDS.h ../Geo/GFace.h \
+  ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
+  ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
+  ../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h \
+  ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/SVector3.h \
+  ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/SPoint2.h ../Geo/SVector3.h \
+  ../Geo/Pair.h ../Post/PView.h
 meshGRegion.o: meshGRegion.cpp meshGRegion.h \
   meshGRegionDelaunayInsertion.h ../Geo/MElement.h \
   ../Common/GmshDefines.h ../Geo/MVertex.h ../Geo/SPoint3.h \
diff --git a/doc/VERSIONS b/doc/VERSIONS
index a8cce0fd98af503ed99002448013ffba9db9019b..15a0297ce796e4f7c4e1eff115649e99df606916 100644
--- a/doc/VERSIONS
+++ b/doc/VERSIONS
@@ -1,6 +1,7 @@
-$Id: VERSIONS,v 1.410 2008-06-19 15:58:42 geuzaine Exp $
+$Id: VERSIONS,v 1.411 2008-06-20 05:51:37 geuzaine Exp $
 
-2.2.2 (): added geometrical transformations on volumes.
+2.2.2 (Jun 20, 2008): added geometrical transformations on volumes;
+fixed bug in high order mesh generation.
 
 2.2.1 (Jun 15, 2008): various small improvements (adaptive views, gui,
 code cleanup) and bug fixes (high order meshes, Netgen interface).
diff --git a/doc/texinfo/opt_geometry.texi b/doc/texinfo/opt_geometry.texi
index 7ca9f6ee738cc2ec6e9f52fb01e90ded8160947e..95733892965fc2bca3090a991727323830c18bb1 100644
--- a/doc/texinfo/opt_geometry.texi
+++ b/doc/texinfo/opt_geometry.texi
@@ -139,6 +139,11 @@ Snapping grid spacing along the Z-axis@*
 Default value: @code{0.1}@*
 Saved in: @code{General.OptionsFileName}
 
+@item Geometry.SphereGeodesic
+Use projection to compute sphere geodesic@*
+Default value: @code{0}@*
+Saved in: @code{General.OptionsFileName}
+
 @item Geometry.Surfaces
 Display geometry surfaces?@*
 Default value: @code{0}@*