diff --git a/Fltk/GUI_Projection.cpp b/Fltk/GUI_Projection.cpp
index ebe0848d713dcb3ecfa79f0ebb7df8922de6bea1..da58e7fafd367dc3d5f7bf48cf902f82e7ba91d3 100644
--- a/Fltk/GUI_Projection.cpp
+++ b/Fltk/GUI_Projection.cpp
@@ -18,16 +18,17 @@ void uvPlot::draw()
   fl_rectf(0, 0, w(), h());
 
   // draw points
-  fl_color(FL_BLACK);
   if(_u.size() != _v.size()) return;
   for(unsigned int i = 0; i < _u.size(); i++){
     int x = (int)(_u[i] * w());
     int y = (int)(_v[i] * h());
+    // fixme: change color depending on f
+    fl_color(FL_BLACK);
     fl_rect(x, y, 3, 3);
   }
 
-  fl_font(FL_HELVETICA, 14);
-  fl_draw("Hello", w() / 2, h() / 2);
+  //fl_font(FL_HELVETICA, 14);
+  //fl_draw("Hello", w() / 2, h() / 2);
 }
 
 projection::projection(FProjectionFace *f, int x, int y, int w, int h, int BB, int BH, 
@@ -65,14 +66,13 @@ projection::projection(FProjectionFace *f, int x, int y, int w, int h, int BB, i
       v->step(CTX.lc / 100.);
       v->label((i == 6) ? "X translation" : (i == 7) ? "Y translation" : 
 	       "Z translation");
-      //v->value(currentParams[i]);
-      v->value(bounds.center()[i]);
+      v->value(bounds.center()[i - 6]);
     }
     else{ // other parameters
       v->maximum(10. * CTX.lc);
       v->minimum(-10. * CTX.lc);
       v->step(CTX.lc / 100.);
-      v->label("My nice label");
+      v->label(strdup(ps->GetLabel(i - 9).c_str()));
       v->value(ps->GetParameter(i - 9));
       v->value(currentParams[i]);
     }
@@ -173,67 +173,10 @@ void browse_cb(Fl_Widget *w, void *data)
     projections[i]->group->hide();
   }
   
-  /*
-<<<<<<< GUI_Projection.cpp
-  for(int i = 0; i < MAX_PROJECTION_PARAMETERS; i++)  
-    e->getValueInput(i)->hide();
-
-  FProjectionFace *f = e->getCurrentProjectionFace();
-  if(f){
-    f->setVisibility(true);
-    ProjectionSurface *ps = f->GetProjectionSurface();
-    for(int i = 0; i < 9; i++){
-      e->getValueInput(i)->show();
-      ps->GetScale(currentParams[0],currentParams[1],currentParams[2]);
-      ps->GetOrigin(currentParams[6],currentParams[7],currentParams[8]);
-      if(i < 3){ // scaling
-	e->getValueInput(i)->maximum(CTX.lc * 10.);
-	e->getValueInput(i)->minimum(CTX.lc / 100.);
-	e->getValueInput(i)->step(CTX.lc / 100.);
-	e->getValueInput(i)->
-	  label((i == 0) ? "X scale" : (i == 1) ? "Y scale" : "Z scale");
-	e->getValueInput(i)->value(currentParams[i]);
-      }
-      else if(i < 6){ //rotation
-	currentParams[i] = 0.;
-	e->getValueInput(i)->maximum(180.);
-	e->getValueInput(i)->minimum(-180.);
-	e->getValueInput(i)->step(0.1);
-	e->getValueInput(i)->
-	  label((i == 3) ? "X Rotation" : (i == 4) ? "Y Rotation" : 
-		"Z Rotation");
-	e->getValueInput(i)->value(currentParams[i]);
-      }
-      else{ // translation
-	e->getValueInput(i)->maximum(bounds.max()[i] + CTX.lc);
-	e->getValueInput(i)->minimum(bounds.min()[i] - CTX.lc);
-	e->getValueInput(i)->step(( CTX.lc) / 100.);
-	e->getValueInput(i)->
-	  label((i == 6) ? "X Translation" : (i == 7) ? "Y Translation" : 
-		"Z Translation");
-	e->getValueInput(i)->value(currentParams[i]);
-      }
-    }
-    for(int i = 9; i < 9 + ps->GetNumParameters(); i++){
-      currentParams[i] = ps->GetParameter(i - 9);
-      e->getValueInput(i)->show();
-      e->getValueInput(i)->maximum(10. * CTX.lc);
-      e->getValueInput(i)->minimum(-10. * CTX.lc);
-      e->getValueInput(i)->step(CTX.lc / 100.);
-      e->getValueInput(i)->label("My nice label");
-      e->getValueInput(i)->value(currentParams[i]);
-    }
-=======
-  */
-
   projection *p = e->getCurrentProjection();
   if(p){
     p->face->setVisibility(true);
     p->group->show();
-
-
-    // >>>>>>> 1.12
-
   }
   Draw();
 }
@@ -262,15 +205,29 @@ void update_cb(Fl_Widget *w, void *data)
     Draw();
     for (int i = 0; i < 9; i++)
       currentParams[i] = p->parameters[i]->value();
-}
-
-  // project all selected points and update u,v display
-  std::vector<double> u, v;
-  for(int i = 0; i < 5000; i++){
-    u.push_back((double)rand() / (double)RAND_MAX);
-    v.push_back((double)rand() / (double)RAND_MAX);
+    
+    // project all selected points and update u,v display
+    std::vector<double> u, v, f;
+    std::vector<GEntity*> &ent(e->getEntities());
+    for(unsigned int i = 0; i < ent.size(); i++){
+      if(ent[i]->getSelection()){
+	GVertex *ve = dynamic_cast<GVertex*>(ent[i]);
+	if(!ve)
+	  Msg(GERROR, "Problem in point selection processing");
+	else{
+	  double uu, vv, xx, yy, zz;
+	  ps->OrthoProjectionOnSurface(ve->x(),ve->y(),ve->z(),uu,vv);
+	  u.push_back(uu);
+	  v.push_back(vv);
+	  ps->F(uu,vv,xx,yy,zz);
+	  double dx = xx - ve->x(), dy= yy - ve->y(), dz = zz - ve->z();
+	  f.push_back(sqrt(dx * dx + dy * dy + dz * dz));
+	}
+      }
+    }
+    // loop over elements and do the same thing
+    e->uv()->fill(u, v, f);
   }
-  e->uv()->fill(u, v);
 }
 
 void select_cb(Fl_Widget *w, void *data)
@@ -364,6 +321,7 @@ void select_cb(Fl_Widget *w, void *data)
       ZeroHighlight();
       break;
     }
+    update_cb(w, data);
   }
 
   CTX.mesh.changed = ENT_ALL;
@@ -403,7 +361,119 @@ void compute_cb(Fl_Widget *w, void *data)
 {
   projectionEditor *e = (projectionEditor*)data;
 
-  Msg(GERROR, "Compute!");
+  projection *p = e->getCurrentProjection();
+  if(p){
+    ProjectionSurface *ps = p->face->GetProjectionSurface();
+
+    // project all selected points and update u,v display
+    std::vector<double> u, v;
+    std::vector<std::complex<double> > f;
+    std::vector<GEntity*> &ent(e->getEntities());
+    for(unsigned int i = 0; i < ent.size(); i++){
+      GVertex *ve = dynamic_cast<GVertex*>(ent[i]);
+      if(!ve){
+	Msg(GERROR, "Problem in point selection processing");
+      }
+      else{
+	double uu, vv;
+	ps->OrthoProjectionOnSurface(ve->x(),ve->y(),ve->z(),uu,vv);
+	u.push_back(uu);
+	v.push_back(vv);
+	double p[3], n[3];
+	ps->F(u[i],v[i],p[0],p[1],p[2]);
+	ps->GetUnitNormal(u[i],v[i],n[0],n[1],n[2]);
+	f.push_back((ve->x() - p[0]) * n[0] + (ve->y() - p[1]) * n[1] + 
+		    (ve->z() - p[2]) * n[2]);
+      }
+    }
+    Patch* patch = new FPatch(0,ps,u,v,f,3);
+    patch->SetMinU(0.0);
+    patch->SetMaxU(0.5);
+
+    double LL[2], LR[2], UL[2], UR[2];
+    LL[0] = 0.0; LL[1] = 0.0;
+    LR[0] = 1.0; LR[1] = 0.0;
+    UL[0] = 0.0; UL[1] = 1.0;
+    UR[0] = 1.0; UR[1] = 1.0;
+
+    int i1, i2;
+    double xx,yy,zz;
+
+    int tagVertex = GMODEL->numVertex();
+    patch->F(LL[0],LL[1],xx,yy,zz);
+    FM_Vertex* vLL = new FM_Vertex(++tagVertex,xx,yy,zz);
+    GMODEL->add(new FVertex(GMODEL,vLL->GetTag(),vLL));
+    patch->F(LR[0],LR[1],xx,yy,zz);
+    FM_Vertex* vLR = new FM_Vertex(++tagVertex,xx,yy,zz);
+    GMODEL->add(new FVertex(GMODEL,vLR->GetTag(),vLR));
+    patch->F(UL[0],UL[1],xx,yy,zz);
+    FM_Vertex* vUL = new FM_Vertex(++tagVertex,xx,yy,zz);
+    GMODEL->add(new FVertex(GMODEL,vUL->GetTag(),vUL));
+    patch->F(UR[0],UR[1],xx,yy,zz);
+    FM_Vertex* vUR = new FM_Vertex(++tagVertex,xx,yy,zz);
+    GMODEL->add(new FVertex(GMODEL,vUR->GetTag(),vUR));
+
+    Curve* curveB = new FCurve(0,patch,LL,LR);
+    Curve* curveR = new FCurve(0,patch,LR,UR);
+    Curve* curveT = new FCurve(0,patch,UR,UL);
+    Curve* curveL = new FCurve(0,patch,UL,LL);
+
+    int tagEdge = GMODEL->numEdge();
+    FM_Edge* eB = new FM_Edge(++tagEdge,curveB,vLL,vLR);
+    i1 = eB->GetStartPoint()->GetTag();
+    i2 = eB->GetEndPoint()->GetTag();
+    GMODEL->add(new FEdge(GMODEL,eB,eB->GetTag(),GMODEL->vertexByTag(i1),
+			  GMODEL->vertexByTag(i2)));
+    FM_Edge* eR = new FM_Edge(++tagEdge,curveR,vLR,vUR); 
+    i1 = eR->GetStartPoint()->GetTag();
+    i2 = eR->GetEndPoint()->GetTag();
+    GMODEL->add(new FEdge(GMODEL,eR,eR->GetTag(),GMODEL->vertexByTag(i1),
+			  GMODEL->vertexByTag(i2))); 
+    FM_Edge* eT = new FM_Edge(++tagEdge,curveT,vUR,vUL);
+    i1 = eT->GetStartPoint()->GetTag();
+    i2 = eT->GetEndPoint()->GetTag();
+    GMODEL->add(new FEdge(GMODEL,eT,eT->GetTag(),GMODEL->vertexByTag(i1),
+			  GMODEL->vertexByTag(i2)));
+    FM_Edge* eL = new FM_Edge(++tagEdge,curveL,vUL,vLL); 
+    i1 = eL->GetStartPoint()->GetTag();
+    i2 = eL->GetEndPoint()->GetTag();
+    GMODEL->add(new FEdge(GMODEL,eL,eL->GetTag(),GMODEL->vertexByTag(i1),
+			  GMODEL->vertexByTag(i2)));
+
+    FM_Face* face = new FM_Face(GMODEL->numFace() + 1,patch);
+    face->AddEdge(eB); face->AddEdge(eR); face->AddEdge(eT); face->AddEdge(eL);
+    std::list<GEdge*> l_edges;
+    for (int j=0;j<face->GetNumEdges();j++) {
+      int tag = face->GetEdge(j)->GetTag(); 
+      l_edges.push_back(GMODEL->edgeByTag(tag));
+    }
+    GMODEL->add(new FFace(GMODEL,face,face->GetTag(),l_edges));
+
+    int nU=64;
+    int nV=64;
+    
+    std::vector<int> color(3);
+    
+    std::vector<std::vector<double> > x(nU,std::vector<double>(nV));
+    std::vector<std::vector<double> > y(nU,std::vector<double>(nV));
+    std::vector<std::vector<double> > z(nU,std::vector<double>(nV));
+    
+    std::vector<std::vector<int> > mask = ones(nU,nV);
+    
+    double hU = 1./(nU-1);
+    double hV = 1./(nV-1);
+    
+    for (int j=0;j<nU;j++) {
+      for (int k=0;k<nV;k++) {
+	double u = j*hU;
+	double v = k*hV;
+	patch->F(u,v,x[j][k],y[j][k],z[j][k]);
+      }
+    }
+    color[0] = 0; color[1] = 0; color[2] = 1;
+    plotSceneViewer(0,"patch.iv",color,x,y,z,nU,nV,mask);
+  }
+  //Msg(GERROR, "Compute!");
 }
 
 void mesh_parameterize_cb(Fl_Widget* w, void* data)
@@ -414,10 +484,10 @@ void mesh_parameterize_cb(Fl_Widget* w, void* data)
   // create one instance of each available projection surface
   std::vector<FProjectionFace*> faces;
   if(faces.empty()){
-    faces.push_back(new FProjectionFace
-		    (GMODEL, 10000, new CylindricalProjectionSurface(0)));
-    faces.push_back(new FProjectionFace
-		    (GMODEL, 10001, new RevolvedParabolaProjectionSurface(1)));
+    faces.push_back(new FProjectionFace(GMODEL, GMODEL->numFace() + 1, 
+					new CylindricalProjectionSurface(0)));
+    faces.push_back(new FProjectionFace(GMODEL, GMODEL->numFace() + 1,
+					new RevolvedParabolaProjectionSurface(1)));
   }
 
   // make each projection surface invisible and 
diff --git a/Fltk/GUI_Projection.h b/Fltk/GUI_Projection.h
index ba90f70d0b1a78588666ef5834abbeb37daef147..b461d68173202384498af601efe1cd960fa31047 100644
--- a/Fltk/GUI_Projection.h
+++ b/Fltk/GUI_Projection.h
@@ -8,11 +8,20 @@
 #include "Shortcut_Window.h"
 #include <FL/Fl_Toggle_Button.H>
 #include <FL/Fl_Round_Button.H>
-
 #if defined(HAVE_FOURIER_MODEL)
 
+#include "Utils.h"
 #include <vector>
+#include <complex>
+#include "FVertex.h"
+#include "FEdge.h"
+#include "FFace.h"
+#include "FPatch.h"
+#include "FCurve.h"
 #include "FProjectionFace.h"
+#include "FM_Vertex.h"
+#include "FM_Edge.h"
+#include "FM_Face.h"
 
 void select_cb(Fl_Widget *w, void *data);
 void browse_cb(Fl_Widget *w, void *data);
@@ -24,11 +33,14 @@ void compute_cb(Fl_Widget *w, void *data);
 
 class uvPlot : public Fl_Window {
  private:
-  std::vector<double> _u, _v;
+  std::vector<double> _u, _v, _f;
   void draw();
  public:
   uvPlot(int x, int y, int w, int h, const char *l=0) : Fl_Window(x, y, w, h, l){}
-  void fill(std::vector<double> &u, std::vector<double> &v){ _u = u; _v = v; redraw(); }
+  void fill(std::vector<double> &u, std::vector<double> &v, std::vector<double> &f)
+  { 
+    _u = u; _v = v; _f = f; redraw();
+  }
 };
 
 class projectionEditor;
diff --git a/Geo/FFace.cpp b/Geo/FFace.cpp
index 5f391f5742962cfa73cf4693d6bd01b69a1827a2..57c1ab85d55e60d18b96b9f4871952c7bde3dd5e 100644
--- a/Geo/FFace.cpp
+++ b/Geo/FFace.cpp
@@ -15,6 +15,17 @@ FFace::FFace(GModel *m, FM_Face *face_, int tag, std::list<GEdge*> l_edges_)
   }
 }
 
+FFace::FFace(GModel *m, FM_Face *face_, int tag, std::list<GEdge*> l_edges_,
+	     std::list<int> l_dirs_) : GFace(m,tag), face(face_)
+{
+  for (std::list<GEdge*>::iterator it = l_edges_.begin();
+       it != l_edges_.end(); it++)
+    l_edges.push_back((*it));  
+  for (std::list<int>::iterator it = l_dirs_.begin();
+       it != l_dirs_.end(); it++)
+    l_dirs.push_back((*it));  
+}
+
 Range<double> FFace::parBounds(int i) const
 {
   return Range<double>(0.,1.);
@@ -38,23 +49,7 @@ SPoint2 FFace::parFromPoint(const SPoint3 &p) const
   double u, v, x, y, z;
   x = p.x(); y = p.y(); z = p.z();
   face->Inverse(x,y,z,u,v);
-  printf("per : %d %d\n",face->IsPeriodic(0),face->IsPeriodic(1));
-  if (face->IsPeriodic(0)) {
-    double s,e;
-    printf("\tper : %d %d\n",face->GetEdge(0)->GetCurveExtent(s,e),
-	   face->GetEdge(1)->GetCurveExtent(s,e));
-    if (face->GetEdge(0)->GetCurveExtent(s,e)) {
-      printf("%g %g %g\n",u,s,e);
-      u -= floor(u - s);
-      printf("%g %g %g\n\n",u,s,e);
-    }
-  }
-  if (face->IsPeriodic(1)) {
-    double s,e;
-    if (face->GetEdge(1)->GetCurveExtent(s,e)) {
-      v -= floor(v - s);
-    } 
-  }
+
   return SPoint2(u, v);
 }
 
diff --git a/Geo/FFace.h b/Geo/FFace.h
index 452d26142d2c6efdfd45f240ea06be4a4f11ac55..d320928ff3ea3434a2622671944e3e01d4b4da5f 100644
--- a/Geo/FFace.h
+++ b/Geo/FFace.h
@@ -15,6 +15,8 @@ class FFace : public GFace {
   FM_Face *face;
  public:
   FFace(GModel *m, FM_Face *face_, int tag, std::list<GEdge*> l_edges_);
+  FFace(GModel *m, FM_Face *face_, int tag, std::list<GEdge*> l_edges_,
+	std::list<int> l_dirs_);
 
   virtual ~FFace() {}
 
diff --git a/Graphics/Geom.cpp b/Graphics/Geom.cpp
index dfe90a6d644801a8d70cf229623c2349a594b4f0..1bbf5ce4746a2c172655835e2cc2e6df349a679e 100644
--- a/Graphics/Geom.cpp
+++ b/Graphics/Geom.cpp
@@ -1,4 +1,4 @@
-// $Id: Geom.cpp,v 1.131 2007-07-26 16:28:27 anand Exp $
+// $Id: Geom.cpp,v 1.132 2007-08-02 20:55:46 anand Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -249,7 +249,7 @@ class drawGFace {
     Range<double> vbounds = f->parBounds(1);
     double umin = ubounds.low(), umax = ubounds.high();
     double vmin = vbounds.low(), vmax = vbounds.high();
-    const int N = 15;
+    const int N = 25;
     
     glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
     //glEnable(GL_LIGHTING);
@@ -285,6 +285,48 @@ class drawGFace {
     glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
   }
 
+  void _drawParametricGFace(GFace *f)
+  {
+    Range<double> ubounds = f->parBounds(0);
+    Range<double> vbounds = f->parBounds(1);
+    double umin = ubounds.low(), umax = ubounds.high();
+    double vmin = vbounds.low(), vmax = vbounds.high();
+    const int N = 40;
+    
+    glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
+    glEnable(GL_LIGHTING);
+    glLightModelf(GL_LIGHT_MODEL_TWO_SIDE, GL_TRUE);
+    glColor4ubv((GLubyte *) & CTX.color.geom.selection);
+    glBegin(GL_QUADS);
+    for(int i = 1; i < N; i++){
+      for(int j = 1; j < N; j++){
+    	double u1 = umin + (double)i/(double)(N-1) * (umax - umin);
+	double v1 = vmin + (double)j/(double)(N-1) * (vmax - vmin);
+    	double u2 = umin + (double)(i-1)/(double)(N-1) * (umax - umin);
+	double v2 = vmin + (double)(j-1)/(double)(N-1) * (vmax - vmin);
+	GPoint p1 = f->point(u1, v1);
+	GPoint p2 = f->point(u2, v1);
+	GPoint p3 = f->point(u2, v2);
+	GPoint p4 = f->point(u1, v2);
+	SVector3 n1 = f->normal(SPoint2(u1, v1));
+	SVector3 n2 = f->normal(SPoint2(u2, v1));
+	SVector3 n3 = f->normal(SPoint2(u2, v2));
+	SVector3 n4 = f->normal(SPoint2(u1, v2));
+	glNormal3d(n1.x(), n1.y(), n1.z());
+	glVertex3d(p1.x(), p1.y(), p1.z());
+	glNormal3d(n2.x(), n2.y(), n2.z());
+	glVertex3d(p2.x(), p2.y(), p2.z());
+	glNormal3d(n3.x(), n3.y(), n3.z());
+	glVertex3d(p3.x(), p3.y(), p3.z());
+	glNormal3d(n4.x(), n4.y(), n4.z());
+	glVertex3d(p4.x(), p4.y(), p4.z());
+      }
+    }
+    glEnd();
+    glDisable(GL_LIGHTING);
+    glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
+  }
+
   void _drawPlaneGFace(GFace *f)
   {
     // We create data here and the routine is not designed to be
@@ -410,6 +452,8 @@ public :
       _drawPlaneGFace(f);
     else if(f->geomType() == GEntity::ProjectionFace)
       _drawProjectionGFace(f);
+    if(f->geomType() == GEntity::ParametricSurface)
+      _drawParametricGFace(f);
     else
       _drawNonPlaneGFace(f);
     
diff --git a/contrib/FourierModel/BlendOperator.cpp b/contrib/FourierModel/BlendOperator.cpp
index acaafd96013a5e06e53a579c9f8b92c538b0f1df..da4758f9257bce71c157c42b3024818904ca95e4 100755
--- a/contrib/FourierModel/BlendOperator.cpp
+++ b/contrib/FourierModel/BlendOperator.cpp
@@ -111,9 +111,9 @@ bool BlendOperator::GetPointOnPatch_
     else {
       double QuPlus = Qu + h;
       double QvPlus = Qv + h;
-      if (GetPatch(j)->_PI->periodic[0] == 1)
+      if (GetPatch(j)->IsUPeriodic())
 	QuPlus -= std::floor(QuPlus);
-      if (GetPatch(j)->_PI->periodic[1] == 1)
+      if (GetPatch(j)->IsVPeriodic())
 	QvPlus -= std::floor(QvPlus);
 
       double QplusU[3], QplusV[3];
@@ -154,9 +154,9 @@ bool BlendOperator::GetPointOnPatch_
 
       //printf("Q1 = (%g,%g)\n",Qu,Qv);
 
-      if (GetPatch(j)->_PI->periodic[0] == 1)
+      if (GetPatch(j)->IsUPeriodic())
 	Qu -= std::floor(Qu);
-      if (GetPatch(j)->_PI->periodic[1] == 1)
+      if (GetPatch(j)->IsVPeriodic())
 	Qv -= std::floor(Qv);
 
       //printf("Q2 = (%g,%g)\n",Qu,Qv);
@@ -391,9 +391,9 @@ bool BlendOperator::GetPointOnPatch
     else {
       double QuPlus = Qu + h;
       double QvPlus = Qv + h;
-      if (GetPatch(patchTag)->_PI->periodic[0] == 1)
+      if (GetPatch(patchTag)->IsUPeriodic())
 	QuPlus -= std::floor(QuPlus);
-      if (GetPatch(patchTag)->_PI->periodic[1] == 1)
+      if (GetPatch(patchTag)->IsVPeriodic())
 	QvPlus -= std::floor(QvPlus);
       
       double QplusU[3], QplusV[3];
diff --git a/contrib/FourierModel/BlendedPatch.cpp b/contrib/FourierModel/BlendedPatch.cpp
index ecfeb0e323c64d4673239b3d04fa1e442c1309aa..30f42f1d1521afc4909500b04bf06d0b6fb9691c 100755
--- a/contrib/FourierModel/BlendedPatch.cpp
+++ b/contrib/FourierModel/BlendedPatch.cpp
@@ -13,8 +13,8 @@ BlendedPatch::BlendedPatch
 
   _derivative = patch_->GetDerivativeBitField();
 
-  _uM = patch_->GetUModes();
-  _vM = patch_->GetVModes();
+  //_uM = patch_->GetUModes();
+  //_vM = patch_->GetVModes();
  
   _periodU = patch_->IsUPeriodic() ? 1. : 2.;
   _periodV = patch_->IsVPeriodic() ? 1. : 2.;
@@ -860,16 +860,16 @@ double  BlendedPatch::GetPou(double u, double v)
 {
   double pouU, pouV;
 
-  if (patch_->_PI->hardEdge[3])
+  if (patch_->IsHardEdge(3))
     pouU = OneSidedPartitionOfUnity(0.,1.,u);
-  else if (patch_->_PI->hardEdge[1])
+  else if (patch_->IsHardEdge(1))
     pouU = 1. - OneSidedPartitionOfUnity(0.,1.,u);
   else
     pouU = PartitionOfUnity(u, 0., 0.3, 0.7, 1.);
 
-  if (patch_->_PI->hardEdge[0])
+  if (patch_->IsHardEdge(0))
     pouV = OneSidedPartitionOfUnity(0.,1.,v);
-  else if (patch_->_PI->hardEdge[2])
+  else if (patch_->IsHardEdge(2))
     pouV = 1. - OneSidedPartitionOfUnity(0.,1.,v);
   else
     pouV = PartitionOfUnity(v, 0., 0.3, 0.7, 1.);
diff --git a/contrib/FourierModel/ContinuationPatch.cpp b/contrib/FourierModel/ContinuationPatch.cpp
index 1d6e982ae74be1917432013cac5a30efb97eef72..1200e4df9ee0ef9e6d4464f3afa2dc981eb661bc 100644
--- a/contrib/FourierModel/ContinuationPatch.cpp
+++ b/contrib/FourierModel/ContinuationPatch.cpp
@@ -4,24 +4,32 @@
 
 ContinuationPatch::ContinuationPatch
 (PatchInfo* PI, ProjectionSurface* ps, int derivative) 
-  : Patch(PI),_coeffOriginalData(0),_coeffData(0),_coeffDerivU(0),
+  : Patch(),_coeffOriginalData(0),_coeffData(0),_coeffDerivU(0),
     _coeffDerivV(0),_coeffDerivUU(0),_coeffDerivVV(0),_coeffDerivUV(0)
 {
-  SetProjectionSurface(ps);
-
+  _PI = PI;
+  _tag = _PI->tag;
   _derivative = derivative;
 
-  _uModes = _PI->nModes[0];
-  _vModes = _PI->nModes[1];
+  SetProjectionSurface(ps);
 
-  _uM = _PI->nM[0];
-  _vM = _PI->nM[1];
+  _uMin = _PI->uMin;
+  _uMax = _PI->uMax;
+  _vMin = _PI->vMin;
+  _vMax = _PI->vMax;
 
-  _vM = 16;
+  _periodicityU = _PI->periodic[0];
+  _periodicityV = _PI->periodic[1];
+
+  for (int i = 0; i < 4; i++)
+    _hardEdge[i] = _PI->hardEdge[i];
 
   _periodU = (1-_PI->periodic[0]) + 1;
   _periodV = (1-_PI->periodic[1]) + 1;
 
+  _uModes = _PI->nModes[0];
+  _vModes = _PI->nModes[1];
+
   if ((_uModes % 2) == 0) {
     _uModesLower = - _uModes/2;
     _uModesUpper = _uModes/2 - 1;
@@ -39,6 +47,9 @@ ContinuationPatch::ContinuationPatch
     _vModesUpper = (_vModes - 1)/2;
   }
 
+  _uM = _PI->nM[0];
+  _vM = _PI->nM[1];
+
   if (IsUPeriodic()) {
     if ((_uM % 2) == 0) {
       _uMLower = - _uM/2;
@@ -882,16 +893,16 @@ double  ContinuationPatch::GetPou(double u, double v)
 {
   double pouU, pouV;
 
-  if (_PI->hardEdge[3])
+  if (IsHardEdge(3))
     pouU = OneSidedPartitionOfUnity(0.,1.,u);
-  else if (_PI->hardEdge[1])
+  else if (IsHardEdge(1))
     pouU = 1. - OneSidedPartitionOfUnity(0.,1.,u);
   else
     pouU = PartitionOfUnity(u, 0., 0.3, 0.7, 1.);
 
-  if (_PI->hardEdge[0])
+  if (IsHardEdge(0))
     pouV = OneSidedPartitionOfUnity(0.,1.,v);
-  else if (_PI->hardEdge[2])
+  else if (IsHardEdge(2))
     pouV = 1. - OneSidedPartitionOfUnity(0.,1.,v);
   else
     pouV = PartitionOfUnity(v, 0., 0.3, 0.7, 1.);
diff --git a/contrib/FourierModel/ContinuationPatch.h b/contrib/FourierModel/ContinuationPatch.h
index 9ad2c8221526b9fa36859201ff30206e9bad2c71..9ab90c60972847cb3d04626612a807db7bfddaa6 100644
--- a/contrib/FourierModel/ContinuationPatch.h
+++ b/contrib/FourierModel/ContinuationPatch.h
@@ -3,6 +3,7 @@
 
 #include <fftw3.h>
 #include "Patch.h"
+#include "FM_Info.h"
 #include "PartitionOfUnity.h"
 #include "Interpolator1D.h"
 
@@ -53,12 +54,10 @@ class ContinuationPatch : public Patch {
     (PatchInfo* PI, ProjectionSurface* ps, int derivative = 0);
   virtual ~ContinuationPatch();
 
-  // These are the virtual functions that must be provided by all
-  // derived patches: GetPou() returns the original smooth
-  // (non-normalized) cutoff on the patch; F() and Inverse() implement
-  // the mapping f: (u,v)->(x,y,z) and its inverse; and the Df*() and Dn*()
-  // functions return the derivatives of the mapping f and unit normal n 
-  // with respect to u and v
+  PatchInfo* _PI;
+
+  // Abstract functions of Patch
+
   virtual double GetPou(double u, double v);
   virtual void F(double u, double v, double &x, double &y, double &z);
   virtual bool Inverse(double x,double y,double z,double &u,double &v);
diff --git a/contrib/FourierModel/Curve.cpp b/contrib/FourierModel/Curve.cpp
index 47ac5b301d10102b76197806d5849abe354b8bab..7c9e419fd047e2e803a625052a5675567a774e08 100644
--- a/contrib/FourierModel/Curve.cpp
+++ b/contrib/FourierModel/Curve.cpp
@@ -1,693 +1,11 @@
-#include <math.h>
-#include "Message.h"
 #include "Curve.h"
 
-Curve::Curve(IntersectionInfo* II, std::vector<Patch*> patches) 
-  : _II(II)
+Curve::Curve() 
 {
-  _along = _II->along;
-  _h = 1.e-6;
-  _tol = 1.e-12;
-
-  if (_II->intersectingPatches[0].patchTag < 0)
-    _patch0 = 0;
-  else {
-    for (int i=0;i<patches.size();i++) {
-      if (patches[i]->GetTag() == _II->intersectingPatches[0].patchTag) {
-	_patch0 = patches[i];
-	break;
-      }
-    }
-  }
-  if (_II->intersectingPatches[1].patchTag < 0)
-    _patch1 = 0;
-  else {
-    for (int i=0;i<patches.size();i++) {
-      if (patches[i]->GetTag() == _II->intersectingPatches[1].patchTag) {
-	_patch1 = patches[i];
-	break;
-      }
-    }
-  }
-
-  if ((_patch0) && (_patch1)) {
-    double u,v,x,y,z;
-    x = _II->SP[0]; y = _II->SP[1]; z = _II->SP[2];
-    _patch0->Inverse(x,y,z,u,v);
-    _SP0[0] = u;
-    _SP0[1] = v;
-    _patch1->Inverse(x,y,z,u,v);
-    _SP1[0] = u;
-    _SP1[1] = v;
-    //printf("SPx = %g; SPy = %g; SPz = %g;\n",x,y,z);
-    
-    x = _II->EP[0]; y = _II->EP[1]; z = _II->EP[2];
-    _patch0->Inverse(x,y,z,u,v);
-    _EP0[0] = u;
-    _EP0[1] = v;
-    _patch1->Inverse(x,y,z,u,v);
-    _EP1[0] = u;
-    _EP1[1] = v;
-    //printf("EPx = %g; EPy = %g; EPz = %g;\n",x,y,z);
-  }
-  else if (_patch0) {
-    double u,v,x,y,z;
-    x = _II->SP[0]; y = _II->SP[1]; z = _II->SP[2];
-    _patch0->Inverse(x,y,z,u,v);
-    _SP0[0] = u;
-    _SP0[1] = v;
-    x = _II->EP[0]; y = _II->EP[1]; z = _II->EP[2];
-    _patch0->Inverse(x,y,z,u,v);
-    _EP0[0] = u;
-    _EP0[1] = v;
-  }
-  else if (_patch1) {
-    double u,v,x,y,z;
-    x = _II->SP[0]; y = _II->SP[1]; z = _II->SP[2];
-    _patch1->Inverse(x,y,z,u,v);
-    _SP1[0] = u;
-    _SP1[1] = v;
-    x = _II->EP[0]; y = _II->EP[1]; z = _II->EP[2];
-    _patch1->Inverse(x,y,z,u,v);
-    _EP1[0] = u;
-    _EP1[1] = v;
-  }
-  else
-    Msg::Error("No patches to find the intersection curve.");
-}
-
-void Curve::F(double t, double &x, double &y, double &z)
-{
-  if ((_patch0) && (_patch1)) {
-    double u0, v0, u1, v1;
-    //printf("patch0 : %d\n",_patch0->GetTag());
-    //printf("patch1 : %d\n",_patch1->GetTag());
-    //printf("%g,%g,%g\n",_II->SP[0],_II->SP[1],_II->SP[2]);
-    //printf("%g,%g,%g\n",_II->EP[0],_II->EP[1],_II->EP[2]);
-    //printf("%g,%g,%g,%g\n",_SP0[0],_SP0[1],_EP0[0],_EP0[1]);
-    //printf("%g,%g,%g,%g\n",_SP1[0],_SP1[1],_EP1[0],_EP1[1]);
-
-    if (_patch0->_PI->periodic[0] == 1) {
-      if (std::abs(_SP0[0]-_EP0[0])<1.e-12)
-	u0 = _SP0[0] + t * (1. + _EP0[0] - _SP0[0]);
-      else
-	u0 = _SP0[0] + t * (_EP0[0] - _SP0[0]);
-      u0 -= floor(u0);
-    }
-    else
-      u0 = _SP0[0] + t * (_EP0[0] - _SP0[0]);
-    
-    if (_patch0->_PI->periodic[1] == 1) {
-      if (std::abs(_SP0[1]-_EP0[1])<1.e-12)
-	v0 = _SP0[1] + t * (1. + _EP0[1] - _SP0[1]);
-      else
-	v0 = _SP0[1] + t * (_EP0[1] - _SP0[1]);
-      v0 -= floor(v0);
-    }
-    else
-      v0 = _SP0[1] + t * (_EP0[1] - _SP0[1]);
-    
-    if (_patch1->_PI->periodic[0] == 1) {
-      if (std::abs(_SP1[0]-_EP1[0])<1.e-12)
-	u1 = _SP1[0] + t * (1. + _EP1[0] - _SP1[0]);
-      else
-	u1 = _SP1[0] + t * (_EP1[0] - _SP1[0]);
-      u1 -= floor(u1);
-    }
-    else
-      u1 = _SP1[0] + t * (_EP1[0] - _SP1[0]);
-    
-    if (_patch1->_PI->periodic[1] == 1) {
-      if (std::abs(_SP1[1]-_EP1[1])<1.e-12)
-	v1 = _SP1[1] + t * (1. + _EP1[1] - _SP1[1]);
-      else
-	v1 = _SP1[1] + t * (_EP1[1] - _SP1[1]);
-      v1 -= floor(v1);
-    }
-    else
-      v1 = _SP1[1] + t * (_EP1[1] - _SP1[1]);
-    
-    if (_along) {
-      double x0, y0, z0;
-      _patch0->F(u0,v0,x0,y0,z0);
-      
-      double x1, y1, z1;
-      _patch1->F(u1,v1,x1,y1,z1);
-      
-      double r,u,v,rPlus,uPlus,vPlus;
-      
-      r = u0;
-      u = u1;
-      v = v1;
-      
-      double F[3];
-      F[0] = x1 - x0;
-      F[1] = y1 - y0;
-      F[2] = z1 - z0;
-      
-      while ((std::abs(F[0])>_tol) || (std::abs(F[1])>_tol) ||
-	     (std::abs(F[2])>_tol)) {
-	rPlus = r + _h;
-	uPlus = u + _h;
-	vPlus = v + _h;
-	if (_patch0->_PI->periodic[0] == 1)
-	  rPlus -= std::floor(rPlus);
-	if (_patch1->_PI->periodic[0] == 1)
-	  uPlus -= std::floor(uPlus);
-	if (_patch1->_PI->periodic[1] == 1)
-	  vPlus -= floor(vPlus);
-	
-	double x0rPlus, y0rPlus, z0rPlus;
-	_patch0->F(rPlus,v0,x0rPlus,y0rPlus,z0rPlus);
-	double x1uPlus, y1uPlus, z1uPlus;
-	_patch1->F(uPlus,v,x1uPlus,y1uPlus,z1uPlus);
-	double x1vPlus, y1vPlus, z1vPlus;
-	_patch1->F(u,vPlus,x1vPlus,y1vPlus,z1vPlus);
-	
-	double Df[3][3];
-	Df[0][0] = - (x0rPlus - x0) / _h;
-	Df[0][1] = (x1uPlus - x1) / _h;
-	Df[0][2] = (x1vPlus - x1) / _h;
-	Df[1][0] = - (y0rPlus - y0) / _h;
-	Df[1][1] = (y1uPlus - y1) / _h;
-	Df[1][2] = (y1vPlus - y1) / _h;
-	Df[2][0] = - (z0rPlus - z0) / _h;
-	Df[2][1] = (z1uPlus - z1) / _h;
-	Df[2][2] = (z1vPlus - z1) / _h;
-	
-	double det = 
-	  Df[0][0] * (Df[1][1] * Df[2][2] - Df[1][2] * Df[2][1]) +
-	  Df[0][1] * (Df[1][2] * Df[2][0] - Df[1][0] * Df[2][2]) +
-	  Df[0][2] * (Df[1][0] * Df[2][1] - Df[1][1] * Df[2][0]);
-	
-	double update[3];
-	update[0] = 
-	  (Df[1][1] * Df[2][2] - Df[1][2] * Df[2][1]) * F[0] +
-	  (Df[0][2] * Df[2][1] - Df[0][1] * Df[2][2]) * F[1] +
-	  (Df[0][1] * Df[1][2] - Df[0][2] * Df[1][1]) * F[2];
-	update[1] =
-	  (Df[1][2] * Df[2][0] - Df[1][0] * Df[2][2]) * F[0] +
-	  (Df[0][0] * Df[2][2] - Df[0][2] * Df[2][0]) * F[1] +
-	  (Df[0][2] * Df[1][0] - Df[0][0] * Df[1][2]) * F[2];
-	update[2] =
-	  (Df[1][0] * Df[2][1] - Df[1][1] * Df[2][0]) * F[0] +
-	  (Df[0][1] * Df[2][0] - Df[0][0] * Df[2][1]) * F[1] +
-	  (Df[0][0] * Df[1][1] - Df[0][1] * Df[1][0]) * F[2];
-	
-	r -= update[0] / det;
-	u -= update[1] / det;
-	v -= update[2] / det;
-	
-	if (_patch0->_PI->periodic[0] == 1)
-	  r -= std::floor(r);
-	if (_patch1->_PI->periodic[0] == 1)
-	  u -= std::floor(u);
-	if (_patch1->_PI->periodic[1] == 1)
-	  v -= floor(v);
-	
-	_patch0->F(r,v0,x0,y0,z0);
-	_patch1->F(u,v,x1,y1,z1);
-	
-	F[0] = x1 - x0;
-	F[1] = y1 - y0;
-	F[2] = z1 - z0;
-      }
-      x = x0; y = y0; z = z0;
-    }
-    else {      
-      double x0, y0, z0;
-      _patch0->F(u0,v0,x0,y0,z0);
-      
-      double x1, y1, z1;
-      _patch1->F(u1,v1,x1,y1,z1);
-      
-      double r,u,v,rPlus,uPlus,vPlus;
-      
-      r = v0;
-      u = u1;
-      v = v1;
-      
-      double F[3];
-      F[0] = x1 - x0;
-      F[1] = y1 - y0;
-      F[2] = z1 - z0;
-      
-      while ((std::abs(F[0])>_tol) || (std::abs(F[1])>_tol) ||
-	     (std::abs(F[2])>_tol)) {
-	rPlus = r + _h;
-	uPlus = u + _h;
-	vPlus = v + _h;
-	if (_patch0->_PI->periodic[1] == 1)
-	  rPlus -= std::floor(rPlus);
-	if (_patch1->_PI->periodic[0] == 1)
-	  uPlus -= std::floor(uPlus);
-	if (_patch1->_PI->periodic[1] == 1)
-	  vPlus -= floor(vPlus);
-	
-	double x0rPlus, y0rPlus, z0rPlus;
-	_patch0->F(u0,rPlus,x0rPlus,y0rPlus,z0rPlus);
-	double x1uPlus, y1uPlus, z1uPlus;
-	_patch1->F(uPlus,v,x1uPlus,y1uPlus,z1uPlus);
-	double x1vPlus, y1vPlus, z1vPlus;
-	_patch1->F(u,vPlus,x1vPlus,y1vPlus,z1vPlus);
-	
-	double Df[3][3];
-	Df[0][0] = - (x0rPlus - x0) / _h;
-	Df[0][1] = (x1uPlus - x1) / _h;
-	Df[0][2] = (x1vPlus - x1) / _h;
-	Df[1][0] = - (y0rPlus - y0) / _h;
-	Df[1][1] = (y1uPlus - y1) / _h;
-	Df[1][2] = (y1vPlus - y1) / _h;
-	Df[2][0] = - (z0rPlus - z0) / _h;
-	Df[2][1] = (z1uPlus - z1) / _h;
-	Df[2][2] = (z1vPlus - z1) / _h;
-	
-	double det = 
-	  Df[0][0] * (Df[1][1] * Df[2][2] - Df[1][2] * Df[2][1]) +
-	  Df[0][1] * (Df[1][2] * Df[2][0] - Df[1][0] * Df[2][2]) +
-	  Df[0][2] * (Df[1][0] * Df[2][1] - Df[1][1] * Df[2][0]);
-	
-	double update[3];
-	update[0] = 
-	  (Df[1][1] * Df[2][2] - Df[1][2] * Df[2][1]) * F[0] +
-	  (Df[0][2] * Df[2][1] - Df[0][1] * Df[2][2]) * F[1] +
-	  (Df[0][1] * Df[1][2] - Df[0][2] * Df[1][1]) * F[2];
-	update[1] =
-	  (Df[1][2] * Df[2][0] - Df[1][0] * Df[2][2]) * F[0] +
-	  (Df[0][0] * Df[2][2] - Df[0][2] * Df[2][0]) * F[1] +
-	  (Df[0][2] * Df[1][0] - Df[0][0] * Df[1][2]) * F[2];
-	update[2] =
-	  (Df[1][0] * Df[2][1] - Df[1][1] * Df[2][0]) * F[0] +
-	  (Df[0][1] * Df[2][0] - Df[0][0] * Df[2][1]) * F[1] +
-	  (Df[0][0] * Df[1][1] - Df[0][1] * Df[1][0]) * F[2];
-	
-	r -= update[0] / det;
-	u -= update[1] / det;
-	v -= update[2] / det;
-	
-	if (_patch0->_PI->periodic[1] == 1)
-	  r -= std::floor(r);
-	if (_patch1->_PI->periodic[0] == 1)
-	  u -= std::floor(u);
-	if (_patch1->_PI->periodic[1] == 1)
-	  v -= floor(v);
-	
-	_patch0->F(u0,r,x0,y0,z0);
-	_patch1->F(u,v,x1,y1,z1);
-	
-	F[0] = x1 - x0;
-	F[1] = y1 - y0;
-	F[2] = z1 - z0;
-      }
-      x = x0; y = y0; z = z0;
-    }
-  }
-  else if (_patch0) {
-    double u, v;
-    if (_II->edgeInfo.edgeType == 0) {
-      double uS = _II->edgeInfo.startValue;
-      double uE = _II->edgeInfo.endValue;
-      v = _II->edgeInfo.constValue;
-      if (_II->edgeInfo.acrossDiscontinuity) {
-	if (uS < uE)
-	  uE -= 1.;
-	else
-	  uE += 1.;
-	u = uS + t * (uE - uS);
-	u -= floor(u);
-	_patch0->F(u,v,x,y,z);
-      }
-      else {
-	u = uS + t * (uE - uS);
-	_patch0->F(u,v,x,y,z);
-      }
-    }
-    else if (_II->edgeInfo.edgeType == 1) {
-      double vS = _II->edgeInfo.startValue;
-      double vE = _II->edgeInfo.endValue;
-      u = _II->edgeInfo.constValue;
-      if (_II->edgeInfo.acrossDiscontinuity) {
-	if (vS < vE)
-	  vE -= 1.;
-	else
-	  vE += 1.;
-	v = vS + t * (vE - vS);
-	v -= floor(v);
-	_patch0->F(u,v,x,y,z);
-      }
-      else {
-	v = vS + t * (vE - vS);
-	_patch0->F(u,v,x,y,z);
-      }
-    }
-    else
-      Msg::Error("Unknown edge type : %d", _II->edgeInfo.edgeType);
-  }
-  else if (_patch1) {
-    double u, v;
-    if (_II->edgeInfo.edgeType == 0) {
-      double uS = _II->edgeInfo.startValue;
-      double uE = _II->edgeInfo.endValue;
-      v = _II->edgeInfo.constValue;
-      if (_II->edgeInfo.acrossDiscontinuity) {
-	if (uS < uE)
-	  uE -= 1.;
-	else
-	  uE += 1.;
-	u = uS + t * (uE - uS);
-	u -= floor(u);
-	_patch1->F(u,v,x,y,z);
-      }
-      else {
-	u = uS + t * (uE - uS);
-	_patch1->F(u,v,x,y,z);
-      }
-    }
-    else if (_II->edgeInfo.edgeType == 1) {
-      double vS = _II->edgeInfo.startValue;
-      double vE = _II->edgeInfo.endValue;
-      u = _II->edgeInfo.constValue;
-      if (_II->edgeInfo.acrossDiscontinuity) {
-	if (vS < vE)
-	  vE -= 1.;
-	else
-	  vE += 1.;
-	v = vS + t * (vE - vS);
-	v -= floor(v);
-	_patch1->F(u,v,x,y,z);
-      }
-      else {
-	v = vS + t * (vE - vS);
-	_patch1->F(u,v,x,y,z);
-      }
-    }
-    else
-      Msg::Error("Unknown edge type : %d", _II->edgeInfo.edgeType);
-  }
-  /*
-  else if (_patch0) {
-    if (_II->edgeNumber == 0)
-      _patch0->F(t,0.,x,y,z);
-    else if (_II->edgeNumber == 1)
-      _patch0->F(1.,t,x,y,z);
-    else if (_II->edgeNumber == 2)
-      _patch0->F(1.-t,1.,x,y,z);
-    else if (_II->edgeNumber == 3)
-      _patch0->F(0.,1.-t,x,y,z);
-    else if (_II->edgeNumber == 4)
-      _patch0->F(0.5*t,0.,x,y,z);
-    else if (_II->edgeNumber == 5)
-      _patch0->F(0.5,0.5*t,x,y,z);
-    else if (_II->edgeNumber == 6)
-      _patch0->F(0.5*(1.-t),0.5,x,y,z);
-    else if (_II->edgeNumber == 7)
-      _patch0->F(0.,0.5*(1.-t),x,y,z);
-    else if (_II->edgeNumber == 8)
-      _patch0->F(0.5*(1.+t),0.,x,y,z);
-    else if (_II->edgeNumber == 9)
-      _patch0->F(1.,0.5*t,x,y,z);
-    else if (_II->edgeNumber == 10)
-      _patch0->F(1.-0.5*t,0.5,x,y,z);
-    else if (_II->edgeNumber == 11)
-      _patch0->F(0.5,0.5*(1.-t),x,y,z);
-    else if (_II->edgeNumber == 12)
-      _patch0->F(0.5*(1.+t),0.5,x,y,z);
-    else if (_II->edgeNumber == 13)
-      _patch0->F(1.,0.5*(1.+t),x,y,z);
-    else if (_II->edgeNumber == 14)
-      _patch0->F(1.-0.5*t,1.,x,y,z);
-    else if (_II->edgeNumber == 15)
-      _patch0->F(0.5,1.-0.5*t,x,y,z);
-    else if (_II->edgeNumber == 16)
-      _patch0->F(0.5*t,0.5,x,y,z);
-    else if (_II->edgeNumber == 17)
-      _patch0->F(0.5,0.5*(1.+t),x,y,z);
-    else if (_II->edgeNumber == 18)
-      _patch0->F(0.5*(1.-t),1.,x,y,z);
-    else if (_II->edgeNumber == 19)
-      _patch0->F(0.,1.-0.5*t,x,y,z);
-    else
-      Msg::Error("Unknown edge number : %d", _II->edgeNumber);
-  }
-  else if (_patch1) {
-    if (_II->edgeNumber == 0)
-      _patch1->F(t,0.,x,y,z);
-    else if (_II->edgeNumber == 1)
-      _patch1->F(1.,t,x,y,z);
-    else if (_II->edgeNumber == 2)
-      _patch1->F(1.-t,1.,x,y,z);
-    else if (_II->edgeNumber == 3)
-      _patch1->F(0.,1.-t,x,y,z);
-    else if (_II->edgeNumber == 4)
-      _patch1->F(0.5*t,0.,x,y,z);
-    else if (_II->edgeNumber == 5)
-      _patch1->F(0.5,0.5*t,x,y,z);
-    else if (_II->edgeNumber == 6)
-      _patch1->F(0.5*(1.-t),0.5,x,y,z);
-    else if (_II->edgeNumber == 7)
-      _patch1->F(0.,0.5*(1.-t),x,y,z);
-    else if (_II->edgeNumber == 8)
-      _patch1->F(0.5*(1.+t),0.,x,y,z);
-    else if (_II->edgeNumber == 9)
-      _patch1->F(1.,0.5*t,x,y,z);
-    else if (_II->edgeNumber == 10)
-      _patch1->F(1.-0.5*t,0.5,x,y,z);
-    else if (_II->edgeNumber == 11)
-      _patch1->F(0.5,0.5*(1.-t),x,y,z);
-    else if (_II->edgeNumber == 12)
-      _patch1->F(0.5*(1.+t),0.5,x,y,z);
-    else if (_II->edgeNumber == 13)
-      _patch1->F(1.,0.5*(1.+t),x,y,z);
-    else if (_II->edgeNumber == 14)
-      _patch1->F(1.-0.5*t,1.,x,y,z);
-    else if (_II->edgeNumber == 15)
-      _patch1->F(0.5,1.-0.5*t,x,y,z);
-    else if (_II->edgeNumber == 16)
-      _patch1->F(0.5*t,0.5,x,y,z);
-    else if (_II->edgeNumber == 17)
-      _patch1->F(0.5,0.5*(1.+t),x,y,z);
-    else if (_II->edgeNumber == 18)
-      _patch1->F(0.5*(1.-t),1.,x,y,z);
-    else if (_II->edgeNumber == 19)
-      _patch1->F(0.,1.-0.5*t,x,y,z);
-    else
-      Msg::Error("Unknown edge number : %d", _II->edgeNumber);
-  }
-  */
-  else
-    Msg::Error("No patches to find the intersection curve.");
-}
-
-bool Curve::Inverse(double x,double y,double z,double &t)
-{
-  if ((_patch0) && (_patch1)) {  
-    double u0, v0;
-    _patch0->Inverse(x,y,z,u0,v0);
-    if (_along) {
-      if ((std::abs(_SP0[1]-_EP0[1])<1.e-12) && 
-	  (_patch0->_PI->periodic[1] == 1)) {
-	t = (v0 - _SP0[1] > 0 ? v0 - _SP0[1] : 1. + v0 - _SP0[1]);
-      }
-      else
-	t = (v0 - _SP0[1]) / (_EP0[1] - _SP0[1]);
-    }
-    else {
-      if ((std::abs(_SP0[0]-_EP0[0])<1.e-12) && 
-	  (_patch0->_PI->periodic[0] == 1)) {
-	t = (u0 - _SP0[0] > 0 ? u0 - _SP0[0] : 1. + u0 - _SP0[0]);
-      }
-      else
-	t = (u0 - _SP0[0]) / (_EP0[0] - _SP0[0]);
-    }
-  }
-  else if (_patch0) {
-    double u,v;
-    _patch0->Inverse(x,y,z,u,v);
-    if (_II->edgeInfo.edgeType == 0) {
-      double uS = _II->edgeInfo.startValue;
-      double uE = _II->edgeInfo.endValue;
-      if (_II->edgeInfo.acrossDiscontinuity) {
-	if (uS < uE)
-	  uE -= 1.;
-	else
-	  uE += 1.;
-	t = (u - uS) / (uE - uS);
-      }
-      else {
-	t = (u - uS) / (uE - uS);
-      }
-    }
-    else if (_II->edgeInfo.edgeType == 1) {
-      double vS = _II->edgeInfo.startValue;
-      double vE = _II->edgeInfo.endValue;
-      if (_II->edgeInfo.acrossDiscontinuity) {
-	if (vS < vE)
-	  vE -= 1.;
-	else
-	  vE += 1.;
-	t = (v - vS) / (vE - vS);
-      }
-      else {
-	t = (v - vS) / (vE - vS);
-      }
-    }
-    else
-      Msg::Error("Unknown edge type : %d", _II->edgeInfo.edgeType); 
-  }
-  else if (_patch1) {
-    double u,v;
-    _patch1->Inverse(x,y,z,u,v);
-    if (_II->edgeInfo.edgeType == 0) {
-      double uS = _II->edgeInfo.startValue;
-      double uE = _II->edgeInfo.endValue;
-      if (_II->edgeInfo.acrossDiscontinuity) {
-	if (uS < uE)
-	  uE -= 1.;
-	else
-	  uE += 1.;
-	t = (u - uS) / (uE - uS);
-      }
-      else {
-	t = (u - uS) / (uE - uS);
-      }
-    }
-    else if (_II->edgeInfo.edgeType == 1) {
-      double vS = _II->edgeInfo.startValue;
-      double vE = _II->edgeInfo.endValue;
-      if (_II->edgeInfo.acrossDiscontinuity) {
-	if (vS < vE)
-	  vE -= 1.;
-	else
-	  vE += 1.;
-	t = (v - vS) / (vE - vS);
-      }
-      else {
-	t = (v - vS) / (vE - vS);
-      }
-    }
-    else
-      Msg::Error("Unknown edge type : %d", _II->edgeInfo.edgeType); 
-  }
-  /*
-  else if (_patch0) {
-     double u,v;
-    _patch0->Inverse(x,y,z,u,v);
-    if (_II->edgeNumber == 0)
-      t = u;
-    else if (_II->edgeNumber == 1)
-      t = v;
-    else if (_II->edgeNumber == 2)
-      t = 1. - u;
-    else if (_II->edgeNumber == 3)
-      t = 1. - v;
-    else if (_II->edgeNumber == 4)
-      t = 2.*u;
-    else if (_II->edgeNumber == 5)
-      t = 2.*v;
-    else if (_II->edgeNumber == 6)
-      t = 1. - 2.*u;
-    else if (_II->edgeNumber == 7)
-      t = 1. - 2.*v;
-    else if (_II->edgeNumber == 8)
-      t = 2.*u - 1.;
-    else if (_II->edgeNumber == 9)
-      t = 2.*v;
-    else if (_II->edgeNumber == 10)
-      t = 2.*(1. - u);
-    else if (_II->edgeNumber == 11)
-      t = 1. - 2.*v;
-    else if (_II->edgeNumber == 12)
-      t = 2.*u - 1.;
-    else if (_II->edgeNumber == 13)
-      t = 2.*v - 1.;
-    else if (_II->edgeNumber == 14)
-      t = 2.*(1.- u);
-    else if (_II->edgeNumber == 15)
-      t = 2.*(1. - v);
-    else if (_II->edgeNumber == 16)
-      t = 2.*u;
-    else if (_II->edgeNumber == 17)
-      t = 2.*v - 1.;
-    else if (_II->edgeNumber == 18)
-      t = 1. - 2.*u;
-    else if (_II->edgeNumber == 19)
-      t = 2.*(1. - v);
-    else
-      Msg::Error("Unknown edge number : %d", _II->edgeNumber);
-  }
-  else if (_patch1) {
-    double u,v;
-    _patch1->Inverse(x,y,z,u,v);
-    if (_II->edgeNumber == 0)
-      t = u;
-    else if (_II->edgeNumber == 1)
-      t = v;
-    else if (_II->edgeNumber == 2)
-      t = 1. - u;
-    else if (_II->edgeNumber == 3)
-      t = 1. - v;
-    else if (_II->edgeNumber == 4)
-      t = 2.*u;
-    else if (_II->edgeNumber == 5)
-      t = 2.*v;
-    else if (_II->edgeNumber == 6)
-      t = 1. - 2.*u;
-    else if (_II->edgeNumber == 7)
-      t = 1. - 2.*v;
-    else if (_II->edgeNumber == 8)
-      t = 2.*u - 1.;
-    else if (_II->edgeNumber == 9)
-      t = 2.*v;
-    else if (_II->edgeNumber == 10)
-      t = 2.*(1. - u);
-    else if (_II->edgeNumber == 11)
-      t = 1. - 2.*v;
-    else if (_II->edgeNumber == 12)
-      t = 2.*u - 1.;
-    else if (_II->edgeNumber == 13)
-      t = 2.*v - 1.;
-    else if (_II->edgeNumber == 14)
-      t = 2.*(1.- u);
-    else if (_II->edgeNumber == 15)
-      t = 2.*(1. - v);
-    else if (_II->edgeNumber == 16)
-      t = 2.*u;
-    else if (_II->edgeNumber == 17)
-      t = 2.*v - 1.;
-    else if (_II->edgeNumber == 18)
-      t = 1. - 2.*u;
-    else if (_II->edgeNumber == 19)
-      t = 2.*(1. - v);
-    else
-      Msg::Error("Unknown edge number : %d", _II->edgeNumber);
-  }
-  */
-  else
-    Msg::Error("No patches to find the intersection curve.");
-
-  if ((t < 0.) || (t > 1.))
-    return false;
-  else
-    return true;
-}
-
-double Curve::GetPou(double t)
-{
-  return 1.;
-}
-
-int Curve::GetTag()
-{
-  return _II->tag;
-}
-
-int Curve::GetEdgeType()
-{
-  return _II->edgeInfo.edgeType;
+  _tag = -1;
 }
 
-void Curve::GetEdgeParExtent(double &start, double &end)
+Curve::Curve(int tag)
 {
-  start = _II->edgeInfo.startValue;
-  end = _II->edgeInfo.endValue;
+  _tag = tag;
 }
diff --git a/contrib/FourierModel/Curve.h b/contrib/FourierModel/Curve.h
index dfbfb4838055f0231f40f3495ed5d5ea9f78b8b9..72e55c8c65f1f4b97ca6a0db301a6c70e4025e2b 100644
--- a/contrib/FourierModel/Curve.h
+++ b/contrib/FourierModel/Curve.h
@@ -1,40 +1,23 @@
 #ifndef _CURVE_H_
 #define _CURVE_H_
 
-#include "Patch.h"
-#include "FM_Info.h"
-
 // The base class for the patches
 class Curve {
- private:
-  int _along;
-  double _h, _tol;
  protected:
-  // Patches
-  Patch* _patch0;
-  Patch* _patch1;
-  // End Points
-  double _SP0[2],_EP0[2];
-  double _SP1[2],_EP1[2];
-
+  int _tag;
  public:
   // Intersection Information
-  IntersectionInfo* _II;
-  Curve(IntersectionInfo* II, std::vector<Patch*> patches);
+  Curve();
+  Curve(int tag);
   virtual ~Curve() {}
 
-  // These are the virtual functions that must be provided by all
-  // derived patches: GetPou() returns the original smooth
-  // (non-normalized) cutoff on the patch; F() and Inverse() implement
-  // the mapping f: (t)->(x,y,z) and its inverse; and the Df*() and Dn*()
-  // functions return the derivatives of the mapping f and unit normal n 
-  // with respect to u and v
-  virtual double GetPou(double t);
-  virtual int GetTag();
-  virtual int GetEdgeType();
-  virtual void GetEdgeParExtent(double &start, double &end);
-  virtual void F(double t, double &x, double &y, double &z);
-  virtual bool Inverse(double x,double y,double z,double &t);
+  inline int GetTag() { return _tag; }
+
+  virtual double GetPou(double t) = 0;
+  virtual void F(double t, double &x, double &y, double &z) = 0;
+  virtual bool Inverse(double x,double y,double z,double &t) = 0;
+  virtual void Dfdt(double t, double &x, double &y, double &z) = 0;
+  virtual void Dfdfdtdt(double t, double &x, double &y, double &z) = 0;
 };
 
 #endif
diff --git a/contrib/FourierModel/CylindricalProjectionSurface.cpp b/contrib/FourierModel/CylindricalProjectionSurface.cpp
index 8307d7162430cadce70f9ab9e64f2c06da6dcb56..065783e053ab99ed615dc222c20f6dba5685f9d3 100755
--- a/contrib/FourierModel/CylindricalProjectionSurface.cpp
+++ b/contrib/FourierModel/CylindricalProjectionSurface.cpp
@@ -234,3 +234,14 @@ GetParameter(int i)
     return Z_;
   }
 }
+
+std::string CylindricalProjectionSurface::
+GetLabel(int i)
+{
+  switch (i) {
+  case 0:
+    return std::string("R");
+  case 1:
+    return std::string("Hight Scale");
+  }
+}
diff --git a/contrib/FourierModel/CylindricalProjectionSurface.h b/contrib/FourierModel/CylindricalProjectionSurface.h
index 01ac33a553d7ccc8d96d65c3386a2c4a5973ebed..627142e9b49aec2fc632d75fac0ed00ef4b54e2c 100755
--- a/contrib/FourierModel/CylindricalProjectionSurface.h
+++ b/contrib/FourierModel/CylindricalProjectionSurface.h
@@ -47,6 +47,8 @@ class CylindricalProjectionSurface : public ProjectionSurface {
     (int i, double x);
   virtual double GetParameter
     (int i);
+  virtual std::string GetLabel
+    (int i);
 
   // Redefinitions for CylindricalProjectionSurface
 
diff --git a/contrib/FourierModel/ExactPatch.cpp b/contrib/FourierModel/ExactPatch.cpp
index e2a0267c54378819a126d16214a632fb2320a06e..4c7085a52327f8520fe2b4e9c1c6aa2d1bb00dd2 100644
--- a/contrib/FourierModel/ExactPatch.cpp
+++ b/contrib/FourierModel/ExactPatch.cpp
@@ -1,8 +1,9 @@
 #include "Message.h"
 #include "ExactPatch.h"
 
-ExactPatch::ExactPatch(PatchInfo *PI, ProjectionSurface* ps) : Patch(PI)
+ExactPatch::ExactPatch(PatchInfo *PI, ProjectionSurface* ps) : Patch()
 {
+  _PI = PI;
   SetProjectionSurface(ps);
 
   if (!strcmp(_PI->type,"plane")) {
diff --git a/contrib/FourierModel/ExactPatch.h b/contrib/FourierModel/ExactPatch.h
index d6f94dfe32f34ebd932ae257ea96ede152f1f03f..21eba4f064b4112c8dc6f31e14e80a6ae31b22cf 100644
--- a/contrib/FourierModel/ExactPatch.h
+++ b/contrib/FourierModel/ExactPatch.h
@@ -2,6 +2,7 @@
 #define _EXACT_PATCH_H_
 
 #include "Patch.h"
+#include "FM_Info.h"
 #include "ProjectionSurface.h"
 
 // The base class for the patches
@@ -10,6 +11,8 @@ class ExactPatch : public Patch {
   ExactPatch(PatchInfo* PI, ProjectionSurface* ps);
   virtual ~ExactPatch() {}
 
+  PatchInfo* _PI;
+
   // These are the virtual functions that must be provided by all
   // derived patches: GetPou() returns the original smooth
   // (non-normalized) cutoff on the patch; F() and Inverse() implement
diff --git a/contrib/FourierModel/FCurve.cpp b/contrib/FourierModel/FCurve.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..67d87913f64412895059a968252dd86db24af056
--- /dev/null
+++ b/contrib/FourierModel/FCurve.cpp
@@ -0,0 +1,94 @@
+#include "FCurve.h"
+
+FCurve::FCurve
+(int tag, Patch* patch) 
+  : Curve(tag)
+{
+  patch_ = patch;
+
+  // default end points
+  SP_[0] = SP_[1] = 0.;
+  EP_[0] = 1.; EP_[1] = 0.;
+}
+
+FCurve::FCurve
+(int tag, Patch* patch, double SP[2], double EP[2]) 
+  : Curve(tag)
+{
+  patch_ = patch;
+
+  SP_[0] = SP[0]; SP_[1] = SP[1];
+  EP_[0] = EP[0]; EP_[1] = EP[1];
+}
+
+double FCurve::GetPou
+(double t)
+{
+  return 1.;
+}
+
+void FCurve::F
+(double t, double &x, double &y, double &z)
+{
+  double u = SP_[0] + t * (EP_[0] - SP_[0]);
+  double v = SP_[1] + t * (EP_[1] - SP_[1]);
+
+  patch_->F(u,v,x,y,z);
+}
+
+bool FCurve::Inverse
+(double x, double y, double z, double &t)
+{
+  double u, v;
+  if (patch_->Inverse(x,y,z,u,v)) {
+    if (std::abs(EP_[0] - SP_[0]) > 1.e-12)
+      t = (u - SP_[0]) / (EP_[0] - SP_[0]);
+    else if (std::abs(EP_[1] - SP_[1]) > 1.e-12)
+      t = (v - SP_[1]) / (EP_[1] - SP_[1]);
+    else
+      t = 0.;
+    return true;
+  }
+  else
+    return false;
+}
+
+void FCurve::Dfdt
+(double t, double &x, double &y, double &z)
+{
+  double u = SP_[0] + t * (EP_[0] - SP_[0]);
+  double v = SP_[1] + t * (EP_[1] - SP_[1]);
+
+  double dfdu[3], dfdv[3];
+  patch_->Dfdu(u,v,dfdu[0],dfdu[1],dfdu[2]);
+  patch_->Dfdv(u,v,dfdv[0],dfdv[1],dfdv[2]);
+
+  x = dfdu[0] * (EP_[0] - SP_[0]) + dfdv[0] * (EP_[1] - SP_[1]);
+  y = dfdu[1] * (EP_[0] - SP_[0]) + dfdv[1] * (EP_[1] - SP_[1]);
+  z = dfdu[2] * (EP_[0] - SP_[0]) + dfdv[2] * (EP_[1] - SP_[1]);
+}
+
+void FCurve::Dfdfdtdt
+(double t, double &x, double &y, double &z)
+{
+  double u = SP_[0] + t * (EP_[0] - SP_[0]);
+  double v = SP_[1] + t * (EP_[1] - SP_[1]);
+
+  double dfdfdudu[3], dfdfdudv[3], dfdfdvdv[3];
+  patch_->Dfdfdudu(u,v,dfdfdudu[0],dfdfdudu[1],dfdfdudu[2]);
+  patch_->Dfdfdudv(u,v,dfdfdudv[0],dfdfdudv[1],dfdfdudv[2]);
+  patch_->Dfdfdvdv(u,v,dfdfdvdv[0],dfdfdvdv[1],dfdfdvdv[2]);  
+
+  x = 
+    dfdfdudu[0] * (EP_[0] - SP_[0]) * (EP_[0] - SP_[0]) +
+    dfdfdudv[0] * (EP_[0] - SP_[0]) * (EP_[1] - SP_[1]) * 2. +
+    dfdfdvdv[0] * (EP_[1] - SP_[1]) * (EP_[1] - SP_[1]);
+  y = 
+    dfdfdudu[1] * (EP_[0] - SP_[0]) * (EP_[0] - SP_[0]) +
+    dfdfdudv[1] * (EP_[0] - SP_[0]) * (EP_[1] - SP_[1]) * 2. +
+    dfdfdvdv[1] * (EP_[1] - SP_[1]) * (EP_[1] - SP_[1]);
+  z = 
+    dfdfdudu[2] * (EP_[0] - SP_[0]) * (EP_[0] - SP_[0]) +
+    dfdfdudv[2] * (EP_[0] - SP_[0]) * (EP_[1] - SP_[1]) * 2. +
+    dfdfdvdv[2] * (EP_[1] - SP_[1]) * (EP_[1] - SP_[1]);
+}
diff --git a/contrib/FourierModel/FCurve.h b/contrib/FourierModel/FCurve.h
new file mode 100644
index 0000000000000000000000000000000000000000..18267398aeaa8ac70e50a79c1677ee84fb914324
--- /dev/null
+++ b/contrib/FourierModel/FCurve.h
@@ -0,0 +1,60 @@
+#ifndef _F_CURVE_H_
+#define _F_CURVE_H_
+
+#include "Curve.h"
+#include "Patch.h"
+
+// The base class for the patches
+class FCurve : public Curve {
+ protected:
+  // Underlying Patch
+  Patch* patch_;
+  // End Points
+  double SP_[2],EP_[2];
+
+ public:
+  FCurve(int tag, Patch* patch);
+  FCurve(int tag, Patch* patch, double SP[2], double EP[2]);  
+  virtual ~FCurve() {}
+
+  void SetStartPoint
+    (double u, double v) 
+  { 
+    SP_[0] = u; SP_[1] = v; 
+  }
+  void GetStartPoint
+    (double &u, double &v) 
+  { 
+    u = SP_[0]; v = SP_[1]; 
+  }
+  void SetEndPoint
+    (double u, double v) 
+  { 
+    EP_[0] = u; EP_[1] = v; 
+  }
+  void GetEndPoint
+    (double &u, double &v) 
+  { 
+    u = EP_[0]; v = EP_[1]; 
+  }
+
+  // These are the virtual functions that must be provided by all
+  // derived patches: GetPou() returns the original smooth
+  // (non-normalized) cutoff on the patch; F() and Inverse() implement
+  // the mapping f: (t)->(x,y,z) and its inverse; and the Df*() and Dn*()
+  // functions return the derivatives of the mapping f and unit normal n 
+  // with respect to u and v
+
+  virtual double GetPou
+    (double t);
+  virtual void F
+    (double t, double &x, double &y, double &z);
+  virtual bool Inverse
+    (double x,double y,double z,double &t);
+  virtual void Dfdt
+    (double t, double &x, double &y, double &z);
+  virtual void Dfdfdtdt
+    (double t, double &x, double &y, double &z);
+};
+
+#endif
diff --git a/contrib/FourierModel/FM_Edge.cpp b/contrib/FourierModel/FM_Edge.cpp
index 138560341bc8a461cbee2657a804ca37ae7087ce..f5641b8c47655981cbc81e17c90aeef0a1c85edd 100644
--- a/contrib/FourierModel/FM_Edge.cpp
+++ b/contrib/FourierModel/FM_Edge.cpp
@@ -1,42 +1,11 @@
+#include <cmath>
 #include "FM_Edge.h"
 #include "Message.h"
 
-bool FM_Edge::GetCurveExtent(double &start, double &end)
-{
-  bool result = false;
-  if (_curve) {
-    printf("\t\tedgeType : %d\n",_curve->GetEdgeType());
-    if (_curve->GetEdgeType() >= 0) {
-      _curve->GetEdgeParExtent(start,end);
-      result = true;
-    }
-  }
-
-  return result;
-}
-
 void FM_Edge::F(double t, double &x, double &y, double &z)
 {
   if (_curve) {
-    double tRescaled;
-    if (_curve->GetEdgeType() < 0) {
-      double tStart, tEnd;
-      _curve->Inverse(_SP->GetX(),_SP->GetY(),_SP->GetZ(),tStart);
-      _curve->Inverse(_EP->GetX(),_EP->GetY(),_EP->GetZ(),tEnd);
-      
-      if (std::abs(tEnd - tStart) < 1.e-12) {
-	tRescaled = tStart + t * (1. + tEnd - tStart);
-	tRescaled -= floor(tRescaled);
-      }
-      else
-	tRescaled = tStart + t * (tEnd - tStart);
-    }
-    else
-      tRescaled = t;
-    _curve->F(tRescaled, x, y, z);
-    //Msg::Info("%g %g %g",_SP->GetX(),_SP->GetY(),_SP->GetZ());
-    //Msg::Info("%g %g %g",_EP->GetX(),_EP->GetY(),_EP->GetZ());
-    //Msg::Info("t : %g %g %g %g",t,tRescaled, tStart, tEnd);
+    _curve->F(t,x,y,z);
   }
   else {
     x = _SP->GetX() + t * (_EP->GetX() - _SP->GetX());
@@ -48,19 +17,7 @@ void FM_Edge::F(double t, double &x, double &y, double &z)
 bool FM_Edge::Inverse(double x,double y,double z,double &t)
 {
   if (_curve) {
-    double tStart, tEnd;
-    _curve->Inverse(_SP->GetX(),_SP->GetY(),_SP->GetZ(),tStart);
-    _curve->Inverse(_EP->GetX(),_EP->GetY(),_EP->GetZ(),tEnd);
-    
-    double tCurve;
-    _curve->Inverse(x, y, z, tCurve);
-
-    if (std::abs(tEnd - tStart) < 1.e-12) {
-      t = (tCurve - tStart) / (1. + tEnd - tStart);
-      t -= floor(t);
-    }
-    else
-      t = (tCurve - tStart) / (tEnd - tStart);
+    _curve->Inverse(x,y,z,t);
   }
   else {
     if (_EP->GetX() - _SP->GetX())
diff --git a/contrib/FourierModel/FM_Face.cpp b/contrib/FourierModel/FM_Face.cpp
index 0443f9e0739da0dcc6b782c8341c9ba97cabcfbd..034ab80fcf81de2dd0fc6edb1d6843c73157db2c 100644
--- a/contrib/FourierModel/FM_Face.cpp
+++ b/contrib/FourierModel/FM_Face.cpp
@@ -113,32 +113,27 @@ bool FM_Face::Inverse(double x,double y,double z,double &u, double &v)
 
 void FM_Face::Dfdu(double u, double v, double &x, double &y, double &z)
 {
-  Msg::Info("Not implemented yet.");
-  x = y = z = 0.;
+  _patch->Dfdu(u,v,x,y,z);
 }
 
 void FM_Face::Dfdv(double u, double v, double &x, double &y, double &z)
 {
-  Msg::Info("Not implemented yet.");
-  x = y = z = 0.;
+  _patch->Dfdv(u,v,x,y,z);
 }
 
 void FM_Face::Dfdfdudu(double u, double v, double &x, double &y, double &z)
 {
-  Msg::Info("Not implemented yet.");
-  x = y = z = 0.;
+  _patch->Dfdfdudu(u,v,x,y,z);
 }
 
 void FM_Face::Dfdfdudv(double u, double v, double &x, double &y, double &z)
 {
-  Msg::Info("Not implemented yet.");
-  x = y = z = 0.;
+  _patch->Dfdfdudv(u,v,x,y,z);
 }
 
 void FM_Face::Dfdfdvdv(double u, double v, double &x, double &y, double &z)
 {
-  Msg::Info("Not implemented yet.");
-  x = y = z = 0.;
+  _patch->Dfdfdvdv(u,v,x,y,z);
 }
 
 void FM_Face::GetNormal(double u, double v, double &x, double &y, double &z)
diff --git a/contrib/FourierModel/FM_Face.h b/contrib/FourierModel/FM_Face.h
index 813ab7e9a3c25edc62dd8a1b798aff542a740e40..213182fc11b2d9573bc8a06bbb2a0aeef3605eb1 100644
--- a/contrib/FourierModel/FM_Face.h
+++ b/contrib/FourierModel/FM_Face.h
@@ -1,6 +1,7 @@
 #ifndef _FM_FACE_H_
 #define _FM_FACE_H_
 
+#include <vector>
 #include "Patch.h"
 #include "FM_Edge.h"
 
@@ -21,7 +22,11 @@ class FM_Face {
 
   inline void SetTag(int tag) { _tag = tag; }
   inline int GetTag() { return _tag; }
-  inline bool IsPeriodic(int dim) { return _patch->_PI->periodic[dim]; }
+  inline bool IsPeriodic(int dim) 
+    { 
+      if (dim) return _patch->IsVPeriodic();
+      else return _patch->IsUPeriodic();
+    }
   inline void AddEdge(FM_Edge* edge) { _edge.push_back(edge); }
   inline int GetNumEdges() { return _edge.size(); }
   inline FM_Edge* GetEdge(int i) { return _edge[i]; }
diff --git a/contrib/FourierModel/FM_Reader.cpp b/contrib/FourierModel/FM_Reader.cpp
index ab87d1a0d46467f175f9325024376bf6c5a3f44c..15f8587d6dcd67bae1e0f4f6db6a00ec593d22f5 100644
--- a/contrib/FourierModel/FM_Reader.cpp
+++ b/contrib/FourierModel/FM_Reader.cpp
@@ -103,7 +103,7 @@ FM_Reader::FM_Reader(const char* fn)
     }
   }
 
-  _intersection.resize(_nIntersections,0);
+  _curve.resize(_nIntersections,0);
 
   for (int i=0;i<_nPatches;i++) {
     _patchList[i]->tag = i;
@@ -123,7 +123,7 @@ FM_Reader::FM_Reader(const char* fn)
 
   for (int i=0;i<_nIntersections;i++) {
     IntersectionInfo* II = _intersectionList[i];
-    _intersection[II->tag] = new Curve(II,_patch);
+    _curve[II->tag] = new IntersectionCurve(II,_patch);
   }
 
   InputFile >> _nVertices;
@@ -140,7 +140,7 @@ FM_Reader::FM_Reader(const char* fn)
     if (edgeTag < 0)
       _edge.push_back(new FM_Edge(i,0,_vertex[svTag],_vertex[evTag]));
     else
-      _edge.push_back(new FM_Edge(i,GetIntersection(edgeTag),
+      _edge.push_back(new FM_Edge(i,GetCurve(edgeTag),
 				  _vertex[svTag],_vertex[evTag]));
   }
 
@@ -165,12 +165,12 @@ Patch* FM_Reader::GetPatch(int tag)
       return _patch[i];
 }
 
-Curve* FM_Reader::GetIntersection(int tag)
+Curve* FM_Reader::GetCurve(int tag)
 {
   Curve* curve = 0;
-  for (int i=0;i<_intersection.size();i++) {
-    if (_intersection[i]->GetTag() == tag) {
-      curve = _intersection[i];
+  for (int i=0;i<_curve.size();i++) {
+    if (_curve[i]->GetTag() == tag) {
+      curve = _curve[i];
       break;
     }
   }
diff --git a/contrib/FourierModel/FM_Reader.h b/contrib/FourierModel/FM_Reader.h
index c00c4c1527526314ac267f251176edcec7b6bafe..9b229767bba1c5fe4da111000863eb60ad3544b7 100644
--- a/contrib/FourierModel/FM_Reader.h
+++ b/contrib/FourierModel/FM_Reader.h
@@ -4,7 +4,9 @@
 #include <vector>
 #include <string>
 #include <complex>
+#include "Point.h"
 #include "Curve.h"
+#include "IntersectionCurve.h"
 #include "ExactPatch.h"
 #include "ContinuationPatch.h"
 #include "CylindricalProjectionSurface.h"
@@ -29,7 +31,7 @@ class FM_Reader {
   std::vector<FM_Edge*> _edge;
   std::vector<FM_Face*> _face;
   std::vector<Patch*> _patch;
-  std::vector<Curve*> _intersection;
+  std::vector<Curve*> _curve;
   std::vector<ProjectionSurface*> _ps;
 
   BlendOperator* _blendOperator;
@@ -83,7 +85,7 @@ class FM_Reader {
   BlendedPatch* GetBlendedPatch
     (int tag);
 
-  Curve* GetIntersection
+  Curve* GetCurve
     (int tag);
 
   ProjectionSurface* GetProjectionSurface
diff --git a/contrib/FourierModel/FPatch.cpp b/contrib/FourierModel/FPatch.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..df3c873346a520e69d8b18cf9efe7053c2b0a6a6
--- /dev/null
+++ b/contrib/FourierModel/FPatch.cpp
@@ -0,0 +1,850 @@
+#include <math.h>
+#include "Message.h"
+#include "FPatch.h"
+
+extern "C" {
+  void zgelss_(int &,int &,int &,std::complex<double> *,int &,
+               std::complex<double> *,int &,double *,double &,int &,
+               std::complex<double> *,int &,double *,int &);
+}
+
+FPatch::FPatch(int tag, ProjectionSurface* ps, 
+	       std::vector<double> &u, std::vector<double> &v,
+	       std::vector< std::complex<double> > &data, int derivative) 
+  : _coeffOriginalData(0),_coeffData(0),_coeffDerivU(0),
+    _coeffDerivV(0),_coeffDerivUU(0),_coeffDerivVV(0),_coeffDerivUV(0)
+{
+  _ps = ps;
+
+  _tag = tag;
+  _derivative = derivative;
+
+  _uModes = 10;
+  _vModes = 8;
+
+  _uM = 16+1;
+  _vM = 16;
+
+  _hardEdge[0] = false;
+  _hardEdge[1] = false;
+  _hardEdge[2] = false;
+  _hardEdge[3] = false;
+
+  if (_ps->IsUPeriodic())
+    _periodicityU = 1;
+  if (_ps->IsVPeriodic())
+    _periodicityV = 1;
+
+  if (_ps->IsUPeriodic())
+    _periodU = 1.;
+  else
+    _periodU = 2.;
+
+  if (_ps->IsVPeriodic())
+    _periodV = 1.;
+  else
+    _periodV = 2.;
+
+  if ((_uModes % 2) == 0) {
+    _uModesLower = - _uModes/2;
+    _uModesUpper = _uModes/2 - 1;
+  }
+  else {
+    _uModesLower = - (_uModes - 1)/2;
+    _uModesUpper = (_uModes - 1)/2;
+  }
+  if ((_vModes % 2) == 0) {
+    _vModesLower = - _vModes/2;
+    _vModesUpper = _vModes/2 - 1;
+  }
+  else {
+    _vModesLower = - (_vModes - 1)/2;
+    _vModesUpper = (_vModes - 1)/2;
+  }
+
+  if (_ps->IsUPeriodic()) {
+    if ((_uM % 2) == 0) {
+      _uMLower = - _uM/2;
+      _uMUpper = _uM/2 - 1;
+    }
+    else {
+      _uMLower = - (_uM - 1)/2;
+      _uMUpper = (_uM - 1)/2;
+    }
+  }
+  else {
+    _uMLower = 0;
+    _uMUpper = _uM;
+  }
+  if (_ps->IsVPeriodic()) {
+    if ((_vM % 2) == 0) {
+      _vMLower = - _vM/2;
+      _vMUpper = _vM/2 - 1;
+    }
+    else {
+      _vMLower = - (_vM - 1)/2;
+      _vMUpper = (_vM - 1)/2;
+    }
+  }
+  else {
+    _vMLower = 0;
+    _vMUpper = _vM;    
+  }
+
+
+  // Initialize Data
+
+  _coeffOriginalData = new std::complex<double>*[_uModes];
+  for(int j = 0; j < _uModes; j++){
+    _coeffOriginalData[j] = new std::complex<double>[_vModes];
+
+    //for(int k = 0; k < _vModes; k++)
+    //_coeffOriginalData[j][k] = _PI->coeff[j][k];
+  }
+
+  _nData = data.size();
+  _nModes = _uModes * _vModes;
+
+  std::complex<double> *LSRhs = new std::complex<double> [_nData];
+  std::complex<double> *LSMatrix = new std::complex<double> [_nModes * _nData];
+
+  for (int i=0;i<_nData;i++)
+    LSRhs[i] = data[i];
+
+  for (int j=0;j<_uModes;j++)
+    for (int k=0;k<_vModes;k++)
+      for (int i=0;i<_nData;i++)
+        LSMatrix[i+_nData*(k+_vModes*j)] = std::complex<double>
+          (cos(2*M_PI * ((j + _uModesLower) * u[i] / _periodU +
+                         (k + _vModesLower) * v[i] / _periodV)),
+           sin(2*M_PI * ((j + _uModesLower) * u[i] / _periodU +
+                         (k + _vModesLower) * v[i] / _periodV)));
+
+  // some parameters for the lease square solvers "zgelss"
+  int info, rank;
+  int lwork = 66*_nData, one = 1;
+  double rcond = -1.;
+  double *s = new double [_nModes];
+  double *rwork = new double [5*_nModes];
+  std::complex<double> *work = new std::complex<double> [lwork];
+
+  zgelss_(_nData,_nModes,one,&LSMatrix[0],_nData,&LSRhs[0],_nData,&s[0],rcond,
+          rank,&work[0],lwork,&rwork[0],info);
+
+  for(int j = 0; j < _uModes; j++)
+    for(int k = 0; k < _vModes; k++)
+      _coeffOriginalData[j][k] = LSRhs[k+_vModes*j];
+
+  // deleting lease square arrays
+  delete[] s, rwork, work;
+  delete[] LSMatrix, LSRhs;
+
+
+  // Initialize interpolation variables
+  _tmpOrigCoeff = std::vector< std::complex<double> >(_vModes);
+  _tmpOrigInterp = std::vector< std::complex<double> >(_uModes);
+
+  if (_derivative)
+    _ReprocessSeriesCoeff();
+
+  // Check if we need to interpolate the derivative(s)
+  if(_derivative){
+    // Initialize _fineDeriv and _fineDeriv2 to zero
+    if(_derivative & 1){
+      _coeffDerivU = new std::complex<double>*[_uM];
+      _coeffDerivV = new std::complex<double>*[_uM];
+      for(int j = 0; j < _uM; j++){
+        _coeffDerivU[j] = new std::complex<double>[_vM];
+        _coeffDerivV[j] = new std::complex<double>[_vM];
+        for(int k = 0; k < _vM; k++){
+          _coeffDerivU[j][k] = 0.;
+          _coeffDerivV[j][k] = 0.;
+        }
+      }
+    }
+
+    if(_derivative & 2){
+      _coeffDerivUU = new std::complex<double>*[_uM];
+      _coeffDerivVV = new std::complex<double>*[_uM];
+      _coeffDerivUV = new std::complex<double>*[_uM];
+      for(int j = 0; j < _uM; j++){
+        _coeffDerivUU[j] = new std::complex<double>[_vM];
+        _coeffDerivVV[j] = new std::complex<double>[_vM];
+        _coeffDerivUV[j] = new std::complex<double>[_vM];
+        for(int k = 0; k < _vM; k++){
+          _coeffDerivUU[j][k] = 0.;
+          _coeffDerivVV[j][k] = 0.;
+          _coeffDerivUV[j][k] = 0.;
+        }
+      }
+    }
+
+    // Copy the Fourier coefficients into _coeffDeriv and _coeffDeriv2
+    std::complex<double> I(0., 1.);
+    for(int j = _uM - 1; j >= 0; j--){
+      for(int k = _vM - 1; k >= 0; k--){
+        if(_derivative & 1){
+	  if (_ps->IsUPeriodic()) {
+	    int J = j+_uMLower;
+	    _coeffDerivU[j][k] = (2 * M_PI * J * I / _periodU) *
+	      _coeffData[j][k];
+	  }
+	  else {
+	    if (j == _uM - 1)
+	      _coeffDerivU[j][k] = 0.;
+	    else if (j == _uM - 2)
+	      _coeffDerivU[j][k] = 2. * (double)(j + 1) * _coeffData[j + 1][k];
+	    else
+	      _coeffDerivU[j][k] = _coeffDerivU[j + 2][k] +
+		2. * (double)(j + 1) * _coeffData[j + 1][k];
+	    //if (j != 0)
+	    //_coeffDerivU[j][k] *= 2.;
+	  }
+	  if (_ps->IsVPeriodic()) {
+	    int K = k+_vMLower;
+	    _coeffDerivV[j][k] = (2 * M_PI * K * I / _periodV) *
+	      _coeffData[j][k];
+	  }
+	  else {
+	    if (k == _vM - 1)
+	      _coeffDerivV[j][k] = 0.;
+	    else if (k == _vM - 2)
+	      _coeffDerivV[j][k] = 2. * (double)(k + 1) * _coeffData[j][k + 1];
+	    else
+	      _coeffDerivV[j][k] = _coeffDerivV[j][k + 2] +
+		2. * (double)(k + 1) * _coeffData[j][k + 1];
+	    //if (k != 0)
+	    //_coeffDerivV[j][k] *= 2.;
+	  }
+        }
+      }
+    }
+    for (int j = 0; j < _uM; j++) {
+      for (int k = 0; k < _vM; k++) {
+	if (_derivative & 1) {
+	  if (!_ps->IsUPeriodic()) {
+	    if (j != 0) {
+	      _coeffDerivU[j][k] *= 2.;
+	    }
+	  }
+	  if (!_ps->IsVPeriodic()) {
+	    if (k != 0) {
+	      _coeffDerivV[j][k] *= 2.;
+	    }
+	  }
+	}
+      }
+    }
+
+    for(int j = _uM - 1; j >= 0; j--) {
+      for(int k = _vM - 1; k >= 0; k--) {
+	if(_derivative & 2) {
+	  if (_ps->IsUPeriodic()) {
+	    int J = j+_uMLower;
+	    _coeffDerivUU[j][k] = (2 * M_PI * J * I / _periodU) * 
+	      _coeffDerivU[j][k];
+	  }
+	  else {
+	    if (j == _uM - 1)
+	      _coeffDerivUU[j][k] = 0.;
+	    else if (j == _uM - 2)
+	      _coeffDerivUU[j][k] = 2. * (double)(j + 1) * 
+		_coeffDerivU[j + 1][k];
+	    else
+	      _coeffDerivUU[j][k] = _coeffDerivUU[j + 2][k] +
+		2. * (double)(j + 1) * _coeffDerivU[j + 1][k];
+	    //if (j != 0)
+	    //_coeffDerivUU[j][k] *= 2.;
+	  }
+	  if (_ps->IsVPeriodic()) {
+	    int K = k+_vMLower;
+	    _coeffDerivVV[j][k] = (2 * M_PI * K * I / _periodV) * 
+	      _coeffDerivV[j][k];
+	    _coeffDerivUV[j][k] = (2 * M_PI * K * I / _periodV) * 
+	      _coeffDerivU[j][k];
+	  }
+	  else {
+	    if (k == _vM - 1) {
+	      _coeffDerivVV[j][k] = 0.;
+	      _coeffDerivUV[j][k] = 0.;
+	    }
+	    else if (k == _vM - 2) {
+	      _coeffDerivVV[j][k] = 2. * (double)(k + 1) * 
+		_coeffDerivV[j][k + 1];
+	      _coeffDerivUV[j][k] = 2. * (double)(k + 1) * 
+		_coeffDerivU[j][k + 1];
+	    }
+	    else {
+	      _coeffDerivVV[j][k] = _coeffDerivVV[j][k + 2] +
+		2. * (double)(k + 1) * _coeffDerivV[j][k + 1];
+	      _coeffDerivUV[j][k] = _coeffDerivUV[j][k + 2] +
+		2. * (double)(k + 1) * _coeffDerivU[j][k + 1];
+	    }
+	  }
+	}
+      }
+    }
+    for (int j = 0; j < _uM; j++) {
+      for (int k = 0; k < _vM; k++) {
+	if (_derivative & 2) {
+	  if (!_ps->IsUPeriodic() && _ps->IsVPeriodic()) {
+	    if (j != 0) {
+	      _coeffDerivUU[j][k] *= 2.;
+	    }
+	  }
+	  if (_ps->IsUPeriodic() && !_ps->IsVPeriodic()) {
+	    if (k != 0) {
+	      _coeffDerivVV[j][k] *= 2.;
+	      _coeffDerivUV[j][k] *= 2.;
+	    }
+	  }
+	}
+      }
+    }
+  }
+
+  // Initialize interpolation variables
+  _tmpCoeff = std::vector< std::complex<double> >(_vM);
+  _tmpInterp = std::vector< std::complex<double> >(_uM);
+
+}
+
+FPatch::~FPatch()
+{
+  for(int j = 0; j < _uModes; j++)
+    delete [] _coeffOriginalData[j];
+  delete [] _coeffOriginalData;
+
+  for(int j = 0; j < _uM; j++)
+    delete [] _coeffData[j];
+  delete [] _coeffData;
+
+  if(_coeffDerivU){
+    for(int j = 0; j < _uM; j++)
+      delete [] _coeffDerivU[j];
+    delete [] _coeffDerivU;
+  }
+  if(_coeffDerivV){
+    for(int j = 0; j < _uM; j++)
+      delete [] _coeffDerivV[j];
+    delete [] _coeffDerivV;
+  }
+  if(_coeffDerivUU){
+    for(int j = 0; j < _uM; j++)
+      delete [] _coeffDerivUU[j];
+    delete [] _coeffDerivUU;
+  }
+  if(_coeffDerivVV){
+    for(int j = 0; j < _uM; j++)
+      delete [] _coeffDerivVV[j];
+    delete [] _coeffDerivVV;
+  }
+  if(_coeffDerivUV){
+    for(int j = 0; j < _uM; j++)
+      delete [] _coeffDerivUV[j];
+    delete [] _coeffDerivUV;
+  }
+}
+
+int FPatch::_forwardSize = 0;
+int FPatch::_backwardSize = 0;
+fftw_plan FPatch::_forwardPlan;
+fftw_plan FPatch::_backwardPlan;
+fftw_complex *FPatch::_forwardData = 0;
+fftw_complex *FPatch::_backwardData = 0;
+
+void FPatch::_SetForwardPlan(int n)
+{
+  if(n != _forwardSize){
+    if(_forwardSize){
+      fftw_destroy_plan(_forwardPlan);
+      fftw_free(_forwardData);
+    }
+    _forwardSize = n;
+    _forwardData = 
+      (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * _forwardSize);
+    _forwardPlan = 
+      fftw_plan_dft_1d(_forwardSize, _forwardData, _forwardData,
+		       FFTW_FORWARD, FFTW_ESTIMATE);
+  }
+}
+
+void FPatch::_SetBackwardPlan(int n)
+{
+  if(n != _backwardSize){
+    if(_backwardSize){
+      fftw_destroy_plan(_backwardPlan);
+      fftw_free(_backwardData);
+    }
+    _backwardSize = n;
+    _backwardData = 
+      (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * _backwardSize);
+    _backwardPlan = 
+      fftw_plan_dft_1d(_backwardSize, _backwardData, _backwardData,
+		       FFTW_BACKWARD, FFTW_ESTIMATE);
+  }
+}
+
+void FPatch::_ForwardFft(int n, std::complex<double> *fftData)
+{
+  // Initialize fftw plan and array (ignoring the last element of
+  // fftData, which should just be the periodic extension)
+  _SetForwardPlan(n - 1);
+  for(int i = 0; i < n - 1; i++){
+    _forwardData[i][0] = fftData[i].real();
+    _forwardData[i][1] = fftData[i].imag();
+  }
+
+  // Perform forward FFT
+  fftw_execute(_forwardPlan);
+
+  // Copy data back into fftData and scale by 1/(n - 1)
+  double s = 1. / (double)(n - 1);
+  for(int i = 0; i < n - 1; i++)
+    fftData[i] = 
+      s * std::complex<double>(_forwardData[i][0], _forwardData[i][1]);
+}
+
+void FPatch::_BackwardFft(int n, std::complex<double> *fftData)
+{
+  // Initialize fftw plan and array (ignoring last element of fftData)
+  _SetBackwardPlan(n - 1);
+  for(int i = 0; i < n - 1; i++){
+    _backwardData[i][0] = fftData[i].real();
+    _backwardData[i][1] = fftData[i].imag();
+  }
+
+  // Perform backward FFT
+  fftw_execute(_backwardPlan);
+
+  // Copy data back into fftData
+  for(int i = 0; i < n - 1; i++)
+    fftData[i] = 
+      std::complex<double>(_backwardData[i][0], _backwardData[i][1]);
+
+  // Fill in last element with copy of first element
+  fftData[n - 1] = fftData[0];
+}
+
+void FPatch::_ReprocessSeriesCoeff()
+{
+  _coeffData = new std::complex<double>*[_uM];
+  for(int j = 0; j < _uM; j++) {
+    _coeffData[j] = new std::complex<double>[_vM];
+    for(int k = 0; k < _vM; k++)
+      _coeffData[j][k] = 0.;
+  }
+  if (_ps->IsUPeriodic() && _ps->IsVPeriodic()) {
+    int uShift = (_uM-_uModes)%2 == 0 ? (_uM-_uModes)/2 : (_uM-_uModes)/2 + 1;
+    int vShift = (_vM-_vModes)%2 == 0 ? (_vM-_vModes)/2 : (_vM-_vModes)/2 + 1;
+    for (int j = 0; j < _uModes; j++)
+      for (int k = 0; k < _vModes; k++)
+	_coeffData[uShift + j][vShift + k] = _coeffOriginalData[j][k];
+  }
+  else if (_ps->IsUPeriodic()) {
+    std::vector<double> u(_uM), v(_vM + 1);
+    for (int j = 0; j < _uM; j++)
+      u[j] = (double)j / (double)(_uM-1);
+    for (int j = 0; j < _vM + 1; j++)
+      v[j] = (double)j / (double)_vM;
+
+    std::complex<double> **dataU = new std::complex<double> *[_vM];
+    for (int k = 0; k < _vM; k++)
+      dataU[k] = new std::complex<double> [_uM];
+    std::complex<double> *dataV = new std::complex<double>[2*_vM + 1];
+    for (int j = 0; j < _uM - 1; j++) {
+      for (int k = 0; k < _vM + 1; k++) {
+	//dataV[k] = 1.;
+	//dataV[k] = 2. * cos(M_PI * v[k]) * cos(M_PI * v[k]) - 1.;
+	//dataV[k] = 2. * (0.5 * cos(M_PI * v[k]) + 0.5) *
+	//(0.5 * cos(M_PI * v[k]) + 0.5) - 1.;
+	//dataV[k] = 2. * (0.5 * cos(M_PI * v[k]) + 0.5) *
+	//(0.5 * cos(M_PI * v[k]) + 0.5) * cos(2 * M_PI * u[j]);
+	//dataV[k] = cos(0.5 * cos(M_PI * v[k]) + 0.5) * 
+	//cos(2 * M_PI * u[j]);
+	dataV[k] = _Interpolate(u[j],0.5 * cos(M_PI * v[k]) + 0.5);
+      }
+      for (int k = 1; k < _vM+1; k++)
+	dataV[_vM + k] = dataV[_vM -k];
+      _BackwardFft(2*_vM + 1, dataV);
+      dataU[0][j] = 0.5 * dataV[0] / (double)_vM;
+      for (int k=1; k<_vM-1; k++)
+	dataU[k][j] = dataV[k] / (double)_vM;
+      dataU[_vM-1][j] = 0.5 * dataV[_vM-1] / (double)_vM;
+    }
+    for (int k = 0; k < _vM; k++) {
+      dataU[k][_uM - 1] = dataU[k][0];
+      _ForwardFft(_uM, dataU[k]);
+    }
+    for (int j = _uMLower; j <= _uMUpper; j++) {
+      for (int k = 0; k < _vM; k++)
+	if ((j == _uMLower) || (j == _uMUpper))
+	  _coeffData[_uMUpper + j][k] = dataU[k][_uMUpper] / 2.;
+	else if ((j >= 0) && (j < _uMUpper))
+	  _coeffData[_uMUpper + j][k] = dataU[k][j];
+	else
+	  _coeffData[_uMUpper + j][k] = dataU[k][_uM + j -1];
+    }
+    for (int k = 0; k < _vM; k++)
+      delete [] dataU[k];
+    delete [] dataU;
+    delete [] dataV;
+  }
+  else {
+    std::vector<double> u(_uM + 1), v(_vM);
+    for (int j = 0; j < _uM + 1; j++)
+      u[j] = (double)j / (double)_uM;
+    for (int j = 0; j < _vM; j++)
+      v[j] = (double)j / (double)(_vM-1);
+
+    std::complex<double> **dataV = new std::complex<double> *[_uM];
+    for (int j = 0; j < _uM; j++)
+      dataV[j] = new std::complex<double> [_vM];
+    std::complex<double> *dataU = new std::complex<double>[2*_uM + 1];
+    for (int k = 0; k < _vM - 1; k++) {
+      for (int j = 0; j < _uM + 1; j++) {
+	//dataU[j] = 1.;
+	//dataU[j] = 2. * cos(M_PI * u[j]) * cos(M_PI * u[j]);
+	dataU[j] = _Interpolate(0.5 * cos(M_PI * u[j]) + 0.5,v[k]);
+      }
+      for (int j = 1; j < _uM+1; j++)
+	dataU[_uM + j] = dataU[_uM -j];
+      _BackwardFft(2*_uM + 1, dataU);
+      dataV[0][k] = 0.5 * dataU[0] / (double)_uM;
+      for (int j=1; j<_uM-1; j++)
+	dataV[j][k] = dataU[j] / (double)_uM;
+      dataV[_uM-1][k] = 0.5 * dataU[_uM-1] / (double)_uM;
+    }
+    for (int j = 0; j < _uM; j++) {
+      dataV[j][_uM - 1] = dataV[j][0];
+      _ForwardFft(_vM, dataV[j]);
+    }
+    for (int k = _vMLower; k <= _vMUpper; k++) {
+      //for (int j = 0; j < _uM - 1; j++) {
+      for (int j = 0; j < _uM; j++)
+	if ((k == _vMLower) || (k == _vMUpper))
+	  _coeffData[j][_vMUpper + k] = dataV[j][_vMUpper] / 2.;
+	else if ((j >= 0) && (j < _vMUpper))
+	  _coeffData[j][_uMUpper + k] = dataV[j][k];
+	else
+	  _coeffData[j][_uMUpper + k] = dataV[j][_uM + k -1];
+    }
+    for (int j = 0; j < _uM; j++)
+      delete [] dataV[j];
+    delete [] dataV;
+    delete [] dataU;
+  }
+}
+
+std::complex<double> FPatch::
+_PolyEval(std::vector< std::complex<double> > _coeff, std::complex<double> x)
+{
+  int _polyOrder = _coeff.size()-1;
+  std::complex<double> out = 0.;
+
+  out = x * _coeff[_polyOrder];
+  for (int i = _polyOrder - 1; i > 0; i--)
+    out = x * (out + _coeff[i]);
+  out = out + _coeff[0];
+
+  return out;
+}
+
+std::complex<double> FPatch::
+  _Interpolate(double u, double v)
+{
+  double epsilon = 1e-12;
+  if (u < 0. - epsilon || u > 1. + epsilon || 
+      v < 0. - epsilon || v > 1. + epsilon) {
+    Msg::Error("Trying to interpolate outside interval: (u,v)=(%.16g,%.16g) "
+               "not in [%g,%g]x[%g,%g]", u, v, 0., 1., 0., 1.); 
+  }
+  
+  // Interpolate to find value at (u,v)
+  for(int j = 0; j < _uModes; j++){
+    for(int k = 0; k < _vModes; k++) {
+      _tmpOrigCoeff[k] = _coeffOriginalData[j][k];
+    }
+    std::complex<double> y(cos(2 * M_PI * v / _periodV),
+                           sin(2 * M_PI * v / _periodV));
+    _tmpOrigInterp[j] = _PolyEval(_tmpOrigCoeff, y);
+    _tmpOrigInterp[j] *= std::complex<double>
+      (cos(2 * M_PI * _vModesLower * v / _periodV),
+       sin(2 * M_PI * _vModesLower * v / _periodV));
+  }
+  std::complex<double> x(cos(2 * M_PI * u / _periodU),
+                         sin(2 * M_PI * u / _periodU));
+  return _PolyEval(_tmpOrigInterp, x) * std::complex<double>
+    (cos(2 * M_PI * _uModesLower * u / _periodU),
+     sin(2 * M_PI * _uModesLower * u / _periodU));
+}
+
+std::complex<double> FPatch::
+  _Interpolate(double u, double v, int uDer, int vDer)
+{
+  //Msg::Info("%d %d %d",uDer,vDer,_derivative);
+  if (((uDer==2 || vDer==2 || (uDer==1 && vDer==1)) && !(_derivative & 2) ) ||
+      ((uDer==1 || vDer==1) && !(_derivative & 1)) ||
+      (uDer<0 || uDer>2 || vDer<0 || vDer>2) ) {
+    Msg::Error("Derivative data not available: check contructor call %d %d %d",
+               uDer,vDer,_derivative);
+    return 0.;
+  }
+
+  double epsilon = 1e-12;
+  if (u < 0. - epsilon || u > 1. + epsilon) {
+    Msg::Error("Trying to interpolate outside interval: (u,v)=(%.16g,%.16g) "
+               "not in [%g,%g]x[%g,%g]", u, v, 0., 1., 0., 1.); 
+  }
+  std::vector<double> uT(_uM,0.);
+  std::vector<double> vT(_vM,0.);
+  if (!_ps->IsUPeriodic()) {
+    for (int j = 0; j < _uM; j++)
+      if (j == 0)
+	uT[j] = 1.;
+      else if (j == 1)
+	uT[j] = 2. * u - 1;
+      else
+	uT[j] = 2. * uT[1] * uT[j-1] - uT[j-2];
+  }
+  if (!_ps->IsVPeriodic()) {
+    for (int k = 0; k < _vM; k++)
+      if (k == 0)
+	vT[k] = 1.;
+      else if (k == 1)
+	vT[k] = 2. * v - 1.;
+      else
+	vT[k] = 2. * vT[1] * vT[k-1] - vT[k-2];
+  }
+  // Interpolate to find value at (u,v)
+  for(int j = 0; j < _uM; j++){
+    _tmpInterp[j] = 0.;
+    for(int k = 0; k < _vM; k++){
+      //printf("i was here %d %d\n",j,k);
+      std::complex<double> tmp;
+      if(uDer == 0 && vDer == 0)
+	tmp = _coeffData[j][k];
+      else if(uDer == 1 && vDer == 0)
+	tmp = _coeffDerivU[j][k];
+      else if(uDer == 0 && vDer == 1)
+	tmp = _coeffDerivV[j][k];
+      else if(uDer == 2 && vDer == 0)
+	tmp = _coeffDerivUU[j][k];
+      else if(uDer == 0 && vDer == 2)
+	tmp = _coeffDerivVV[j][k];
+      else
+	tmp = _coeffDerivUV[j][k];
+      _tmpCoeff[k] = tmp;
+    }
+    //printf("i was here 00\n");
+    if (_ps->IsVPeriodic()) {
+      std::complex<double> y(cos(2 * M_PI * v / _periodV),
+			     sin(2 * M_PI * v / _periodV));
+      _tmpInterp[j] = _PolyEval(_tmpCoeff, y);
+      _tmpInterp[j] *= std::complex<double>
+	(cos(2 * M_PI * _vMLower * v / _periodV),
+	 sin(2 * M_PI * _vMLower * v / _periodV));
+    }
+    else {
+      //printf("i was here 0\n");
+      for(int k = 0; k < _vM; k++)
+	_tmpInterp[j] += _tmpCoeff[k] * vT[k];
+    }
+  }
+  //printf("i was here\n");
+  if (_ps->IsUPeriodic()) {
+    std::complex<double> x(cos(2 * M_PI * u / _periodU),
+			   sin(2 * M_PI * u / _periodU));
+    return _PolyEval(_tmpInterp, x) * std::complex<double>
+      (cos(2 * M_PI * _uMLower * u / _periodU),
+       sin(2 * M_PI * _uMLower * u / _periodU));
+  }
+  else {
+    std::complex<double> tmp = 0.;
+    for(int j = 0; j < _uM; j++)
+      tmp += _tmpInterp[j] * uT[j];
+    return tmp;
+  }
+}
+
+void FPatch::
+F(double u, double v, double &x, double &y, double &z)
+{
+  double px, py, pz, nx, ny, nz, d;
+
+  u = RescaleU(u);
+  v = RescaleV(v);
+
+  _ps->F(u,v,px,py,pz);
+  _ps->GetUnitNormal(u,v,nx,ny,nz);
+
+  d = _Interpolate(u, v).real();
+
+  x = px + d * nx;
+  y = py + d * ny;
+  z = pz + d * nz;
+}
+
+bool FPatch::
+Inverse(double x,double y,double z,double &u,double &v)
+{
+  bool result = _ps->OrthoProjectionOnSurface(x,y,z,u,v);
+
+  u = UnscaleU(u);
+  v = UnscaleV(v);
+
+  double tol = 1.e-12;
+  if ((u > - tol) && (u < 1. + tol) && (v > - tol) && (v < 1. + tol))
+    result = true;
+  else
+    result = false;
+
+  return result;
+}
+
+void FPatch::
+Dfdu(double u, double v, double &x, double &y, double &z)
+{
+  double px, py, pz, nx, ny, nz, d;
+  double pxu, pyu, pzu, nxu, nyu, nzu, du;
+
+  u = RescaleU(u);
+  v = RescaleV(v);
+
+  _ps->F(u,v,px,py,pz);
+  _ps->GetUnitNormal(u,v,nx,ny,nz);
+  _ps->Dfdu(u,v,pxu,pyu,pzu);
+  _ps->Dndu(u,v,nxu, nyu, nzu);
+
+  d = _Interpolate(u, v, 0, 0).real();
+  du = _Interpolate(u, v, 1, 0).real();
+
+  x = pxu + du * nx + d * nxu;
+  y = pyu + du * ny + d * nyu;
+  z = pzu + du * nz + d * nzu;
+}
+
+void FPatch::
+Dfdv(double u, double v, double &x, double &y, double &z)
+{
+  double px, py, pz, nx, ny, nz, d;
+  double pxv, pyv, pzv, nxv, nyv, nzv, dv;
+
+  u = RescaleU(u);
+  v = RescaleV(v);
+
+  _ps->F(u,v,px,py,pz);
+  _ps->GetUnitNormal(u,v,nx,ny,nz);
+  _ps->Dfdv(u,v,pxv,pyv,pzv);
+  _ps->Dndv(u,v,nxv,nyv,nzv);
+
+  d = _Interpolate(u, v, 0, 0).real();
+  dv = _Interpolate(u, v, 0, 1).real();
+
+  x = pxv + dv * nx + d * nxv;
+  y = pyv + dv * ny + d * nyv;
+  z = pzv + dv * nz + d * nzv;
+}
+
+void FPatch::
+Dfdfdudu(double u, double v, double &x, double &y, double &z)
+{
+  double px, py, pz, nx, ny, nz, d;
+  double pxu, pyu, pzu, nxu, nyu, nzu, du;
+  double pxuu, pyuu, pzuu, nxuu, nyuu, nzuu, duu;
+
+  u = RescaleU(u);
+  v = RescaleV(v);
+
+  _ps->F(u,v,px,py,pz);
+  _ps->GetUnitNormal(u,v,nx,ny,nz);
+  _ps->Dfdu(u,v,pxu,pyu,pzu);
+  _ps->Dndu(u,v,nxu,nyu,nzu);
+  _ps->Dfdfdudu(u,v,pxuu,pyuu,pzuu);
+  _ps->Dndndudu(u,v,nxuu,nyuu,nzuu);
+
+  d = _Interpolate(u, v, 0, 0).real();
+  du = _Interpolate(u, v, 1, 0).real();
+  duu = _Interpolate(u, v, 2, 0).real();
+
+  x = pxuu + duu * nx + du * nxu + du * nxu + d * nxuu;
+  y = pyuu + duu * ny + du * nyu + du * nyu + d * nyuu;
+  z = pzuu + duu * nz + du * nzu + du * nzu + d * nzuu;
+}
+
+void FPatch::
+Dfdfdvdv(double u, double v, double &x, double &y, double &z)
+{
+  double px, py, pz, nx, ny, nz, d;
+  double pxv, pyv, pzv, nxv, nyv, nzv, dv;
+  double pxvv, pyvv, pzvv, nxvv, nyvv, nzvv, dvv;
+
+  u = RescaleU(u);
+  v = RescaleV(v);
+
+  _ps->F(u,v,px,py,pz);
+  _ps->GetUnitNormal(u,v,nx,ny,nz);
+  _ps->Dfdv(u,v,pxv,pyv,pzv);
+  _ps->Dndv(u,v,nxv,nyv,nzv);
+  _ps->Dfdfdvdv(u,v,pxvv,pyvv,pzvv);
+  _ps->Dndndvdv(u,v,nxvv,nyvv,nzvv);
+
+  d = _Interpolate(u, v, 0, 0).real();
+  dv = _Interpolate(u, v, 0, 1).real();
+  dvv = _Interpolate(u, v, 0, 2).real();
+
+  x = pxvv + dvv * nx + dv * nxv + dv * nxv + d * nxvv;
+  y = pyvv + dvv * ny + dv * nyv + dv * nyv + d * nyvv;
+  z = pzvv + dvv * nz + dv * nzv + dv * nzv + d * nzvv;
+}
+
+void FPatch::
+Dfdfdudv(double u, double v, double &x, double &y, double &z)
+{
+  double px, py, pz, nx, ny, nz, d;
+  double pxu, pyu, pzu, nxu, nyu, nzu, du;
+  double pxv, pyv, pzv, nxv, nyv, nzv, dv;
+  double pxuv, pyuv, pzuv, nxuv, nyuv, nzuv, duv;
+
+  u = RescaleU(u);
+  v = RescaleV(v);
+
+  _ps->F(u,v,px,py,pz);
+  _ps->GetUnitNormal(u,v,nx,ny,nz);
+  _ps->Dfdu(u,v,pxu,pyu,pzu);
+  _ps->Dndu(u,v,nxu,nyu,nzu);
+  _ps->Dfdv(u,v,pxv,pyv,pzv);
+  _ps->Dndv(u,v,nxv,nyv,nzv);
+  _ps->Dfdfdudv(u,v,pxuv,pyuv,pzuv);
+  _ps->Dndndudv(u,v,nxuv,nyuv,nzuv);
+
+  d = _Interpolate(u, v, 0, 0).real();
+  du = _Interpolate(u, v, 1, 0).real();
+  dv = _Interpolate(u, v, 1, 0).real();
+  duv = _Interpolate(u, v, 1, 1).real();
+
+  x = pxuv + duv * nx + du * nxv + dv * nxu + d * nxuv;
+  y = pyuv + duv * ny + du * nyv + dv * nyu + d * nyuv;
+  z = pzuv + duv * nz + du * nzv + dv * nzu + d * nzuv;
+}
+
+double  FPatch::GetPou(double u, double v)
+{
+  double pouU, pouV;
+
+  if (_hardEdge[3])
+    pouU = OneSidedPartitionOfUnity(0.,1.,u);
+  else if (_hardEdge[1])
+    pouU = 1. - OneSidedPartitionOfUnity(0.,1.,u);
+  else
+    pouU = PartitionOfUnity(u, 0., 0.3, 0.7, 1.);
+
+  if (_hardEdge[0])
+    pouV = OneSidedPartitionOfUnity(0.,1.,v);
+  else if (_hardEdge[2])
+    pouV = 1. - OneSidedPartitionOfUnity(0.,1.,v);
+  else
+    pouV = PartitionOfUnity(v, 0., 0.3, 0.7, 1.);
+
+  return pouU * pouV;
+}
diff --git a/contrib/FourierModel/FPatch.h b/contrib/FourierModel/FPatch.h
new file mode 100644
index 0000000000000000000000000000000000000000..1a1abd2045d25e5292fe69d112fb765f20ac682d
--- /dev/null
+++ b/contrib/FourierModel/FPatch.h
@@ -0,0 +1,87 @@
+#ifndef _F_PATCH_H_
+#define _F_PATCH_H_
+
+#include <vector>
+#include <fftw3.h>
+#include <complex>
+#include "Patch.h"
+#include "PartitionOfUnity.h"
+#include "ProjectionSurface.h"
+
+// The base class for the patches
+class FPatch : public Patch {
+ protected:
+  // Data Size
+  int _nData;
+
+  // Number of Fourier Modes
+  int _uModes, _vModes, _nModes;
+
+  // Number of Modes in reprocessed Fourier/Chebyshev series
+  int _uM, _vM;
+
+   // Period Information
+  double _periodU, _periodV;
+
+  // Limits of Fourier Series
+  int _uModesLower, _uModesUpper, _vModesLower, _vModesUpper;
+
+  // Limits in the new series
+  int _uMLower, _uMUpper, _vMLower, _vMUpper;
+
+  // data (and its first 2 derivatives)
+  std::complex<double> **_coeffOriginalData;
+  std::complex<double> **_coeffData, **_coeffDerivU, **_coeffDerivV;
+  std::complex<double> **_coeffDerivUU, **_coeffDerivVV, **_coeffDerivUV;
+
+  // temporary interpolation variables
+  std::vector< std::complex<double> > _tmpOrigCoeff, _tmpOrigInterp;
+  std::vector< std::complex<double> > _tmpCoeff, _tmpInterp;
+
+  // polynomial evaluator
+  std::complex<double> _PolyEval(std::vector< std::complex<double> > _coeff,
+                                 std::complex<double> x);
+  // interpolation wrapper
+  std::complex<double> _Interpolate(double u,double v);
+  std::complex<double> _Interpolate(double u,double v,int uDer,int vDer);
+
+  // other internal functions
+  void _ReprocessSeriesCoeff();
+
+  // persistent fftw data (used to avoid recomputing the FFTW plans
+  // and reallocating memory every time)
+  static int _forwardSize, _backwardSize;
+  static fftw_plan _forwardPlan, _backwardPlan;
+  static fftw_complex *_forwardData, *_backwardData;
+  void _SetForwardPlan(int n);
+  void _SetBackwardPlan(int n);
+
+  // This routine computes the forward FFT (ignoring the last element
+  // of the data)
+  void _ForwardFft(int n, std::complex<double> *fftData);
+
+  // This routine computes the inverse FFT (ignoring the last element
+  // in the array), returns values in entire array (by using the
+  // periodic extension)
+  void _BackwardFft(int n, std::complex<double> *fftData);
+
+ public:
+  FPatch
+    (int tag, ProjectionSurface* ps,
+     std::vector<double> &u, std::vector<double> &v,
+     std::vector< std::complex<double> > &data,int derivative = 0);
+  virtual ~FPatch();
+
+  // Abstract functions of Patch
+
+  virtual double GetPou(double u, double v);
+  virtual void F(double u, double v, double &x, double &y, double &z);
+  virtual bool Inverse(double x,double y,double z,double &u,double &v);
+  virtual void Dfdu(double u, double v, double &x, double &y, double &z);
+  virtual void Dfdv(double u, double v, double &x, double &y, double &z);
+  virtual void Dfdfdudu(double u,double v,double &x,double &y,double &z);
+  virtual void Dfdfdudv(double u,double v,double &x,double &y,double &z);
+  virtual void Dfdfdvdv(double u,double v,double &x,double &y,double &z);
+};
+
+#endif
diff --git a/contrib/FourierModel/Interpolator2D.cpp b/contrib/FourierModel/Interpolator2D.cpp
index 422708376cdea7b66341a964477a5eea8a4bc301..e3fe9c6dad2864932e71bf6f9fa1f11cfc02f9c4 100755
--- a/contrib/FourierModel/Interpolator2D.cpp
+++ b/contrib/FourierModel/Interpolator2D.cpp
@@ -3,11 +3,17 @@
 #include <complex>
 #include "Interpolator2D.h"
 
+extern "C" {
+  void zgelss_(int &,int &,int &,std::complex<double> *,int &,
+	       std::complex<double> *,int &,double *,double &,int &,
+	       std::complex<double> *,int &,double *,int &);
+}
+
 FftPolyInterpolator2D::
 FftPolyInterpolator2D(const std::vector<double> &u, 
 		      const std::vector<double> &v,
-		      const std::vector< std::vector< std::complex<double> > > &data,
-		      int refineFactor, int polyOrder, int derivative)
+		      const std::vector< std::vector< std::complex<double> > > 
+		      &data, int refineFactor, int polyOrder, int derivative)
   : _refineFactor(refineFactor), _polyOrder(polyOrder), _derivative(derivative),
     _uFine(0), _vFine(0), _fineData(0), _fineDerivU(0), _fineDerivV(0), 
     _fineDerivUU(0), _fineDerivVV(0), _fineDerivUV(0)
@@ -457,3 +463,292 @@ std::complex<double> FftPolyInterpolator2D::Dfdfdudv(double u, double v)
 {
   return _Interpolate(u, v, 1, 1);
 }
+
+
+
+FourierContinuationInterpolator2D::FourierContinuationInterpolator2D
+(const std::vector<double> &u, const std::vector<double> &v,
+ const std::vector< std::complex<double> > &data, 
+ int derivative, int uModes, int vModes, double periodU, double periodV)
+  : _derivative(derivative), _uModes(uModes), _vModes(vModes), 
+    _periodU(periodU), _periodV(periodV),
+    _coeffData(0), _coeffDerivU(0), _coeffDerivV(0), _coeffDerivUU(0), 
+    _coeffDerivVV(0), _coeffDerivUV(0)
+{ 
+  // sanity check
+  if (!((u.size() == v.size()) && (v.size() == data.size()))) {
+    Msg::Fatal("Input data of unequal size.");
+  }
+  _nData = data.size();
+  _nModes = _uModes * _vModes;
+
+  for (int i=0;i<_nData;i++) {
+    _u.push_back(u[i]);
+    _v.push_back(v[i]);
+  }
+
+  // Initialize _Data to zero
+  _coeffData = new std::complex<double>*[_uModes];
+  for(int j = 0; j < _uModes; j++){
+    _coeffData[j] = new std::complex<double>[_vModes];
+    for(int k = 0; k < _vModes; k++)
+      _coeffData[j][k] = 0.;
+  }
+
+  if ((_uModes % 2) == 0) {
+    _uModesLower = - _uModes/2;
+    _uModesUpper = _uModes/2 - 1;
+  }
+  else {
+    _uModesLower = - (_uModes - 1)/2;
+    _uModesUpper = (_uModes - 1)/2;
+  }
+  if ((_vModes % 2) == 0) {
+    _vModesLower = - _vModes/2;
+    _vModesUpper = _vModes/2 - 1;
+  }
+  else {
+    _vModesLower = - (_vModes - 1)/2;
+    _vModesUpper = (_vModes - 1)/2;
+  }
+
+  std::complex<double> *LSRhs = new std::complex<double> [_nData];
+  std::complex<double> *LSMatrix = new std::complex<double> [_nModes * _nData];
+
+  for (int i=0;i<_nData;i++)
+    LSRhs[i] = data[i];
+
+  for (int j=0;j<_uModes;j++)
+    for (int k=0;k<_vModes;k++)
+      for (int i=0;i<_nData;i++)
+	LSMatrix[i+_nData*(k+_vModes*j)] = std::complex<double>
+	  (cos(2*M_PI * ((j + _uModesLower) * u[i] / _periodU + 
+			 (k + _vModesLower) * v[i] / _periodV)),
+	   sin(2*M_PI * ((j + _uModesLower) * u[i] / _periodU + 
+			 (k + _vModesLower) * v[i] / _periodV)));
+
+  // some parameters for the lease square solvers "zgelss"
+  int info, rank;
+  int lwork = 66*_nData, one = 1;
+  double rcond = -1.;
+  double *s = new double [_nModes];
+  double *rwork = new double [5*_nModes];
+  std::complex<double> *work = new std::complex<double> [lwork];
+
+  zgelss_(_nData,_nModes,one,&LSMatrix[0],_nData,&LSRhs[0],_nData,&s[0],rcond,
+	  rank,&work[0],lwork,&rwork[0],info);
+
+  for(int j = 0; j < _uModes; j++)
+    for(int k = 0; k < _vModes; k++)
+      _coeffData[j][k] = LSRhs[k+_vModes*j];
+
+  // deleting lease square arrays
+  delete[] s, rwork, work;
+  delete[] LSMatrix, LSRhs;
+
+  // Check if we need to interpolate the derivative(s)
+  if(_derivative){
+    // Initialize _fineDeriv and _fineDeriv2 to zero
+    if(_derivative & 1){
+      _coeffDerivU = new std::complex<double>*[_uModes];
+      _coeffDerivV = new std::complex<double>*[_uModes];
+      for(int j = 0; j < _uModes; j++){
+	_coeffDerivU[j] = new std::complex<double>[_vModes];
+	_coeffDerivV[j] = new std::complex<double>[_vModes];
+	for(int k = 0; k < _vModes; k++){
+	  _coeffDerivU[j][k] = 0.;
+	  _coeffDerivV[j][k] = 0.;
+	}
+      }
+    }
+    
+    if(_derivative & 2){
+      _coeffDerivUU = new std::complex<double>*[_uModes];
+      _coeffDerivVV = new std::complex<double>*[_uModes];
+      _coeffDerivUV = new std::complex<double>*[_uModes];
+      for(int j = 0; j < _uModes; j++){
+	_coeffDerivUU[j] = new std::complex<double>[_vModes];
+	_coeffDerivVV[j] = new std::complex<double>[_vModes];
+	_coeffDerivUV[j] = new std::complex<double>[_vModes];
+	for(int k = 0; k < _vModes; k++){
+	  _coeffDerivUU[j][k] = 0.;
+	  _coeffDerivVV[j][k] = 0.;
+	  _coeffDerivUV[j][k] = 0.;
+	}
+      }
+    }
+    
+    // Copy the Fourier coefficients into _coeffDeriv and _coeffDeriv2
+    std::complex<double> I(0., 1.);
+    for(int j = 0; j < _uModes; j++){
+      for(int k = 0; k < _vModes; k++){
+	int J = j+_uModesLower;
+	int K = k+_vModesLower;
+	if(_derivative & 1){
+	  _coeffDerivU[j][k] = (2 * M_PI * J * I / _periodU) *_coeffData[j][k];
+	  _coeffDerivV[j][k] = (2 * M_PI * K * I / _periodV) *_coeffData[j][k];
+	}
+	if(_derivative & 2){
+	  _coeffDerivUU[j][k] = - (4 * M_PI * M_PI * J * J / 
+				   (_periodU * _periodU)) * _coeffData[j][k];
+	  _coeffDerivVV[j][k] = - (4 * M_PI * M_PI * K * K /
+				   (_periodV * _periodV)) * _coeffData[j][k];
+	  _coeffDerivUV[j][k] = - (4 * M_PI * M_PI * J * K /
+				   (_periodU * _periodV)) * _coeffData[j][k];
+	}
+      }
+    }
+  }    
+  // Initialize interpolation variables
+  _tmpCoeff = std::vector< std::complex<double> >(_vModes);
+  _tmpInterp = std::vector< std::complex<double> >(_uModes);
+}
+
+FourierContinuationInterpolator2D::~FourierContinuationInterpolator2D()
+{
+  for(int j = 0; j < _uModes; j++)
+    delete [] _coeffData[j];
+  delete [] _coeffData;
+
+  if(_coeffDerivU){
+    for(int j = 0; j < _uModes; j++)
+      delete [] _coeffDerivU[j];
+    delete [] _coeffDerivU;
+  }
+  if(_coeffDerivV){
+    for(int j = 0; j < _uModes; j++)
+      delete [] _coeffDerivV[j];
+    delete [] _coeffDerivV;
+  }
+  if(_coeffDerivUU){
+    for(int j = 0; j < _uModes; j++)
+      delete [] _coeffDerivUU[j];
+    delete [] _coeffDerivUU;
+  }
+  if(_coeffDerivVV){
+    for(int j = 0; j < _uModes; j++)
+      delete [] _coeffDerivVV[j];
+    delete [] _coeffDerivVV;
+  }
+  if(_coeffDerivUV){
+    for(int j = 0; j < _uModes; j++)
+      delete [] _coeffDerivUV[j];
+    delete [] _coeffDerivUV;
+  }
+}
+
+std::complex<double> FourierContinuationInterpolator2D::
+_PolyEval(std::vector< std::complex<double> > _coeff, std::complex<double> x)
+{
+  int _polyOrder = _coeff.size()-1;
+  std::complex<double> out = 0.;
+
+  out = x * _coeff[_polyOrder];
+  for (int i = _polyOrder - 1; i > 0; i--)
+    out = x * (out + _coeff[i]);
+  out = out + _coeff[0];
+  
+  return out;
+}
+
+std::complex<double> FourierContinuationInterpolator2D::
+  _Interpolate(double u, double v, int uDer, int vDer)
+{
+  //Msg::Info("%d %d %d",uDer,vDer,_derivative);
+  if (((uDer==2 || vDer==2 || (uDer==1 && vDer==1)) && !(_derivative & 2) ) || 
+      ((uDer==1 || vDer==1) && !(_derivative & 1)) ||
+      (uDer<0 || uDer>2 || vDer<0 || vDer>2) ) {
+    Msg::Error("Derivative data not available: check contructor call %d %d %d",
+	       uDer,vDer,_derivative);
+    return 0.;
+  }
+
+  double epsilon = 1e-12;
+  if(u < 0. - epsilon || u > 1. + epsilon){
+    Msg::Error("Trying to interpolate outside interval: (u,v)=(%.16g,%.16g) "
+	       "not in [%g,%g]x[%g,%g]", u, v, 0., 1., 0., 1.);
+    return 0.;
+  }
+
+  // Interpolate to find value at (u,v)
+  for(int j = 0; j < _uModes; j++){
+    for(int k = 0; k < _vModes; k++){
+      std::complex<double> tmp;
+      if(uDer == 0 && vDer == 0)
+	tmp = _coeffData[j][k];
+      else if(uDer == 1 && vDer == 0)
+	tmp = _coeffDerivU[j][k];
+      else if(uDer == 0 && vDer == 1)
+	tmp = _coeffDerivV[j][k];
+      else if(uDer == 2 && vDer == 0)
+	tmp = _coeffDerivUU[j][k];
+      else if(uDer == 0 && vDer == 2)
+	tmp = _coeffDerivVV[j][k];
+      else
+	tmp = _coeffDerivUV[j][k];
+      _tmpCoeff[k] = tmp;
+    }
+    std::complex<double> y(cos(2 * M_PI * v / _periodV), 
+			   sin(2 * M_PI * v / _periodV));
+    _tmpInterp[j] = _PolyEval(_tmpCoeff, y);
+    _tmpInterp[j] *= std::complex<double>
+      (cos(2 * M_PI * _vModesLower * v / _periodV), 
+       sin(2 * M_PI * _vModesLower * v / _periodV));
+  }
+  std::complex<double> x(cos(2 * M_PI * u / _periodU), 
+			 sin(2 * M_PI * u / _periodU));
+  return _PolyEval(_tmpInterp, x) * std::complex<double>
+    (cos(2 * M_PI * _uModesLower * u / _periodU),
+     sin(2 * M_PI * _uModesLower * u / _periodU));
+}
+
+std::complex<double> FourierContinuationInterpolator2D::F(double u, double v)
+{
+  return _Interpolate(u, v, 0, 0);
+}
+
+std::complex<double> FourierContinuationInterpolator2D::Dfdu(double u, double v)
+{
+  return _Interpolate(u, v, 1, 0);
+}
+
+std::complex<double> FourierContinuationInterpolator2D::Dfdv(double u, double v)
+{
+  return _Interpolate(u, v, 0, 1);
+}
+
+std::complex<double> FourierContinuationInterpolator2D::Dfdfdudu(double u, double v)
+{
+  return _Interpolate(u, v, 2, 0);
+}
+
+std::complex<double> FourierContinuationInterpolator2D::Dfdfdvdv(double u, double v)
+{
+  return _Interpolate(u, v, 0, 2);
+}
+
+std::complex<double> FourierContinuationInterpolator2D::Dfdfdudv
+(double u, double v)
+{
+  return _Interpolate(u, v, 1, 1);
+}
+
+bool FourierContinuationInterpolator2D::IsPointInInterpolationRange
+(double u, double v)
+{
+  double nbdRadius = 1.e-2;
+  bool status = false;
+  for (int i=0;i<_nData;i++) {
+    double tmp = sqrt((u - _u[i])*(u - _u[i])+(v - _v[i])*(v - _v[i]));
+    if (tmp < nbdRadius) {
+      status = true;
+      break;
+    }
+  }
+  return status;
+}
+
+int FourierContinuationInterpolator2D::GetDataSize()
+{
+  return _nData;
+}
diff --git a/contrib/FourierModel/Interpolator2D.h b/contrib/FourierModel/Interpolator2D.h
index df44f4c7d3db85aa06230e42ccdc829d8978954d..1b55eb5fc156b74dabbda72a383af2967f1b88cd 100755
--- a/contrib/FourierModel/Interpolator2D.h
+++ b/contrib/FourierModel/Interpolator2D.h
@@ -43,6 +43,18 @@ class Interpolator2D {
     Msg::Error("Second derivative not implemented for this interpolator");
     return 0.; 
   }
+  // checks if the interpolation point is good enough
+  virtual bool IsPointInInterpolationRange(double u, double v)
+  {
+    Msg::Error("Goodness check not implemented for this interpolator");
+    return 0; 
+  }
+  // returns the size of interpolation data
+  virtual int GetDataSize()
+  {
+    Msg::Error("GetDataSize not implemented for this interpolator");
+    return 0;
+  }
 };
 
 // FFT + polynomial interpolation on refined grid (assumes that the
@@ -95,8 +107,9 @@ class FftPolyInterpolator2D : public Interpolator2D {
  public:
   FftPolyInterpolator2D(const std::vector<double> &u,
 			const std::vector<double> &v,
-			const std::vector< std::vector< std::complex<double> > > &data,
-			int refineFactor=16, int polyOrder=3, int derivative=0);
+			const std::vector<std::vector<std::complex<double> > >
+			&data, int refineFactor=16, int polyOrder=3, 
+			int derivative=0);
   ~FftPolyInterpolator2D();
   virtual std::complex<double> F(double u, double v);
   virtual std::complex<double> Dfdu(double u, double v);
@@ -106,4 +119,48 @@ class FftPolyInterpolator2D : public Interpolator2D {
   virtual std::complex<double> Dfdfdudv(double u, double v);
 };
 
+// Fourier Continuation interpolation
+class FourierContinuationInterpolator2D : public Interpolator2D {
+ private:
+  friend class Body;
+  // u and v data
+  std::vector<double> _u, _v;
+  // temporary interpolation variables
+  std::vector< std::complex<double> > _tmpCoeff, _tmpInterp;
+ protected:
+  // bitfield telling if we also interpolate the derivative(s)
+  int _derivative;
+  // Data Size
+  int _nData;
+  // Number of Fourier Modes
+  int _uModes, _vModes, _nModes;
+  // Period Information
+  double _periodU, _periodV;
+  // Limits of Fourier Series
+  int _uModesLower, _uModesUpper, _vModesLower, _vModesUpper;
+  // data (and its first 2 derivatives)
+  std::complex<double> **_coeffData, **_coeffDerivU, **_coeffDerivV;
+  std::complex<double> **_coeffDerivUU, **_coeffDerivVV, **_coeffDerivUV;
+  // polynomial evaluator
+  std::complex<double> _PolyEval(std::vector< std::complex<double> > _coeff, 
+				 std::complex<double> x);
+  // interpolation wrapper
+  std::complex<double> _Interpolate(double u,double v,int uDer=0,int vDer=0);
+ public:
+  FourierContinuationInterpolator2D
+    (const std::vector<double> &u, const std::vector<double> &v,
+     const std::vector< std::complex<double> > &data,
+     int derivative = 0, int uModes = 1, int vModes = 1, 
+     double periodU = 2, double periodV = 2);
+  ~FourierContinuationInterpolator2D();
+  virtual std::complex<double> F(double u, double v);
+  virtual std::complex<double> Dfdu(double u, double v);
+  virtual std::complex<double> Dfdv(double u, double v);
+  virtual std::complex<double> Dfdfdudu(double u, double v);
+  virtual std::complex<double> Dfdfdvdv(double u, double v);
+  virtual std::complex<double> Dfdfdudv(double u, double v);
+  virtual bool IsPointInInterpolationRange(double u, double v);
+  virtual int GetDataSize();
+};
+
 #endif
diff --git a/contrib/FourierModel/Makefile b/contrib/FourierModel/Makefile
index f8de9d751edb21c6de9e8eff87724f0deafbb8d5..917b74c370ea44ce3e99e24ebe3c5ae618c196b3 100644
--- a/contrib/FourierModel/Makefile
+++ b/contrib/FourierModel/Makefile
@@ -7,12 +7,16 @@ CFLAGS = ${OPTIM} ${FLAGS}
 SRC = ProjectionSurface.cpp \
 	CylindricalProjectionSurface.cpp \
 	RevolvedParabolaProjectionSurface.cpp \
-      ContinuationPatch.cpp \
-      ExactPatch.cpp\
-      Patch.cpp\
+      Patch.h \
+	FPatch.cpp \
+	ContinuationPatch.cpp \
+	ExactPatch.cpp\
+      Curve.cpp \
+	FCurve.cpp\
+	IntersectionCurve.cpp \
+      Point.cpp \
       BlendedPatch.cpp \
       BlendOperator.cpp \
-      Curve.cpp\
       FM_Edge.cpp\
       FM_Face.cpp\
       FM_Info.cpp\
@@ -46,36 +50,3 @@ depend:
 	rm -f Makefile.new
 
 # DO NOT DELETE THIS LINE
-ProjectionSurface.o: ProjectionSurface.cpp ProjectionSurface.h
-CylindricalProjectionSurface.o: CylindricalProjectionSurface.cpp \
-  CylindricalProjectionSurface.h ProjectionSurface.h
-RevolvedParabolaProjectionSurface.o:  \
- RevolvedParabolaProjectionSurface.cpp \
-  RevolvedParabolaProjectionSurface.h Utils.h ProjectionSurface.h
-ContinuationPatch.o: ContinuationPatch.cpp Message.h ContinuationPatch.h \
-  Patch.h FM_Info.h ProjectionSurface.h PartitionOfUnity.h \
-  Interpolator1D.h
-ExactPatch.o: ExactPatch.cpp Message.h ExactPatch.h Patch.h FM_Info.h \
-  ProjectionSurface.h
-Patch.o: Patch.cpp Patch.h FM_Info.h ProjectionSurface.h
-BlendedPatch.o: BlendedPatch.cpp BlendedPatch.h Message.h Patch.h \
-  FM_Info.h ProjectionSurface.h BlendOperator.h PartitionOfUnity.h
-BlendOperator.o: BlendOperator.cpp BlendOperator.h FM_Info.h Patch.h \
-  ProjectionSurface.h
-Curve.o: Curve.cpp Message.h Curve.h Patch.h FM_Info.h \
-  ProjectionSurface.h
-FM_Edge.o: FM_Edge.cpp FM_Edge.h Curve.h Patch.h FM_Info.h \
-  ProjectionSurface.h FM_Vertex.h Message.h
-FM_Face.o: FM_Face.cpp FM_Face.h Patch.h FM_Info.h ProjectionSurface.h \
-  FM_Edge.h Curve.h FM_Vertex.h Message.h
-FM_Info.o: FM_Info.cpp FM_Info.h
-FM_Reader.o: FM_Reader.cpp Message.h FM_Reader.h Curve.h Patch.h \
-  FM_Info.h ProjectionSurface.h ExactPatch.h ContinuationPatch.h \
-  PartitionOfUnity.h Interpolator1D.h CylindricalProjectionSurface.h \
-  RevolvedParabolaProjectionSurface.h Utils.h FM_Face.h FM_Edge.h \
-  FM_Vertex.h BlendOperator.h BlendedPatch.h
-FM_Vertex.o: FM_Vertex.cpp
-Message.o: Message.cpp Message.h
-Utils.o: Utils.cpp Utils.h
-PartitionOfUnity.o: PartitionOfUnity.cpp PartitionOfUnity.h
-Interpolator1D.o: Interpolator1D.cpp Interpolator1D.h Message.h
diff --git a/contrib/FourierModel/Patch.cpp b/contrib/FourierModel/Patch.cpp
index 6d9ef939f19b4101ba8dc2c9dc3fa514993af092..0db5438c8c4705cbd75f8b93246b78e43d2329c1 100644
--- a/contrib/FourierModel/Patch.cpp
+++ b/contrib/FourierModel/Patch.cpp
@@ -1,18 +1,9 @@
 #include <cmath>
 #include "Patch.h"
 
-Patch::Patch() :_PI(0), _ps(0), _uMin(0.), _uMax(1.), _vMin(0.), _vMax(1.),
-		_periodicityU(0), _periodicityV(0), _derivative(0) {}
-
-Patch::Patch(PatchInfo* PI) : _PI(PI), _ps(0), _derivative(0) 
-{
-  _uMin = _PI->uMin;
-  _uMax = _PI->uMax;
-  _vMin = _PI->vMin;
-  _vMax = _PI->vMax;
-  _periodicityU = _PI->periodic[0];
-  _periodicityV = _PI->periodic[1];
-}
+Patch::Patch() :_ps(0), _uMin(0.), _uMax(1.), _vMin(0.), _vMax(1.),
+		_periodicityU(0), _periodicityV(0), _derivative(0), 
+		_tag(-1) {}
 
 void Patch::GetNormal(double u, double v, double &x, double &y, double &z)
 {
@@ -90,8 +81,3 @@ void Patch::Dndv(double u, double v, double &x, double &y, double &z)
   z = ((normSquared - n[2]*n[2])*dndv[2] - n[2]*n[0]*dndv[0] -
        n[2]*n[1]*dndv[1]) / normCubed;
 }
-
-int Patch::GetTag()
-{
-  return _PI->tag;
-}
diff --git a/contrib/FourierModel/Patch.h b/contrib/FourierModel/Patch.h
index e798422362091f1b69561c432f9534f1f0ba7d25..80c3ad898f1d4b6e25ea9b9716ac4b1169a9009e 100644
--- a/contrib/FourierModel/Patch.h
+++ b/contrib/FourierModel/Patch.h
@@ -2,25 +2,25 @@
 #define _PATCH_H_
 
 #include <cmath>
-#include "FM_Info.h"
 #include "ProjectionSurface.h"
 
 // The base class for the patches
 class Patch {
  protected:
   // bitfield telling if we also interpolate the derivative(s)
+  int _tag;
   int _derivative;
   double _uMin, _uMax;
   double _vMin, _vMax;
   int _periodicityU, _periodicityV;
   ProjectionSurface* _ps;
+  // Hard edges
+  int _hardEdge[4];
  public:
-  PatchInfo* _PI;
   Patch();
-  Patch(PatchInfo* PI);
   virtual ~Patch() {}
 
-  int GetTag();
+  inline int GetTag() { return _tag; }
 
   inline double GetMinU() { return _uMin; }
   inline double GetMaxU() { return _uMax; }
@@ -36,6 +36,9 @@ class Patch {
   inline bool IsVPeriodic() 
     { if (_periodicityV) return true; else return false; }
 
+  inline bool IsHardEdge(int i)
+    { if (_hardEdge[i]) return true; else return false; }
+
   inline void SetPeriodicity(int pU, int pV) 
     { _periodicityU = pU; _periodicityV = pV; }
 
@@ -81,10 +84,7 @@ class Patch {
       return v;
     }
 
-  inline int GetUModes() { return _PI->nM[0]; }
-  inline int GetVModes() { return _PI->nM[1]; }
   inline int GetDerivativeBitField() { return _derivative; }
-
   inline ProjectionSurface* GetProjectionSurface() { return _ps; }
   inline void SetProjectionSurface(ProjectionSurface* ps) { _ps = ps; }
 
diff --git a/contrib/FourierModel/ProjectionSurface.h b/contrib/FourierModel/ProjectionSurface.h
index 72a5573ad0976bc5bfe1a5e26d3437cd11fb5bf4..0bd3c149aca693f8e6a5bf08ba04a8a0d014c40a 100755
--- a/contrib/FourierModel/ProjectionSurface.h
+++ b/contrib/FourierModel/ProjectionSurface.h
@@ -26,6 +26,9 @@ class ProjectionSurface {
 
   inline int GetNumParameters() { return numParameters_; }
 
+  inline bool IsUPeriodic() { return uPeriodic_; }
+  inline bool IsVPeriodic() { return vPeriodic_; }
+
   // Public access functions
 
   inline void GetOrigin
@@ -85,6 +88,8 @@ class ProjectionSurface {
     (int i, double x) = 0;
   virtual double GetParameter
     (int i) = 0;
+  virtual std::string GetLabel
+    (int i) = 0;
 
   // These functions may also be provided by the derived 
   // projection surfaces (usually for better performance), 
diff --git a/contrib/FourierModel/RevolvedParabolaProjectionSurface.cpp b/contrib/FourierModel/RevolvedParabolaProjectionSurface.cpp
index 1e55997506c3aafafd404d226247e4a2ad2b4dc2..caaa55106120e2d7e687abd2fc950883a5640189 100755
--- a/contrib/FourierModel/RevolvedParabolaProjectionSurface.cpp
+++ b/contrib/FourierModel/RevolvedParabolaProjectionSurface.cpp
@@ -264,3 +264,16 @@ GetParameter(int i)
     return K_[1];
   }
 }
+
+std::string RevolvedParabolaProjectionSurface::
+GetLabel(int i)
+{
+  switch (i) {
+  case 0:
+    return std::string("Radius");
+  case 1:
+    return std::string("Parabola Scale 1");
+  case 2:
+    return std::string("Parabola Scale 2");
+  }
+}
diff --git a/contrib/FourierModel/RevolvedParabolaProjectionSurface.h b/contrib/FourierModel/RevolvedParabolaProjectionSurface.h
index 73e7a861ae02b1498ed0aa6c47e5a272d50ca733..1beb120659f01b251b63c67f00ea1e193b64ee14 100755
--- a/contrib/FourierModel/RevolvedParabolaProjectionSurface.h
+++ b/contrib/FourierModel/RevolvedParabolaProjectionSurface.h
@@ -50,6 +50,8 @@ class RevolvedParabolaProjectionSurface : public ProjectionSurface {
     (int i, double x);
   virtual double GetParameter
     (int i);
+  virtual std::string GetLabel
+    (int i);
 
   // Redefinitions for RevolvedParabolaProjectionSurface
 
diff --git a/contrib/FourierModel/Utils.cpp b/contrib/FourierModel/Utils.cpp
index 956c8f8d9e1927d829b1056a5e755da0856ecf57..344a4e209312084abbe8802c418245e5b676f87a 100755
--- a/contrib/FourierModel/Utils.cpp
+++ b/contrib/FourierModel/Utils.cpp
@@ -1,4 +1,5 @@
 #include "Utils.h"
+#include "Message.h"
 
 std::vector<double> SolveCubic(double a, double b, double c)
 {
@@ -80,8 +81,178 @@ std::vector<double> SolveCubic(double a, double b, double c)
     }
   }
 
-  //  cout << a << ' ' << b << ' ' << c << ' ' <<
-  //     num << ' ' << root[0] << ' ' <<root[1] << ' ' << root[2] << endl; 
-
   return root;
 }
+
+void  find(std::vector<int> &a, int length, std::vector<int> &q, int &num){
+
+  int i,j=0;
+
+  for (i=0 ; i < length ; i++){
+    if (a[i]!=0){
+      q[j]=i+1;
+      j=j+1;
+    }
+  }
+  num=j;
+}
+
+int minVec(std::vector<int> &a,int n){
+//   vector<int> a(n);
+//   copyVec(aa,n,a);
+  int i;
+  int s=a[0];
+
+  for (i=1;i<n ;i++){
+    if (s>a[i]){
+      s=a[i];
+    }
+  }
+  return s;
+}
+
+int maxVec(std::vector<int> &a,int n){
+//   vector<int> a(n);
+//   copyVec(aa,n,a);
+  int i;
+  int s=a[0];
+
+  for (i=1;i<n ;i++){
+    if (s<a[i]){
+      s=a[i];
+    }
+  }
+  return s;
+}
+
+std::vector<std::vector<int> > ones(int row, int col){
+  std::vector<std::vector<int> > tmp(row, std::vector<int> (col));
+  for (int i=0 ; i < row ; i++){
+    for (int j=0 ; j < col ; j++){
+      tmp[i][j]=1;
+    }
+  }
+  return tmp;
+}
+
+void plotSceneViewer(int app, char* bffer, std::vector<int> &color,
+                      std::vector<std::vector<double> > &x,
+                      std::vector<std::vector<double> > &y,
+                      std::vector<std::vector<double> > &z, int ROW, int COL,
+                      std::vector<std::vector<int> > &mask){
+
+  std::fstream outfile;
+  if (app==0){
+    outfile.open(bffer, std::ios::out);
+  }else{
+    outfile.open(bffer, std::ios::out | std::ios::app);
+  }
+
+  if (outfile.fail())
+    {
+      Msg::Error("Could not open readfile.txt");
+      exit(1);
+    }
+  if (app==0){//not appending
+    outfile  << "#Inventor V2.1 ascii" << std::endl;
+    outfile  << "#created by allplot.cpp" << std::endl;
+  }
+
+  outfile  << "Separator {" << std::endl;
+  outfile  << "  Material {" << std::endl;
+  outfile  << "    diffuseColor [" << std::endl;
+  outfile  << "      " << ' ' << color[0] << ' ' << color[1] << ' ' << color[2] << std::endl;
+  outfile  << "    ]" << std::endl;
+  outfile  << "  }" << std::endl;
+  outfile  << "  IndexedTriangleStripSet {" << std::endl;;
+  outfile  << "    vertexProperty VertexProperty {"<< std::endl;;
+  outfile  << "      vertex [" << std::endl;
+  //outfile  << "        " << ' ' << x(0,0) << ' ' << y(0,0) << ' ' << z(0,0);
+
+  int i, j;
+
+  for (j=0 ; j < COL ; j++){
+    for (i=0 ; i < ROW ; i++){
+      if (i==0 & j==0 ){
+          outfile  << "        " << ' ' << x[i][j] << ' ' << y[i][j] << ' ' << z[i][j];
+      }else{
+        outfile  <<"," << ' ' << ' ' << std::endl;
+        outfile  << "        " << ' ' << x[i][j] << ' ' << y[i][j] << ' ' << z[i][j];
+      }
+    }
+  }
+
+  outfile  <<' ' << std::endl;
+  outfile  << "      ]" << std::endl;
+  outfile  << "    }" << std::endl;
+  outfile  << "    coordIndex [" << std::endl;
+
+  int jump=ROW;
+  int col, row;
+  int startrow, startoffset, endrow, endoffset;
+  std::vector<int> qq(ROW), qw(ROW);
+  int num1,num2,min1,min2,max1,max2;
+
+  //  for (col=0 ; col < ROW-1 ; col++){
+  for (col=0 ; col < COL-1 ; col++){
+    //copyColToVec(mask, ROW, col, qq);
+    for (int i=0 ; i < ROW ; i++){
+      qq[i]=mask[i][col];
+    }
+
+    find(qq, ROW, qw,  num1);//cout << num1 << endl;
+
+    if (num1!=0){
+      min1=minVec(qw,num1);
+      max1=maxVec(qw,num1);
+    }
+
+    //  copyColToVec(mask, ROW, col+1, qq);
+    for (int i=0 ; i < ROW ; i++){
+      qq[i]=mask[i][col+1];
+    }
+
+    find(qq, ROW, qw,  num2);
+    if (num2!=0){
+      min2=minVec(qw,num2);
+      max2=maxVec(qw,num2);
+    }
+
+    if (num1!=0 & num2!=0){
+      startrow=std::max(min1,min2)-1;
+      startoffset=min1-min2;
+
+      if (startoffset>0){
+        outfile  << "     " << ' ' << (col+1)*jump+startrow << ' ' <<"," << ' ' <<
+          col*jump+startrow << ' ' <<"," << ' ' << (col+1)*jump+startrow-1  << ' '<<"," << ' ' << "-1,"<< std::endl;
+      }else if (startoffset<0){
+        outfile  << "     " << ' ' << (col+1)*jump+startrow << ' ' <<"," << ' '<< 
+          col*jump+startrow << ' ' <<"," << ' ' << col*jump+startrow-1 <<' ' << "," << ' ' << "-1,"<< std::endl;
+      }
+
+      endrow=std::min(max1,max2)-1;
+      endoffset=max1-max2;
+
+      if (endoffset>0){
+        outfile  << "     " << ' ' <<(col+1)*jump+endrow<< ' ' <<"," << ' ' <<
+          col*jump+endrow+1 << ' ' <<"," << ' ' << col*jump+endrow << ' '<< "," << ' ' << "-1,"<<std::endl;
+      }else if (endoffset<0){
+        outfile  << "     " << ' ' <<(col+1)*jump+endrow << ' ' << "," << ' ' <<
+          (col+1)*jump+endrow+1 << ' ' <<"," << ' ' << col*jump+endrow <<' ' <<"," << ' ' << "-1,"<< std::endl;
+      }
+
+      outfile << "     " << ' ' << col*jump+startrow << ' ' <<"," << ' ' << (col+1)*jump+startrow<< ' ' << "," ;
+
+      for (row= startrow+1; row <= endrow ; row++){
+        outfile << col*jump+row << ' '  <<"," << ' '<< (col+1)*jump+row<< ' ' << ",";
+      }
+    }
+      outfile << " -1," << std::endl;
+  }
+
+  outfile  << "    ]" << std::endl;
+  outfile  << "  }" << std::endl;
+  outfile  << "}" << std::endl;
+
+  outfile.close();
+}
diff --git a/contrib/FourierModel/Utils.h b/contrib/FourierModel/Utils.h
index 9a51797d7e5f154d8dd84d19dc180235ac97974a..7fbfabbcd50d327fc91e045e691d5a3756ea6cc0 100755
--- a/contrib/FourierModel/Utils.h
+++ b/contrib/FourierModel/Utils.h
@@ -3,7 +3,18 @@
 
 #include <cmath>
 #include <vector>
+#include <fstream>
+#include <iostream>
 
 std::vector<double> SolveCubic(double a, double b, double c);
+void  find(std::vector<int> &a, int length, std::vector<int> &q, int &num);
+int minVec(std::vector<int> &a,int n);
+int maxVec(std::vector<int> &a,int n);
+std::vector<std::vector<int> > ones(int row, int col);
+void plotSceneViewer(int app, char* bffer, std::vector<int>& color,
+                      std::vector<std::vector<double> > &x,
+                      std::vector<std::vector<double> > &y,
+                      std::vector<std::vector<double> > &z, int ROW, int COL,
+                      std::vector<std::vector<int> > &mask);
 
 #endif