diff --git a/Geo/FEdge.cpp b/Geo/FEdge.cpp
index 545c57024d7169bc486b60d399efac5b32c0ece8..d164625492c1a2cb72d8b2e8cd42f25782b93a59 100644
--- a/Geo/FEdge.cpp
+++ b/Geo/FEdge.cpp
@@ -1,3 +1,4 @@
+#include "Message.h"
 #include "FEdge.h"
 
 #if defined(HAVE_FOURIER_MODEL)
@@ -10,16 +11,63 @@ Range<double> FEdge::parBounds(int i) const
 GPoint FEdge::point(double p) const 
 {
   double x, y, z;
-
-  edge->F(p,x,y,z);
+  if (edge)
+    edge->F(p,x,y,z);
+  else {
+    if (edgeNum == 0) {
+      double p0, p1, q0, q1;
+      face->Inverse(v0->x(),v0->y(),v0->z(),p0,q0);
+      face->Inverse(v1->x(),v1->y(),v1->z(),p1,q1);
+      face->F(p0+p*(p1-p0),0,x,y,z);
+    }
+    else if (edgeNum == 1) {
+      double p0, p1, q0, q1;
+      face->Inverse(v0->x(),v0->y(),v0->z(),q0,p0);
+      face->Inverse(v1->x(),v1->y(),v1->z(),q1,p1);
+      face->F(1.,p0+p*(p1-p0),x,y,z);
+    }
+    else if (edgeNum == 2) {
+      double p0, p1, q0, q1;
+      face->Inverse(v0->x(),v0->y(),v0->z(),p0,q0);
+      face->Inverse(v1->x(),v1->y(),v1->z(),p1,q1);
+      face->F(p0+p*(p1-p0),1.,x,y,z);
+    }
+    else if (edgeNum == 3) {
+      double p0, p1, q0, q1;
+      face->Inverse(v0->x(),v0->y(),v0->z(),q0,p0);
+      face->Inverse(v1->x(),v1->y(),v1->z(),q1,p1);
+      face->F(0.,p0+p*(p1-p0),x,y,z);
+    }
+    else
+      Msg(INFO,"Invalid edge number.");
+  }
   return GPoint(x,y,z);
 }
 
 double FEdge::parFromPoint(const SPoint3 &pt) const
 {
-  double t;
-  edge->Inverse(pt.x(),pt.y(),pt.z(),t);
-  return t;
+  double p;
+  if (edge)
+    edge->Inverse(pt.x(),pt.y(),pt.z(),p);
+  else {
+    if ((edgeNum == 0) || (edgeNum == 2)) {
+      double p0, p1, q0, q1, q;
+      face->Inverse(v0->x(),v0->y(),v0->z(),p0,q0);
+      face->Inverse(v1->x(),v1->y(),v1->z(),p1,q1);
+      face->Inverse(pt.x(),pt.y(),pt.z(),p,q);
+      p = p0 + p * (p1 - p0);
+    }
+    else if ((edgeNum == 1) || (edgeNum == 3)) {
+      double p0, p1, q0, q1, q;
+      face->Inverse(v0->x(),v0->y(),v0->z(),q0,p0);
+      face->Inverse(v1->x(),v1->y(),v1->z(),q1,p1);
+      face->Inverse(pt.x(),pt.y(),pt.z(),q,p);
+      p = p0 + p * (p1 - p0);
+    }
+    else
+      Msg(INFO,"Invalid edge number.");
+  }
+  return p;
 }
 
 #endif
diff --git a/Geo/FEdge.h b/Geo/FEdge.h
index ba59a9378c6d6e4bbb95188c6158afc108b3d882..ce08ccdd274c3e8f8be874e737391041ce6fa0c2 100644
--- a/Geo/FEdge.h
+++ b/Geo/FEdge.h
@@ -9,13 +9,19 @@
 #if defined(HAVE_FOURIER_MODEL)
 
 #include "FM_Edge.h"
+#include "FM_Face.h"
 
 class FEdge : public GEdge {
  protected:
   FM_Edge* edge;
+  FM_Face* face;
+  int edgeNum;
  public:
-  FEdge(GModel *model, FM_Edge* edge_, int tag, GVertex *v1, GVertex *v2) : 
-    GEdge(model, tag, v1, v2), edge(edge_) {}
+  FEdge(GModel *model, FM_Edge* edge_, int tag, GVertex *v0, GVertex *v1) : 
+    GEdge(model, tag, v0, v1), edge(edge_), face(0), edgeNum(-1) {}
+  FEdge(GModel *model, FM_Face* face_, int edgeNum_, int tag, GVertex *v0, 
+	GVertex *v1) : GEdge(model, tag, v0, v1), edge(0), face(face_), 
+    edgeNum(edgeNum_) {}
   virtual ~FEdge() {}
   double period() const { throw ; }
   virtual bool periodic(int dim=0) const { return false; }
diff --git a/Geo/FFace.cpp b/Geo/FFace.cpp
index 6640975444b6e2e04202c49ec8084acad4eb9811..406675a75ad719b959c2717a10fac6ff2ae45f99 100644
--- a/Geo/FFace.cpp
+++ b/Geo/FFace.cpp
@@ -1,17 +1,56 @@
+#include <list>
 #include "FVertex.h"
 #include "FFace.h"
+#include "Message.h"
 
 #if defined(HAVE_FOURIER_MODEL)
 
 FFace::FFace(GModel *m, FM_Face *face_, int tag) : GFace(m,tag), face(face_) 
 {
-  int curCorner = 0;
-  for (int i=0;i<face->GetNumEdges();i++) {
-    GVertex* p0 = new FVertex(m,curCorner,face->GetEdge(i)->GetStartPoint());
-    curCorner = (curCorner + 1) % face->GetNumEdges();
-    GVertex* p1 = new FVertex(m,curCorner,face->GetEdge(i)->GetEndPoint());
-    l_edges.push_back(new FEdge(m,face->GetEdge(i),i,p0,p1));
-    l_dirs.push_back(1);
+  if (face->GetNumEdges()) {
+    std::list<GVertex*> corners;
+    std::list<GVertex*>::iterator itStart, itEnd;
+
+    corners.push_back(new FVertex(m,0,face->GetEdge(0)->
+				  GetStartPoint()));
+    corners.push_back(new FVertex(m,1,face->GetEdge(1)->
+				  GetStartPoint()));
+    corners.push_back(new FVertex(m,2,face->GetEdge(2)->
+				  GetStartPoint()));
+    corners.push_back(new FVertex(m,3,face->GetEdge(3)->
+				  GetStartPoint()));
+
+    itStart = itEnd = corners.begin(); itEnd++;
+    for (int i=0;i<face->GetNumEdges();i++) {
+      l_edges.push_back(new FEdge(m,face->GetEdge(i),i,*itStart,*itEnd));
+      l_dirs.push_back(1);      
+      itStart++; itEnd++;
+      if (itEnd == corners.end()) {
+	itEnd = corners.begin();
+      }
+    }
+    for (std::list<GEdge*>::iterator it = l_edges.begin();it != l_edges.end();
+	 it++) {
+      GVertex *start = (*it)->getBeginVertex();
+      GVertex *end = (*it)->getEndVertex();
+      Msg(INFO,"(%g,%g,%g) --- (%g,%g,%g)",start->x(),start->y(),start->z(),
+	  end->x(),end->y(),end->z());
+    }
+  }
+  else {
+    double x,y,z;
+    face->F(0.,0.,x,y,z);
+    GVertex* p0 = new FVertex(m,0,new FM_Vertex(x,y,z));
+    face->F(1.,0.,x,y,z);
+    GVertex* p1 = new FVertex(m,1,new FM_Vertex(x,y,z));
+    face->F(1.,1.,x,y,z);
+    GVertex* p2 = new FVertex(m,2,new FM_Vertex(x,y,z));
+    face->F(0.,1.,x,y,z);
+    GVertex* p3 = new FVertex(m,3,new FM_Vertex(x,y,z));
+    l_edges.push_back(new FEdge(m,face,0,0,p0,p1));
+    l_edges.push_back(new FEdge(m,face,1,1,p1,p2));
+    l_edges.push_back(new FEdge(m,face,2,2,p2,p3));
+    l_edges.push_back(new FEdge(m,face,3,3,p3,p0));
   }
 }
 
diff --git a/Geo/FVertex.h b/Geo/FVertex.h
index 481e16fe033451ddb057802863309065bc6831fc..14f610a5fae1ecab219fd73ab1c231374845e26d 100644
--- a/Geo/FVertex.h
+++ b/Geo/FVertex.h
@@ -35,7 +35,7 @@ class FVertex : public GVertex {
     return v->GetZ();
   }
   ModelType getNativeType() const { return FourierModel; }
-  virtual SPoint2 reparamOnFace ( GFace *gf , int) const { throw; }
+  //virtual SPoint2 reparamOnFace ( GFace *gf , int) const { throw; }
 };
 
 #endif
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 983d32f8c64c0a9792ceea1f41120842bc904294..32ef3bdd2993e1939d83726879e7f58af62ef1ad 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -1,4 +1,4 @@
-// $Id: GFace.cpp,v 1.32 2007-02-02 23:50:33 geuzaine Exp $
+// $Id: GFace.cpp,v 1.33 2007-05-05 01:04:40 anand Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -507,6 +507,8 @@ void GFace::XYZtoUV(const double X, const double Y, const double Z,
 SPoint2 GFace::parFromPoint(const SPoint3 &p) const
 {
   double U,V;
+
+  Msg(INFO,"Instead, I died here");
   
   XYZtoUV(p.x(),p.y(),p.z(),U,V,1.0);
 
diff --git a/Geo/GModelIO_F.cpp b/Geo/GModelIO_F.cpp
index c096b4a63fd4ab5351cda257383147dae28590ad..c621df437aa3f33eb5305e21ced0edc498d74350 100644
--- a/Geo/GModelIO_F.cpp
+++ b/Geo/GModelIO_F.cpp
@@ -12,20 +12,6 @@
 
 extern Context_T CTX;
 
-class meshGmsh{
-public:
-  void operator() (GFace *gf)
-  {
-    FFace *ff = dynamic_cast<FFace*>(gf);
-    if(!ff) {
-      Msg(GERROR, "face %d is not Fourier", gf->tag());
-      return;
-    }
-    meshGFace mesh;
-    mesh(ff);
-  }
-};
-
 int GModel::readF(const std::string &fn)
 {
   CTX.terminal = 1; 
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 1a10eac9d2d3fa689b333caf63f131af979c2856..9f82222034be889bbe49aa004fec2204492b8cc1 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.74 2007-04-21 22:08:30 geuzaine Exp $
+// $Id: meshGFace.cpp,v 1.75 2007-05-05 01:04:40 anand Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -44,13 +44,22 @@ void computeEdgeLoops(const GFace *gf,
 		      std::vector<int> &indices)
 {
   std::list<GEdge*> edges = gf->edges();
+  Msg(INFO,"Number of Edges = %d",edges.size());
   std::list<int> ori = gf->orientations();
+  Msg(INFO,"Number of Orientations = %d",ori.size());
   std::list<GEdge*>::iterator it = edges.begin();
   std::list<int>::iterator ito = ori.begin();
     
   indices.push_back(0);
   GVertex *start = ((*ito) == 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
   GVertex *v_end = ((*ito) != 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
+
+  Msg(INFO,"(%g,%g,%g) --- (%g,%g,%g)",start->x(),start->y(),start->z(),
+      v_end->x(),v_end->y(),v_end->z());  
+  Msg(INFO,"Mesh Size = (%d,%d)",start->mesh_vertices.size(),
+      v_end->mesh_vertices.size());
+  Msg(INFO,"Edge Mesh Size = %d",(*it)->mesh_vertices.size());
+
   all_mvertices.push_back(start->mesh_vertices[0]);
   if (*ito == 1)
     for (unsigned int i=0;i<(*it)->mesh_vertices.size();i++)
@@ -63,6 +72,7 @@ void computeEdgeLoops(const GFace *gf,
   while(1){		
     ++it;
     ++ito;
+
     if(v_end == start){
       indices.push_back(all_mvertices.size());
       if(it == edges.end ()) break;
@@ -70,9 +80,11 @@ void computeEdgeLoops(const GFace *gf,
       v_end = ((*ito) != 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
       v_start = start;
     }
-    else{	
+    else{
       if(it == edges.end ()) throw;
       v_start = ((*ito) == 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
+      Msg(INFO,"(%g,%g,%g) --- (%g,%g,%g)",v_start->x(),v_start->y(),
+	  v_start->z(),v_end->x(),v_end->y(),v_end->z());  
       if(v_start != v_end) throw;
       v_end = ((*ito) != 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
     }
@@ -647,6 +659,8 @@ bool gmsh2DMeshGenerator ( GFace *gf , bool debug = true)
 
   v_container::iterator itv = all_vertices.begin();
 
+  Msg(INFO,"all_vertices size = %d",all_vertices.size());
+
   //  FILE *fdeb = fopen("debug.dat","w");
   //  fprintf(fdeb,"surface %d\n" ,gf->tag());
   int count = 0;
@@ -656,15 +670,18 @@ bool gmsh2DMeshGenerator ( GFace *gf , bool debug = true)
     SPoint2 param;
     if(here->onWhat()->dim() == 0){
       GVertex *gv = (GVertex*)here->onWhat();
+      Msg(INFO,"gvEdgeSize = %d",gv->edges().size());
       if (gv->edges().size() == 1)
 	{
 	  GEdge *ge = *(gv->edges().begin());
 	  Range<double> bb = ge->parBounds(0);
 	  param = ge->reparamOnFace(gf, bb.low(), 1);	  
 	}
-      else
+      else {
+	Msg(INFO,"Before dyng here");
 	param = gv->reparamOnFace(gf,1);
-      
+	Msg(INFO,"param = (%g,%g)",param.x(),param.y());
+      }
     }
     else if(here->onWhat()->dim() == 1){
       GEdge *ge = (GEdge*)here->onWhat();
@@ -757,11 +774,15 @@ bool gmsh2DMeshGenerator ( GFace *gf , bool debug = true)
   m->scalingU = 1;
   m->scalingV = 1;
 
+  Msg(INFO,"doc = %d",doc.numPoints);
+
   for(int i = 0; i < doc.numPoints; i++) 
     {
       MVertex *here = (MVertex *)doc.points[i].data;
       int num = here->getNum();
 
+      Msg(INFO,"num = %d",num);
+
       double U, V;
       // This test was missing in 2.0.0 and led to the seemingly
       // random Windows/Mac slowdowns (we were passing random numbers
@@ -777,7 +798,7 @@ bool gmsh2DMeshGenerator ( GFace *gf , bool debug = true)
       }
 
       BDS_Point *pp = m->add_point ( num, U,V, gf);
-      //      printf("here->onWhat = %p dim = %d\n",here->onWhat(),here->onWhat()->dim());
+      //printf("here->onWhat = %p dim = %d\n",here->onWhat(),here->onWhat()->dim());
 
        if (here->onWhat()->dim() == 1)
  	{
@@ -811,6 +832,8 @@ bool gmsh2DMeshGenerator ( GFace *gf , bool debug = true)
   // and compute characteristic lenghts using mesh edge spacing
 
   BDS_GeomEntity CLASS_F (1,2);
+
+  Msg(INFO,"came out 0");
    
   it = edges.begin();
   while(it != edges.end())
@@ -851,7 +874,7 @@ bool gmsh2DMeshGenerator ( GFace *gf , bool debug = true)
 	++ite;
       }
   }
-
+  Msg(INFO,"came out 1");
 
   it = emb_edges.begin();
   while(it != emb_edges.end())
@@ -897,6 +920,8 @@ bool gmsh2DMeshGenerator ( GFace *gf , bool debug = true)
   m->del_point(m->find_point(-3));
   m->del_point(m->find_point(-4));
 
+  Msg(INFO,"came out 2");
+
   // start mesh generation
   //  if (CTX.mesh.algo2d == ALGO_2D_MESHADAPT)
     {
@@ -916,6 +941,8 @@ bool gmsh2DMeshGenerator ( GFace *gf , bool debug = true)
 //     outputScalarField(m->triangles, name,0);
   // fill the small gmsh structures
 
+  Msg(INFO,"came out 2.5");
+
   {
     std::set<BDS_Point*,PointLessThan>::iterator itp =  m->points.begin(); 
     while (itp != m->points.end())
@@ -931,6 +958,8 @@ bool gmsh2DMeshGenerator ( GFace *gf , bool debug = true)
       }
   }
 
+  Msg(INFO,"came out 3");
+
   {
     std::list<BDS_Face*>::iterator itt = m->triangles.begin();
     while (itt != m->triangles.end())
@@ -978,6 +1007,8 @@ bool gmsh2DMeshGenerator ( GFace *gf , bool debug = true)
   delete [] U_;
   delete [] V_;
 
+  Msg(INFO,"came out");
+
   return true;
 
 }
@@ -1525,7 +1556,8 @@ void deMeshGFace::operator() (GFace *gf)
 }
 
 void meshGFace::operator() (GFace *gf) 
-{  
+{
+  Msg(INFO,"Meshing surface %d (%s)",gf->tag(), gf->getTypeString().c_str());
   if(gf->geomType() == GEntity::DiscreteSurface) return;
   if(gf->geomType() == GEntity::BoundaryLayerSurface) return;
   if(gf->geomType() == GEntity::ProjectionSurface) return;
@@ -1548,6 +1580,8 @@ void meshGFace::operator() (GFace *gf)
 
   // temp fix until we create MEdgeLoops in gmshFace
   Msg(DEBUG1, "Generating the mesh");
+
+  Msg(INFO,"ISEMPT : %d %d",gf->edgeLoops.empty(),gf->getNativeType() == GEntity::GmshModel);
   if(gf->getNativeType() == GEntity::GmshModel || gf->edgeLoops.empty()){
     gmsh2DMeshGenerator(gf,false);
   }
diff --git a/contrib/FourierModel/FM_Edge.cpp b/contrib/FourierModel/FM_Edge.cpp
index a49f7d5fad1b44f3bb09226bed220396045ee43b..02eb2ad32a9786b63327d5dca7362379c6977a55 100644
--- a/contrib/FourierModel/FM_Edge.cpp
+++ b/contrib/FourierModel/FM_Edge.cpp
@@ -40,7 +40,7 @@ bool FM_Edge::Inverse(double x,double y,double z,double &t)
     else if (_EP->GetZ() - _SP->GetZ())
       t = (z - _SP->GetZ()) / (_EP->GetZ() - _SP->GetZ());
     else {
-      Msg::Warning("Cannnot inver the curve");
+      Msg::Warning("Cannnot invert the curve");
       t = 0.;
     }
   }
diff --git a/contrib/FourierModel/FM_Face.cpp b/contrib/FourierModel/FM_Face.cpp
index a25f26ee02230acf9c23eaba2a2be3f8b02e0d12..e38e4b3760f749b2c8f88cdb80cf5ac1dc735658 100644
--- a/contrib/FourierModel/FM_Face.cpp
+++ b/contrib/FourierModel/FM_Face.cpp
@@ -19,7 +19,7 @@ void FM_Face::F(double u, double v, double &x, double &y, double &z) {
     double xU, yU, zU, uU, vU;
     _edge[0]->F(v,xD,yD,zD);
     _patch->Inverse(xD,yD,zD,uD,vD);
-    _edge[2]->F(v,xU,yU,zU);
+    _edge[2]->F(1.-v,xU,yU,zU);
     _patch->Inverse(xU,yU,zU,uU,vU);
    
     double U = uD + u * (uU - uD);
@@ -82,7 +82,7 @@ bool FM_Face::Inverse(double x,double y,double z,double &u, double &v)
     v = (V-v0)/(v1-v0);
     _edge[0]->F(v,xD,yD,zD);
     _patch->Inverse(xD,yD,zD,uD,vD);
-    _edge[2]->F(v,xU,yU,zU);
+    _edge[2]->F(1.-v,xU,yU,zU);
     _patch->Inverse(xU,yU,zU,uU,vU);
     
     u = (U - uD) / (uU - uD);