*** empty log message ***

Geo/GFace.h

@@ -102,7 +102,7 @@ class GFace : public GEntity
   virtual GPoint point(double par1, double par2) const = 0;
   virtual GPoint point(const SPoint2 &pt) const { return point(pt.x(), pt.y()); }
 
-  // computes, in parametric space, the interpolation from pt1 to pt2 alon the geodesic
+  // computes, in parametric space, the interpolation from pt1 to pt2 along a geodesic
   // of the surface
   virtual SPoint2 geodesic(const SPoint2 &pt1, const SPoint2 &pt2, double t);
 

Geo/GeoInterpolation.cpp

-// $Id: GeoInterpolation.cpp,v 1.36 2008-05-04 08:31:13 geuzaine Exp $
+// $Id: GeoInterpolation.cpp,v 1.37 2008-06-10 12:59:12 remacle Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -455,6 +455,55 @@ void TransfiniteSph(Vertex S, Vertex center, Vertex *T)
   T->Pos.Z = center.Pos.Z + r * dirz;
 }
 
+bool iSRuledSurfaceASphere(Surface *s, SPoint3 &center, double &radius){
+  if(s->Typ != MSH_SURF_REGL && s->Typ != MSH_SURF_TRIC)return false;
+
+  bool isSphere = true;
+  Vertex *O = 0, OO;
+  Curve *C[4] = {0, 0, 0, 0};
+  for(int i = 0; i < std::min(List_Nbr(s->Generatrices), 4); i++)
+    List_Read(s->Generatrices, i, &C[i]);
+
+  if(List_Nbr(s->RuledSurfaceOptions) == 3) {
+    // it's on a sphere: get the center
+    List_Read(s->RuledSurfaceOptions, 0, & ((double *)center)[0]);
+    List_Read(s->RuledSurfaceOptions, 1, & ((double *)center)[1]);
+    List_Read(s->RuledSurfaceOptions, 2, & ((double *)center)[2]);
+    O = &OO;
+  }
+  else{
+    // try to be intelligent (hum)
+    for(int i = 0; i < std::min(List_Nbr(s->Generatrices), 4); i++) {
+      if(C[i]->Typ != MSH_SEGM_CIRC && C[i]->Typ != MSH_SEGM_CIRC_INV){
+	isSphere = false;
+      }
+      else if(isSphere){
+	if(!i){
+	  List_Read(C[i]->Control_Points, 1, &O);
+	  ((double *)center)[0]= O->Pos.X;
+	  ((double *)center)[1]= O->Pos.Y;
+	  ((double *)center)[2]= O->Pos.Z;
+	}
+	else{
+	  Vertex *tmp;
+	  List_Read(C[i]->Control_Points, 1, &tmp);
+	  if(compareVertex(&O, &tmp))
+	    isSphere = false;
+	}
+      }
+    }
+  }
+  if (isSphere){
+    Vertex *p = C[0]->beg;
+    radius = sqrt ((p->Pos.X - center.x())+
+		   (p->Pos.Y - center.y())+
+		   (p->Pos.Z - center.z()));
+  }
+
+  return isSphere;
+}
+
+
 Vertex InterpolateRuledSurface(Surface *s, double u, double v)
 {
   Curve *C[4] = {0, 0, 0, 0};

Geo/GeoInterpolation.h

@@ -22,6 +22,7 @@
 
 #include "Geo.h"
 
+bool iSRuledSurfaceASphere(Surface *s, SPoint3 &center, double &radius);
 Vertex InterpolateCurve(Curve *Curve, double u, int derivee);
 Vertex InterpolateSurface(Surface *s, double u, double v, int derivee, int u_v);
 SPoint2 InterpolateCubicSpline(Vertex * v[4], double t, double mat[4][4],

Geo/MElement.cpp

-// $Id: MElement.cpp,v 1.71 2008-06-10 08:37:33 remacle Exp $
+// $Id: MElement.cpp,v 1.72 2008-06-10 12:59:12 remacle Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -891,12 +891,8 @@ void MTriangleN::getFaceRep(int num, double *x, double *y, double *z, SVector3 *
   y[0] = pnt1.y(); y[1] = pnt2.y(); y[2] = pnt3.y();
   z[0] = pnt1.z(); z[1] = pnt2.z(); z[2] = pnt3.z();
   
-
-
-
 }
 
-
 int MTriangleN::getNumEdgesRep(){ return 3 * numSubEdges; }
 
 void MTriangleN::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)

Geo/gmshFace.cpp

-// $Id: gmshFace.cpp,v 1.59 2008-06-10 08:37:34 remacle Exp $
+// $Id: gmshFace.cpp,v 1.60 2008-06-10 12:59:12 remacle Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -85,6 +85,7 @@ gmshFace::gmshFace(GModel *m, Surface *face)
         Msg::Error("Unknown point %d", v->Num);
     }
   }
+  isSphere = iSRuledSurfaceASphere (s,center,radius);
 }
 
 double gmshFace::getMetricEigenvalue(const SPoint2 &pt)
@@ -242,6 +243,24 @@ 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());    
+    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());
+  }
+  else{
+    return pt1 + (pt2 - pt1) * t;
+  }
+}
+
+
 int gmshFace::containsPoint(const SPoint3 &pt) const
 { 
   if(geomType() == Plane){

Geo/gmshFace.h

@@ -26,12 +26,15 @@
 class gmshFace : public GFace {
  protected:
   Surface *s;
-
+  bool isSphere;
+  SPoint3 center;
+  double radius;
  public:
   gmshFace(GModel *m, Surface *face);
   virtual ~gmshFace(){}
   Range<double> parBounds(int i) const; 
   void setModelEdges(std::list<GEdge*>&);
+  virtual SPoint2 geodesic(const SPoint2 &pt1, const SPoint2 &pt2, double t);
   virtual GPoint point(double par1, double par2) const; 
   virtual GPoint closestPoint(const SPoint3 & queryPoint, const double initialGuess[2]) const; 
   virtual int containsPoint(const SPoint3 &pt) const;  

Mesh/meshGRegion.cpp

-// $Id: meshGRegion.cpp,v 1.49 2008-06-10 08:37:34 remacle Exp $
+// $Id: meshGRegion.cpp,v 1.50 2008-06-10 12:59:12 remacle Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -188,12 +188,11 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
 	  // This is not 100 % safe yet, so we reserve that operation for high order
 	  // meshes.
 	  GPoint gp = gf->closestPoint (SPoint3(v[j]->x(), v[j]->y(), v[j]->z()),guess);
-	  //	printf("ah que coucou\n");
-	  //	const double U(0),V(0);
-	  //        MVertex *v1b = new MFaceVertex(v[j]->x(), v[j]->y(), v[j]->z(), gf,U,V);
+
 	  Msg::Debug("A new point has been inserted in mesh face %d by the 3D mesher",gf->tag());
 	  Msg::Debug("The point has been projected back to the surface (%g %g %g) -> (%g %g %g)",
 		     v[j]->x(), v[j]->y(), v[j]->z(),gp.x(),gp.y(),gp.z());
+
 	  // To be safe, we should ensure that this mesh motion does not lead to an invalid mesh !!!!
 	  v[j]->x() = gp.x();
 	  v[j]->y() = gp.y();
@@ -201,7 +200,7 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
 	  v1b = new MFaceVertex(v[j]->x(), v[j]->y(), v[j]->z(), gf,gp.u(),gp.v());
 	}
 	else{
-	  v1b = new MVertex(v[j]->x(), v[j]->y(), v[j]->z());
+	  v1b = new MVertex(v[j]->x(), v[j]->y(), v[j]->z(),gf);
 	}
 
         gf->mesh_vertices.push_back(v1b);