diff --git a/Geo/GFace.h b/Geo/GFace.h
index 18c8825b870669d441804e3030e2cd9baccf9549..f482106f67bd5af4f12af029c43c8f34687baa68 100644
--- a/Geo/GFace.h
+++ b/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);
diff --git a/Geo/GeoInterpolation.cpp b/Geo/GeoInterpolation.cpp
index 571f300739d8bb368a6d0495ec762f1b4cb9837b..06ff057b61b4f46117c6c6ab19773fea921fe893 100644
--- a/Geo/GeoInterpolation.cpp
+++ b/Geo/GeoInterpolation.cpp
@@ -1,4 +1,4 @@
-// $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 ¢er, 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};
diff --git a/Geo/GeoInterpolation.h b/Geo/GeoInterpolation.h
index e878cbe6ed601887718ed0543b5e3eed502948d3..aa00d52797f7ccb3b1fa2d5fbec24fa68dae9cc3 100644
--- a/Geo/GeoInterpolation.h
+++ b/Geo/GeoInterpolation.h
@@ -22,6 +22,7 @@
#include "Geo.h"
+bool iSRuledSurfaceASphere(Surface *s, SPoint3 ¢er, 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],
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index e185b08086983898453c0e45e312644db8172020..a39f7744184c71326a2ee090fc9f2c7c81745de9 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1,4 +1,4 @@
-// $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)
diff --git a/Geo/gmshFace.cpp b/Geo/gmshFace.cpp
index 6c86c92fb930d8df7a8172eded112441e5433ab5..1f5c9b1430d9e65e83ada2913c3d0c31b735199a 100644
--- a/Geo/gmshFace.cpp
+++ b/Geo/gmshFace.cpp
@@ -1,4 +1,4 @@
-// $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){
diff --git a/Geo/gmshFace.h b/Geo/gmshFace.h
index 80f8ba5efe0de82d974796ef5d20d4292664e0c2..073caa7cbacb53ad53114e0137dc240c5f76f025 100644
--- a/Geo/gmshFace.h
+++ b/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;
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 44b84b07fd74a3713e4b9c314baec36eb905d550..c9dfebc0e83961f5208e5ba410caf21400afb28b 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -1,4 +1,4 @@
-// $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);