diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index a0698fbf2f14d3030010c49e562561283027e320..423db30402ad3f43180fc5fcde7fb1d7f7846add 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -445,12 +445,12 @@ void GEdge::replaceEndingPoints(GVertex *replOfv0, GVertex *replOfv1)
 {
   replaceEndingPointsInternals(replOfv0, replOfv1);
   if (replOfv0 != v0){
-    v0->delEdge(this);
+    if (v0) v0->delEdge(this);
     replOfv0->addEdge(this);
     v0 = replOfv0;
   }
   if (replOfv1 != v1){
-    v1->delEdge(this);
+    if (v1) v1->delEdge(this);
     replOfv1->addEdge(this);
     v1 = replOfv1;
   }
diff --git a/Geo/GEdge.h b/Geo/GEdge.h
index a9867ba24f9a3914239bdefb1b3a0cfa4b8a9bd9..5260301d8859d92abc17e8d7209b6deea64e7458 100644
--- a/Geo/GEdge.h
+++ b/Geo/GEdge.h
@@ -49,8 +49,8 @@ class GEdge : public GEntity {
   virtual void deleteMesh();
 
   // get the start/end vertices of the edge
-  GVertex *getBeginVertex() const { return v0; }
-  GVertex *getEndVertex() const { return v1; }
+  virtual GVertex *getBeginVertex() const { return v0; }
+  virtual GVertex *getEndVertex() const { return v1; }
 
   void reverse();
 
@@ -207,7 +207,7 @@ class GEdge : public GEntity {
 
   std::vector<MLine*> lines;
 
-  void addLine(MLine *line){ lines.push_back(line); }
+  virtual void addLine(MLine *line){ lines.push_back(line); }
 
   virtual void discretize(double tol, std::vector<SPoint3> &dpts, std::vector<double> &ts);
   SPoint3 closestPoint (SPoint3 &p, double tolerance);
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 0d105850b85a791fbd44512926a19d2c34c83c84..2363410532526632a129da4271a4414e2844a171 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -22,7 +22,7 @@
 #if defined(HAVE_MESH)
 #include "meshGFaceOptimize.h"
 #include "meshGFaceLloyd.h"
-#include "BackgroundMesh.h"
+#include "BackgroundMeshTools.h"
 #endif
 
 #if defined(HAVE_BFGS)
diff --git a/Geo/GModelIO_GEO.cpp b/Geo/GModelIO_GEO.cpp
index fa732069f8c6b839c7d6f32bb39aeb8e8a36f094..fea9a70242a536832d5222ad1229bae977a76874 100644
--- a/Geo/GModelIO_GEO.cpp
+++ b/Geo/GModelIO_GEO.cpp
@@ -58,14 +58,14 @@ int GModel::exportDiscreteGEOInternals()
 
   for(viter it = firstVertex(); it != lastVertex(); it++){
     Vertex *v = Create_Vertex((*it)->tag(), (*it)->x(), (*it)->y(), (*it)->z(),
-                              (*it)->prescribedMeshSizeAtVertex(), 1.0);
+        (*it)->prescribedMeshSizeAtVertex(), 1.0);
     Tree_Add(_geo_internals->Points, &v);
   }
 
   for(eiter it = firstEdge(); it != lastEdge(); it++){
     if((*it)->geomType() == GEntity::DiscreteCurve){
       Curve *c = Create_Curve((*it)->tag(), MSH_SEGM_DISCRETE, 1,
-                              NULL, NULL, -1, -1, 0., 1.);
+          NULL, NULL, -1, -1, 0., 1.);
       List_T *points = Tree2List(_geo_internals->Points);
       GVertex *gvb = (*it)->getBeginVertex();
       GVertex *gve = (*it)->getEndVertex();
@@ -147,7 +147,7 @@ int GModel::importGEOInternals()
       Curve *c;
       List_Read(curves, i, &c);
       if(c->Num >= 0){
-	GEdge *e = getEdgeByTag(c->Num);
+        GEdge *e = getEdgeByTag(c->Num);
         if(!e && c->Typ == MSH_SEGM_COMPOUND){
           std::vector<GEdge*> comp;
           for(unsigned int j = 0; j < c->compound.size(); j++){
@@ -155,18 +155,18 @@ int GModel::importGEOInternals()
             if(ge) comp.push_back(ge);
           }
           e = new GEdgeCompound(this, c->Num, comp);
-	  e->meshAttributes.method = c->Method;
-	  e->meshAttributes.nbPointsTransfinite = c->nbPointsTransfinite;
-	  e->meshAttributes.typeTransfinite = c->typeTransfinite;
-	  e->meshAttributes.coeffTransfinite = c->coeffTransfinite;
-	  e->meshAttributes.extrude = c->Extrude;
-	  e->meshAttributes.reverseMesh = c->ReverseMesh;
+          e->meshAttributes.method = c->Method;
+          e->meshAttributes.nbPointsTransfinite = c->nbPointsTransfinite;
+          e->meshAttributes.typeTransfinite = c->typeTransfinite;
+          e->meshAttributes.coeffTransfinite = c->coeffTransfinite;
+          e->meshAttributes.extrude = c->Extrude;
+          e->meshAttributes.reverseMesh = c->ReverseMesh;
           add(e);
         }
         else if(!e && c->beg && c->end){
           e = new gmshEdge(this, c,
-                           getVertexByTag(c->beg->Num),
-                           getVertexByTag(c->end->Num));
+              getVertexByTag(c->beg->Num),
+              getVertexByTag(c->end->Num));
           add(e);
         }
         else if(!e){
@@ -199,24 +199,24 @@ int GModel::importGEOInternals()
           if(gf)
             comp.push_back(gf);
         }
-	std::list<GEdge*> b[4];
+        std::list<GEdge*> b[4];
         for(int j = 0; j < 4; j++){
-	  for(unsigned int k = 0; k < s->compoundBoundary[j].size(); k++){
-	    GEdge *ge = getEdgeByTag(s->compoundBoundary[j][k]);
-	    if(ge) b[j].push_back(ge);
-	  }
-	}
+          for(unsigned int k = 0; k < s->compoundBoundary[j].size(); k++){
+            GEdge *ge = getEdgeByTag(s->compoundBoundary[j][k]);
+            if(ge) b[j].push_back(ge);
+          }
+        }
         int param = CTX::instance()->mesh.remeshParam;
-	GFaceCompound::typeOfCompound typ = GFaceCompound::HARMONIC_CIRCLE;
-	if (param == 1) typ =  GFaceCompound::CONFORMAL_SPECTRAL;
-	if (param == 2) typ =  GFaceCompound::RADIAL_BASIS;
-	if (param == 3) typ =  GFaceCompound::HARMONIC_PLANE;
-	if (param == 4) typ =  GFaceCompound::CONVEX_CIRCLE;
-	if (param == 5) typ =  GFaceCompound::CONVEX_PLANE;
-	if (param == 6) typ =  GFaceCompound::HARMONIC_SQUARE;
-	if (param == 7) typ =  GFaceCompound::CONFORMAL_FE;
+        GFaceCompound::typeOfCompound typ = GFaceCompound::HARMONIC_CIRCLE;
+        if (param == 1) typ =  GFaceCompound::CONFORMAL_SPECTRAL;
+        if (param == 2) typ =  GFaceCompound::RADIAL_BASIS;
+        if (param == 3) typ =  GFaceCompound::HARMONIC_PLANE;
+        if (param == 4) typ =  GFaceCompound::CONVEX_CIRCLE;
+        if (param == 5) typ =  GFaceCompound::CONVEX_PLANE;
+        if (param == 6) typ =  GFaceCompound::HARMONIC_SQUARE;
+        if (param == 7) typ =  GFaceCompound::CONFORMAL_FE;
         int algo = CTX::instance()->mesh.remeshAlgo;
-	f = new GFaceCompound(this, s->Num, comp, b[0], b[1], b[2], b[3], typ, algo);
+        f = new GFaceCompound(this, s->Num, comp, b[0], b[1], b[2], b[3], typ, algo);
         f->meshAttributes.recombine = s->Recombine;
         f->meshAttributes.recombineAngle = s->RecombineAngle;
         f->meshAttributes.method = s->Method;
@@ -320,14 +320,14 @@ int GModel::importGEOInternals()
       GEntity *ge = 0;
       int tag = CTX::instance()->geom.orientedPhysicals ? abs(num) : num;
       switch(p->Typ){
-      case MSH_PHYSICAL_POINT:   ge = getVertexByTag(tag); break;
-      case MSH_PHYSICAL_LINE:    ge = getEdgeByTag(tag); break;
-      case MSH_PHYSICAL_SURFACE: ge = getFaceByTag(tag); break;
-      case MSH_PHYSICAL_VOLUME:  ge = getRegionByTag(tag); break;
+        case MSH_PHYSICAL_POINT:   ge = getVertexByTag(tag); break;
+        case MSH_PHYSICAL_LINE:    ge = getEdgeByTag(tag); break;
+        case MSH_PHYSICAL_SURFACE: ge = getFaceByTag(tag); break;
+        case MSH_PHYSICAL_VOLUME:  ge = getRegionByTag(tag); break;
       }
       int pnum = CTX::instance()->geom.orientedPhysicals ? (sign(num) * p->Num) : p->Num;
       if(ge && std::find(ge->physicals.begin(), ge->physicals.end(), pnum) ==
-         ge->physicals.end())
+          ge->physicals.end())
         ge->physicals.push_back(pnum);
     }
   }
@@ -342,7 +342,7 @@ int GModel::importGEOInternals()
     }
   }
   for(std::map<int,int>::iterator it = _geo_internals->periodicFaces.begin();
-       it != _geo_internals->periodicFaces.end(); ++it){
+      it != _geo_internals->periodicFaces.end(); ++it){
     GFace *gf = getFaceByTag(abs(it->first));
     if (gf)gf->setMeshMaster(it->second * (it->first > 0 ? 1 : -1));
   }
@@ -369,79 +369,79 @@ int GModel::importGEOInternals()
 }
 
 class writeFieldOptionGEO {
- private :
-  FILE *geo;
-  Field *field;
- public :
-  writeFieldOptionGEO(FILE *fp,Field *_field) { geo = fp ? fp : stdout; field=_field; }
-  void operator() (std::pair<std::string, FieldOption *> it)
-  {
-    std::string v;
-    it.second->getTextRepresentation(v);
-    fprintf(geo, "Field[%i].%s = %s;\n", field->id, it.first.c_str(), v.c_str());
-  }
+  private :
+    FILE *geo;
+    Field *field;
+  public :
+    writeFieldOptionGEO(FILE *fp,Field *_field) { geo = fp ? fp : stdout; field=_field; }
+    void operator() (std::pair<std::string, FieldOption *> it)
+    {
+      std::string v;
+      it.second->getTextRepresentation(v);
+      fprintf(geo, "Field[%i].%s = %s;\n", field->id, it.first.c_str(), v.c_str());
+    }
 };
 
 class writeFieldGEO {
- private :
-  FILE *geo;
- public :
-  writeFieldGEO(FILE *fp) { geo = fp ? fp : stdout; }
-  void operator() (std::pair<const int, Field *> it)
-  {
-    fprintf(geo, "Field[%i] = %s;\n", it.first, it.second->getName());
-    std::for_each(it.second->options.begin(), it.second->options.end(),
-                  writeFieldOptionGEO(geo, it.second));
-  }
+  private :
+    FILE *geo;
+  public :
+    writeFieldGEO(FILE *fp) { geo = fp ? fp : stdout; }
+    void operator() (std::pair<const int, Field *> it)
+    {
+      fprintf(geo, "Field[%i] = %s;\n", it.first, it.second->getName());
+      std::for_each(it.second->options.begin(), it.second->options.end(),
+          writeFieldOptionGEO(geo, it.second));
+    }
 };
 
 class writePhysicalGroupGEO {
- private :
-  FILE *geo;
-  int dim;
-  bool printLabels;
-  std::map<int, std::string> &oldLabels;
-  std::map<std::pair<int, int>, std::string> &newLabels;
- public :
-  writePhysicalGroupGEO(FILE *fp, int i, bool labels,
-                        std::map<int, std::string> &o,
-                        std::map<std::pair<int, int>, std::string> &n)
-    : dim(i), printLabels(labels), oldLabels(o), newLabels(n)
-  {
-    geo = fp ? fp : stdout;
-  }
-  void operator () (std::pair<const int, std::vector<GEntity *> > &g)
-  {
-    std::string oldName, newName;
-    if(printLabels){
-      if(newLabels.count(std::pair<int, int>(dim, g.first))) {
-        newName = newLabels[std::pair<int, int>(dim, g.first)];
-      }
-      else if(oldLabels.count(g.first)) {
-        oldName = oldLabels[g.first];
-        fprintf(geo, "%s = %d;\n", oldName.c_str(), g.first);
-      }
+  private :
+    FILE *geo;
+    int dim;
+    bool printLabels;
+    std::map<int, std::string> &oldLabels;
+    std::map<std::pair<int, int>, std::string> &newLabels;
+  public :
+    writePhysicalGroupGEO(FILE *fp, int i, bool labels,
+        std::map<int, std::string> &o,
+        std::map<std::pair<int, int>, std::string> &n)
+      : dim(i), printLabels(labels), oldLabels(o), newLabels(n)
+    {
+      geo = fp ? fp : stdout;
     }
+    void operator () (std::pair<const int, std::vector<GEntity *> > &g)
+    {
+      std::string oldName, newName;
+      if(printLabels){
+        if(newLabels.count(std::pair<int, int>(dim, g.first))) {
+          newName = newLabels[std::pair<int, int>(dim, g.first)];
+        }
+        else if(oldLabels.count(g.first)) {
+          oldName = oldLabels[g.first];
+          fprintf(geo, "%s = %d;\n", oldName.c_str(), g.first);
+        }
+      }
 
-    switch (dim) {
-    case 0: fprintf(geo, "Physical Point"); break;
-    case 1: fprintf(geo, "Physical Line"); break;
-    case 2: fprintf(geo, "Physical Surface"); break;
-    case 3: fprintf(geo, "Physical Volume"); break;
-    }
+      switch (dim) {
+        case 0: fprintf(geo, "Physical Point"); break;
+        case 1: fprintf(geo, "Physical Line"); break;
+        case 2: fprintf(geo, "Physical Surface"); break;
+        case 3: fprintf(geo, "Physical Volume"); break;
+      }
 
-    if(oldName.size())
-      fprintf(geo, "(%s) = {", oldName.c_str());
-    else if(newName.size())
-      fprintf(geo, "(\"%s\") = {", newName.c_str());
-    else
-      fprintf(geo, "(%d) = {", g.first);
-    for(unsigned int i = 0; i < g.second.size(); i++) {
-      if(i) fprintf(geo, ", ");
-      fprintf(geo, "%d", g.second[i]->tag());
+      if(oldName.size())
+        fprintf(geo, "(%s) = {", oldName.c_str());
+      else if(newName.size())
+        fprintf(geo, "(\"%s\") = {", newName.c_str());
+      else
+        fprintf(geo, "(%d) = {", g.first);
+      for(unsigned int i = 0; i < g.second.size(); i++) {
+        if(i) fprintf(geo, ", ");
+        fprintf(geo, "%d", g.second[i]->tag());
+      }
+      fprintf(geo, "};\n");
     }
-    fprintf(geo, "};\n");
-  }
 };
 
 static bool skipRegion(GRegion *gr)
@@ -527,7 +527,7 @@ int GModel::writeGEO(const std::string &name, bool printLabels, bool onlyPhysica
   getPhysicalGroups(groups);
   for(int i = 0; i < 4; i++)
     std::for_each(groups[i].begin(), groups[i].end(),
-                  writePhysicalGroupGEO(fp, i, printLabels, labels, physicalNames));
+        writePhysicalGroupGEO(fp, i, printLabels, labels, physicalNames));
 
   std::for_each(getFields()->begin(), getFields()->end(), writeFieldGEO(fp));
   if(getFields()->getBackgroundField() > 0)
diff --git a/Geo/GVertex.h b/Geo/GVertex.h
index 6210cadcba04e02c645c665546e9e44b0f0c748a..0326dd835993589b6a48ac7990b5dc33757ca145 100644
--- a/Geo/GVertex.h
+++ b/Geo/GVertex.h
@@ -64,7 +64,7 @@ class GVertex : public GEntity
   virtual GeomType geomType() const { return Point; }
 
   // get/set the prescribed mesh size at the vertex
-  inline double prescribedMeshSizeAtVertex() const { return meshSize; }
+  virtual inline double prescribedMeshSizeAtVertex() const { return meshSize; }
   virtual void setPrescribedMeshSizeAtVertex(double l) { meshSize = l; }
 
   // get the bounding box
diff --git a/Geo/GenericEdge.cpp b/Geo/GenericEdge.cpp
index 91b2cef330347a3b13f067dfc6ca65f4a0911dd9..c636eb83a35a7992590b8af754ec048feda11a5c 100644
--- a/Geo/GenericEdge.cpp
+++ b/Geo/GenericEdge.cpp
@@ -1,4 +1,4 @@
-// Gmsh - Copyright (C) 1997-2015 C. Geuzaine, J.-F. Remacle
+// Gmsh - Copyright (C) 1997-2014 C. Geuzaine, J.-F. Remacle
 //
 // See the LICENSE.txt file for license information. Please report all
 // bugs and problems to the public mailing list <gmsh@geuz.org>.
@@ -11,6 +11,8 @@
 #include "GenericFace.h"
 #include "Context.h"
 
+//------------------------------------------------------------------------
+
 GenericEdge::ptrfunction_int_double_refvector GenericEdge::EdgeEvalXYZFromT = NULL;
 GenericEdge::ptrfunction_int_refdouble_refdouble GenericEdge::EdgeEvalParBounds = NULL;
 GenericEdge::ptrfunction_int_refstring GenericEdge::EdgeGeomType = NULL;
@@ -21,34 +23,41 @@ GenericEdge::ptrfunction_int_refvector_refdouble_refvector_refbool GenericEdge::
 GenericEdge::ptrfunction_int_refbool GenericEdge::EdgeIs3D = NULL;
 GenericEdge::ptrfunction_int_int_double_int_refvector GenericEdge::EdgeReparamOnFace = NULL;
 
-GenericEdge::GenericEdge(GModel *m, int num, int _native_id, GVertex *v1, GVertex *v2)
-  : GEdge(m, num, v1, v2), id(_native_id)
-{
+//------------------------------------------------------------------------
+
+GenericEdge::GenericEdge(GModel *m, int num, int _native_id, GVertex *v1, GVertex *v2, bool _isseam):GEdge(m, num, v1, v2), id(_native_id),is_seam(_isseam){
   if ((!EdgeEvalParBounds)||(!EdgeEvalXYZFromT)) Msg::Error("GenericEdge::ERROR: Callback not set");
   bool ok = EdgeEvalParBounds(id,s0,s1);
   if (!ok) Msg::Error("GenericEdge::ERROR from EdgeEvalParBounds ! " );
 }
 
-GenericEdge::~GenericEdge()
-{
+//------------------------------------------------------------------------
+
+GenericEdge::~GenericEdge(){
 }
 
-Range<double> GenericEdge::parBounds(int i) const
-{
+//------------------------------------------------------------------------
+
+Range<double> GenericEdge::parBounds(int i) const{
   return Range<double>(s0, s1);
 }
 
-SPoint2 GenericEdge::reparamOnFace(const GFace *face, double par, int dir) const
-{
+//------------------------------------------------------------------------
+
+SPoint2 GenericEdge::reparamOnFace(const GFace *face, double par, int dir) const{
   vector<double> res(2,0.);
   if (!EdgeReparamOnFace) Msg::Error("GenericEdge::ERROR: Callback EdgeReparamOnFace not set");
   bool ok = EdgeReparamOnFace(id,face->getNativeInt(),par, dir,res);
-  if (!ok) Msg::Error("GenericEdge::ERROR from EdgeReparamOnFace ! " );
+  if (!ok){
+    cout << "GenericEdge::ERROR from EdgeReparamOnFace ! Edge Native id " << getNativeInt() << endl;
+    Msg::Error("GenericEdge::ERROR from EdgeReparamOnFace ! Edge Native id %d",getNativeInt() );
+  }
   return SPoint2(res[0],res[1]);;
 }
 
-GPoint GenericEdge::closestPoint(const SPoint3 &qp, double &param) const
-{
+//------------------------------------------------------------------------
+
+GPoint GenericEdge::closestPoint(const SPoint3 &qp, double &param) const{
   vector<double> queryPoint(3,0.);
   for (int i=0;i<3;i++) queryPoint[i] = qp[i];
   vector<double> res(3,0.);
@@ -60,15 +69,16 @@ GPoint GenericEdge::closestPoint(const SPoint3 &qp, double &param) const
   return GPoint(res[0], res[1], res[2], this, param);
 }
 
-// True if the edge is a seam for the given face
-bool GenericEdge::isSeam(const GFace *face) const
-{
-  Msg::Error("GenericEdge::isSeam: not implemented yet ! ");
-  return false;
+//------------------------------------------------------------------------
+
+bool GenericEdge::isSeam(const GFace *face) const{
+//  return false;
+  return is_seam;
 }
 
-GPoint GenericEdge::point(double par) const
-{
+//------------------------------------------------------------------------
+
+GPoint GenericEdge::point(double par) const{
   vector<double> res(3,0.);
   if (!EdgeEvalXYZFromT) Msg::Error("GenericEdge::ERROR: Callback EdgeEvalXYZFromT not set");
   bool ok = EdgeEvalXYZFromT(id,par,res);
@@ -76,8 +86,9 @@ GPoint GenericEdge::point(double par) const
   return GPoint(res[0], res[1], res[2], this, par);
 }
 
-SVector3 GenericEdge::firstDer(double par) const
-{
+//------------------------------------------------------------------------
+
+SVector3 GenericEdge::firstDer(double par) const{
   vector<double> res(3,0.);
   if (!EdgeEvalFirstDer) Msg::Error("GenericEdge::ERROR: Callback EdgeEvalFirstDer not set");
   bool ok = EdgeEvalFirstDer(id,par,res);
@@ -85,13 +96,14 @@ SVector3 GenericEdge::firstDer(double par) const
   return SVector3(res[0],res[1],res[2]);
 }
 
-GEntity::GeomType GenericEdge::geomType() const
-{
+//------------------------------------------------------------------------
+
+GEntity::GeomType GenericEdge::geomType() const{
   string s;
   if (!EdgeGeomType) Msg::Error("GenericEdge::ERROR: Callback EdgeGeomType not set");
   bool ok = EdgeGeomType(id,s);
   if (!ok){
-    Msg::Error("GenericEdge::ERROR from EdgeGeomType ! " );
+    Msg::Error("GenericEdge::ERROR from EdgeGeomType ! Edge Native id %d",getNativeInt() );
     return Unknown;
   }
 
@@ -111,14 +123,17 @@ GEntity::GeomType GenericEdge::geomType() const
     return BSpline;
   else if(s.compare("TrimmedCurve")==0)
     return TrimmedCurve;
+  else if(s.compare("Intersection curve")==0)
+    return BSpline;
 
-  Msg::Warning("GenericEdge::geomType:: unknown type from callback: ", s.c_str());
+  Msg::Warning("GenericEdge::geomType:: unknown type from callback: %s", s.c_str());
 
   return Unknown;
 }
 
-double GenericEdge::curvature(double par) const
-{
+//------------------------------------------------------------------------
+
+double GenericEdge::curvature(double par) const{
   double res;
   if (!EdgeEvalCurvature) Msg::Error("GenericEdge::ERROR: Callback EdgeEvalCurvature not set");
   bool ok = EdgeEvalCurvature(id,par,res);
@@ -126,8 +141,9 @@ double GenericEdge::curvature(double par) const
   return res;
 }
 
-bool GenericEdge::is3D() const
-{
+//------------------------------------------------------------------------
+
+bool GenericEdge::is3D() const{
   bool res;
   if (!EdgeIs3D) Msg::Error("GenericEdge::ERROR: Callback EdgeIs3D not set");
   bool ok = EdgeIs3D(id,res);
@@ -135,11 +151,99 @@ bool GenericEdge::is3D() const
   return res;
 }
 
-bool GenericEdge::degenerate(int) const
-{
-  bool res;
+//------------------------------------------------------------------------
+
+bool GenericEdge::degenerate(int) const{
+  bool res=false;
   if (!EdgeDegenerated) Msg::Error("GenericEdge::ERROR: Callback EdgeDegenerated not set");
   bool ok = EdgeDegenerated(id,res);
   if (!ok) Msg::Error("GenericEdge::ERROR from EdgeDegenerated ! " );
   return res;
 }
+
+//------------------------------------------------------------------------
+
+int GenericEdge::minimumDrawSegments() const
+{
+  if(geomType() == Line)
+    return GEdge::minimumDrawSegments();
+  else
+    return CTX::instance()->geom.numSubEdges * GEdge::minimumDrawSegments();
+}
+
+//------------------------------------------------------------------------
+
+
+
+
+
+//------------------------------------------------------------------------
+//------------------------------------------------------------------------
+//------------------------------------------------------------------------
+//------------------------------------------------------------------------
+//------------------------------------------------------------------------
+
+LinearSeamEdge::LinearSeamEdge(GModel *m, int num, GVertex *v1, GVertex *v2):GEdge(m, num, v1, v2){
+  s0=0.;
+  s1=v1->xyz().distance(v2->xyz());
+  first_der = SVector3(v1->xyz(),v2->xyz());
+  first_der.normalize();
+}
+
+//------------------------------------------------------------------------
+
+LinearSeamEdge::~LinearSeamEdge(){
+}
+
+//------------------------------------------------------------------------
+
+Range<double> LinearSeamEdge::parBounds(int i) const{
+  return Range<double>(s0, s1);
+}
+
+//------------------------------------------------------------------------
+
+GPoint LinearSeamEdge::point(double par) const{
+  SVector3 res = v0->xyz() + par*first_der;
+  return GPoint(res[0], res[1], res[2], this, par);
+}
+
+//------------------------------------------------------------------------
+
+SVector3 LinearSeamEdge::firstDer(double par) const{
+  return first_der;
+}
+
+//------------------------------------------------------------------------
+
+GEntity::GeomType LinearSeamEdge::geomType() const{
+  return Line;
+}
+
+//------------------------------------------------------------------------
+
+double LinearSeamEdge::curvature(double par) const{
+  return 0.;
+}
+
+//------------------------------------------------------------------------
+
+bool LinearSeamEdge::is3D() const{
+  return false;
+}
+
+//------------------------------------------------------------------------
+
+bool LinearSeamEdge::degenerate(int) const{
+  return false;
+}
+
+//------------------------------------------------------------------------
+
+GPoint LinearSeamEdge::closestPoint(const SPoint3 &q, double &t) const
+{
+  return GEdge::closestPoint(q,t);
+}
+
+//------------------------------------------------------------------------
+
diff --git a/Geo/GenericEdge.h b/Geo/GenericEdge.h
index 158e458084d9e019d3a0eb2c7b06fcaaf420a41d..e2ffd71b8ec95cffc7d81b7c8f9b6a2022067abf 100644
--- a/Geo/GenericEdge.h
+++ b/Geo/GenericEdge.h
@@ -1,4 +1,4 @@
-// Gmsh - Copyright (C) 1997-2015 C. Geuzaine, J.-F. Remacle
+// Gmsh - Copyright (C) 1997-2014 C. Geuzaine, J.-F. Remacle
 //
 // See the LICENSE.txt file for license information. Please report all
 // bugs and problems to the public mailing list <gmsh@geuz.org>.
@@ -18,67 +18,104 @@ using namespace std;
 class GenericFace;
 
 /*The set of Generic Entities is a generic interface to any other modeler.
-  Callbacks (function pointers) are given, sending requests, enquiries, to the
-  native modeler. */
+  Callbacks (function pointers) are given, sending requests, enquiries, to the native modeler. */
 
 class GenericEdge : public GEdge {
- protected:
-  double s0, s1;
-  int id;
- public:
-  // callbacks typedef
-  typedef bool (*ptrfunction_int_double_refvector)(int, double, vector<double>&);
-  typedef bool (*ptrfunction_int_refdouble_refdouble)(int, double&, double&);
-  typedef bool (*ptrfunction_int_double_refdouble)(int, double, double&);
-  typedef bool (*ptrfunction_int_refstring)(int, string&);
-  typedef bool (*ptrfunction_int_refbool)(int, bool&);
-  typedef bool (*ptrfunction_int_refvector_refdouble_refvector_refbool)
-    (int, const vector<double> &, double &, vector<double>&, bool &);
-  typedef bool (*ptrfunction_int_int_double_int_refvector)(int, int, double, int, vector<double> &);
-
-  GenericEdge(GModel *model, int num,int _native_id, GVertex *v1, GVertex *v2);
-  virtual ~GenericEdge();
-
-  virtual Range<double> parBounds(int i) const;
-  virtual GeomType geomType() const;
-  virtual bool degenerate(int) const;
-  virtual GPoint point(double p) const;
-  virtual SVector3 firstDer(double par) const;
-  virtual double curvature (double par) const;
-  virtual SPoint2 reparamOnFace(const GFace *face, double epar, int dir) const;
-  virtual GPoint closestPoint(const SPoint3 &queryPoint, double &param) const;
-
-  ModelType getNativeType() const { return GenericModel; }
-  virtual int getNativeInt()const{return id;};
-
-  bool is3D() const;
-  bool isSeam(const GFace *) const;
-
-  // sets the callbacks, to be given by the user
-  static void setEdgeEvalXYZFromT(ptrfunction_int_double_refvector fct){EdgeEvalXYZFromT = fct;};
-  static void setEdgeEvalParBounds(ptrfunction_int_refdouble_refdouble fct){EdgeEvalParBounds = fct;};
-  static void setEdgeGeomType(ptrfunction_int_refstring fct){EdgeGeomType = fct;};
-  static void setEdgeDegenerated(ptrfunction_int_refbool fct){EdgeDegenerated = fct;};
-  static void setEdgeEvalFirstDer(ptrfunction_int_double_refvector fct){EdgeEvalFirstDer = fct;};
-  static void setEdgeEvalCurvature(ptrfunction_int_double_refdouble fct){EdgeEvalCurvature = fct;};
-  static void setEdgeClosestPoint(ptrfunction_int_refvector_refdouble_refvector_refbool fct){EdgeClosestPoint = fct;};
-  static void setEdgeIs3D(ptrfunction_int_refbool fct){EdgeIs3D = fct;};
-  static void setEdgeReparamOnFace(ptrfunction_int_int_double_int_refvector fct){EdgeReparamOnFace = fct;};
-
- private:
-  // the callbacks:
-  static ptrfunction_int_double_refvector EdgeEvalXYZFromT;
-  static ptrfunction_int_refdouble_refdouble EdgeEvalParBounds;
-  static ptrfunction_int_refstring EdgeGeomType;
-  static ptrfunction_int_refbool EdgeDegenerated;
-  static ptrfunction_int_double_refvector EdgeEvalFirstDer;
-  static ptrfunction_int_double_refdouble EdgeEvalCurvature;
-  // the first vector is a query point xyz, fills the second vector with closest
-  // point on edge using orthogonal projection.  Fills double param with
-  // parametric coordinate of end point projection.
-  static ptrfunction_int_refvector_refdouble_refvector_refbool EdgeClosestPoint;
-  static ptrfunction_int_refbool EdgeIs3D;
-  static ptrfunction_int_int_double_int_refvector EdgeReparamOnFace;
+  protected:
+    double s0, s1;
+    int id;
+    const bool is_seam;
+  public:
+    // callbacks typedef
+    typedef bool (*ptrfunction_int_double_refvector)(int, double, vector<double>&);
+    typedef bool (*ptrfunction_int_refdouble_refdouble)(int, double&, double&);
+    typedef bool (*ptrfunction_int_double_refdouble)(int, double, double&);
+    typedef bool (*ptrfunction_int_refstring)(int, string&);
+    typedef bool (*ptrfunction_int_refbool)(int, bool&);
+    typedef bool (*ptrfunction_int_refvector_refdouble_refvector_refbool)(int, const vector<double> &, double &, vector<double>&, bool &);
+    typedef bool (*ptrfunction_int_int_double_int_refvector)(int, int, double, int, vector<double> &);
+
+    GenericEdge(GModel *model, int num,int _native_id, GVertex *v1, GVertex *v2, bool _isseam=false);
+    virtual ~GenericEdge();
+
+    virtual Range<double> parBounds(int i) const;
+    virtual GeomType geomType() const;
+    virtual bool degenerate(int) const;
+    virtual GPoint point(double p) const;
+    virtual SVector3 firstDer(double par) const;
+    virtual double curvature (double par) const;
+    virtual SPoint2 reparamOnFace(const GFace *face, double epar, int dir) const;
+    virtual GPoint closestPoint(const SPoint3 &queryPoint, double &param) const;
+
+    ModelType getNativeType() const { return GenericModel; }
+    virtual int getNativeInt()const{return id;};
+  
+    virtual int minimumDrawSegments () const;// for output
+
+    bool is3D() const;
+    bool isSeam(const GFace *) const;
+
+    // sets the callbacks, to be given by the user
+    static void setEdgeEvalXYZFromT(ptrfunction_int_double_refvector fct){EdgeEvalXYZFromT = fct;};
+    static void setEdgeEvalParBounds(ptrfunction_int_refdouble_refdouble fct){EdgeEvalParBounds = fct;};
+    static void setEdgeGeomType(ptrfunction_int_refstring fct){EdgeGeomType = fct;};
+    static void setEdgeDegenerated(ptrfunction_int_refbool fct){EdgeDegenerated = fct;};
+    static void setEdgeEvalFirstDer(ptrfunction_int_double_refvector fct){EdgeEvalFirstDer = fct;};
+    static void setEdgeEvalCurvature(ptrfunction_int_double_refdouble fct){EdgeEvalCurvature = fct;};
+    static void setEdgeClosestPoint(ptrfunction_int_refvector_refdouble_refvector_refbool fct){EdgeClosestPoint = fct;};
+    static void setEdgeIs3D(ptrfunction_int_refbool fct){EdgeIs3D = fct;};
+    static void setEdgeReparamOnFace(ptrfunction_int_int_double_int_refvector fct){EdgeReparamOnFace = fct;};
+
+  private:
+    // the callbacks:
+    // --------------
+    static ptrfunction_int_double_refvector EdgeEvalXYZFromT;
+    static ptrfunction_int_refdouble_refdouble EdgeEvalParBounds;
+    static ptrfunction_int_refstring EdgeGeomType;
+    static ptrfunction_int_refbool EdgeDegenerated;
+    static ptrfunction_int_double_refvector EdgeEvalFirstDer;
+    static ptrfunction_int_double_refdouble EdgeEvalCurvature;
+    // the first vector is a query point xyz, fills the second vector with closest point 
+    // on edge using orthogonal projection.  Fills double param with parametric coordinate of end point projection.
+    static ptrfunction_int_refvector_refdouble_refvector_refbool EdgeClosestPoint;
+    static ptrfunction_int_refbool EdgeIs3D;
+    static ptrfunction_int_int_double_int_refvector EdgeReparamOnFace;
+
+
 };
 
+
+
+
+
+
+
+
+
+
+
+class LinearSeamEdge : public GEdge {
+  protected:
+    double s0, s1;
+    SVector3 first_der;
+  public:
+    LinearSeamEdge(GModel *model, int num, GVertex *v1, GVertex *v2);
+    virtual ~LinearSeamEdge();
+
+    virtual Range<double> parBounds(int i) const;
+    virtual GeomType geomType() const;
+    virtual bool degenerate(int) const;
+    virtual GPoint point(double p) const;
+    virtual SVector3 firstDer(double par) const;
+    virtual double curvature (double par) const;
+    virtual bool isSeam(const GFace *face) const { return true; }
+    bool is3D() const;
+    virtual GPoint closestPoint(const SPoint3 &queryPoint, double &param) const;
+
+    ModelType getNativeType() const { return GenericModel; }
+
+
+};
+
+
 #endif
diff --git a/Geo/GenericFace.cpp b/Geo/GenericFace.cpp
index f2546ca2d181643612c23ad7a1d5dea4e944bff2..7baaf1211940d6c3775f10651fdf88cf08bc7854 100644
--- a/Geo/GenericFace.cpp
+++ b/Geo/GenericFace.cpp
@@ -1,4 +1,4 @@
-// Gmsh - Copyright (C) 1997-2015 C. Geuzaine, J.-F. Remacle
+// Gmsh - Copyright (C) 1997-2014 C. Geuzaine, J.-F. Remacle
 //
 // See the LICENSE.txt file for license information. Please report all
 // bugs and problems to the public mailing list <gmsh@geuz.org>.
@@ -14,6 +14,9 @@
 #include "Context.h"
 #include <math.h>
 
+
+//------------------------------------------------------------------------
+
 GenericFace::ptrfunction_int_refstring GenericFace::FaceGeomType = NULL;
 GenericFace::ptrfunction_int_refvector_refvector GenericFace::FaceUVFromXYZ = NULL;
 GenericFace::ptrfunction_int_refvector_refvector_refvector GenericFace::FaceClosestPoint = NULL;
@@ -25,63 +28,66 @@ GenericFace::ptrfunction_int_refvector_refvector_refvector_refdouble_refdouble G
 GenericFace::ptrfunction_int_refvector_refvector GenericFace::FaceEvalNormal = NULL;
 GenericFace::ptrfunction_int_refvector_refvector_refvector GenericFace::FaceFirstDer = NULL;
 GenericFace::ptrfunction_int_refvector_refvector_refvector_refvector GenericFace::FaceSecondDer = NULL;
+GenericFace::ptrfunction_int_refbool_refbool_refdouble_refdouble GenericFace::FacePeriodicInfo = NULL;
+
+//------------------------------------------------------------------------
 
-GenericFace::GenericFace(GModel *m, int num, int _native_id):GFace(m, num), id(_native_id)
-{
+GenericFace::GenericFace(GModel *m, int num, int _native_id):GFace(m, num), id(_native_id){
   if (!FaceParBounds) Msg::Fatal("Genericface::ERROR: Callback FaceParBounds not set");
-  Range<double> rangeu = parBounds(0);
-  Range<double> rangev = parBounds(1);
-  umin = rangeu.low();
-  umax = rangeu.high();
-  vmin = rangev.low();
-  vmax = rangev.high();
-
-  // TODO: set periodic or not !!!
-  //bool _periodic[2];// is periodic in u, v
-  //  throw;
+
+  bool ok = FaceParBounds(id,0,umin,umax);
+  if (!ok) Msg::Error("GenericEdge::ERROR from EdgeEvalParBounds ! " );
+  ok = FaceParBounds(id,1,vmin,vmax);
+    
+  computePeriodicity();
 }
 
-GenericFace::~GenericFace()
-{
+//------------------------------------------------------------------------
+
+GenericFace::~GenericFace(){
 }
 
-Range<double> GenericFace::parBounds(int i) const
-{
+//------------------------------------------------------------------------
+
+Range<double> GenericFace::parBounds(int i) const{
   if(i == 0) return Range<double>(umin, umax);
   return Range<double>(vmin, vmax);
 }
 
-SVector3 GenericFace::normal(const SPoint2 &param) const
-{
+//------------------------------------------------------------------------
+
+SVector3 GenericFace::normal(const SPoint2 &param) const{
   vector<double> res(3,0.);
   vector<double> par(2,0.);
-  for (int i=0;i<2;i++) par[i] = param[i];
+  for (int i=0;i<3;i++) par[i] = param[i];
   if (!FaceEvalNormal) Msg::Fatal("Genericface::ERROR: Callback FaceEvalNormal not set");
   bool ok = FaceEvalNormal(id,par,res);
   if (!ok) Msg::Error("GenericFace::ERROR from FaceEvalNormal ! " );
   return SVector3(res[0],res[1],res[2]);
 }
 
-Pair<SVector3,SVector3> GenericFace::firstDer(const SPoint2 &param) const
-{
+//------------------------------------------------------------------------
+
+Pair<SVector3,SVector3> GenericFace::firstDer(const SPoint2 &param) const{
   if (!FaceFirstDer) Msg::Fatal("Genericface::ERROR: Callback FaceFirstDer not set");
   vector<double> deru(3,0.);
   vector<double> derv(3,0.);
   vector<double> par(2,0.);
-  for (int i=0;i<2;i++) par[i] = param[i];
+  for (int i=0;i<3;i++) par[i] = param[i];
   bool ok = FaceFirstDer(id,par,deru,derv);
   if (!ok) Msg::Error("GenericFace::ERROR from FaceFirstDer ! " );
   return Pair<SVector3,SVector3>(SVector3(deru[0],deru[1],deru[2]),
                                  SVector3(derv[0],derv[1],derv[2]));
 }
 
-void GenericFace::secondDer(const SPoint2 &param,SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const
-{
+//------------------------------------------------------------------------
+
+void GenericFace::secondDer(const SPoint2 &param,SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const{
   vector<double> deruu(3,0.);
   vector<double> dervv(3,0.);
   vector<double> deruv(3,0.);
   vector<double> par(2,0.);
-  for (int i=0;i<2;i++) par[i] = param[i];
+  for (int i=0;i<3;i++) par[i] = param[i];
   if (!FaceSecondDer) Msg::Fatal("Genericface::ERROR: Callback FaceSecondDer not set");
   bool ok = FaceSecondDer(id,par,deruu,dervv,deruv);
   if (!ok) Msg::Error("GenericFace::ERROR from FaceSecondDer ! " );
@@ -91,21 +97,28 @@ void GenericFace::secondDer(const SPoint2 &param,SVector3 *dudu, SVector3 *dvdv,
   return;
 }
 
-GPoint GenericFace::point(double par1, double par2) const
-{
+//------------------------------------------------------------------------
+
+GPoint GenericFace::point(double par1, double par2) const{
   vector<double> uv(2,0.);
   uv[0] = par1;
   uv[1] = par2;
   vector<double> xyz(3,0.);
+  double pp[2] = {par1, par2};
   if (!FaceXYZFromUV) Msg::Fatal("Genericface::ERROR: Callback FaceXYZFromUV not set");
   bool ok = FaceXYZFromUV(id,uv,xyz);
-  if (!ok) Msg::Error("GenericFace::ERROR from FaceXYZFromUV ! " );
-  double pp[2] = {par1, par2};
+  if (!ok){
+    Msg::Error("GenericFace::ERROR from FaceXYZFromUV ! " );
+    GPoint p(xyz[0], xyz[1], xyz[2], this, pp);
+    p.setNoSuccess();
+    return p;
+  }
   return GPoint(xyz[0], xyz[1], xyz[2], this, pp);
 }
 
-GPoint GenericFace::closestPoint(const SPoint3 &qp, const double initialGuess[2]) const
-{
+//------------------------------------------------------------------------
+
+GPoint GenericFace::closestPoint(const SPoint3 &qp, const double initialGuess[2]) const{
   vector<double> uvres(2,0.);
   vector<double> xyzres(3,0.);
   vector<double> queryPoint(3,0.);
@@ -129,8 +142,9 @@ GPoint GenericFace::closestPoint(const SPoint3 &qp, const double initialGuess[2]
   return GPoint(xyzres[0], xyzres[1], xyzres[2], this, pp);
 }
 
-SPoint2 GenericFace::parFromPoint(const SPoint3 &qp, bool onSurface) const
-{
+//------------------------------------------------------------------------
+
+SPoint2 GenericFace::parFromPoint(const SPoint3 &qp, bool onSurface) const{
   vector<double> uvres(2,0.);
   vector<double> xyzres(3,0.);
   vector<double> queryPoint(3,0.);
@@ -138,7 +152,7 @@ SPoint2 GenericFace::parFromPoint(const SPoint3 &qp, bool onSurface) const
   bool ok=true;
   if (onSurface){
     if (!FaceUVFromXYZ) Msg::Fatal("Genericface::ERROR: Callback FaceUVFromXYZ not set");
-    ok = FaceUVFromXYZ(id,queryPoint,uvres);// assuming point is on surface
+    ok = FaceUVFromXYZ(id,queryPoint,uvres);// assuming point is on surface 
     if (!ok) Msg::Error("GenericFace::ERROR from FaceUVFromXYZ ! " );
   }
   if ((!onSurface)||(!ok)){// if not on surface
@@ -148,14 +162,15 @@ SPoint2 GenericFace::parFromPoint(const SPoint3 &qp, bool onSurface) const
     bool on_the_face;
     if (!FaceContainsPointFromXYZ) Msg::Fatal("Genericface::ERROR: Callback FaceContainsPointFromXYZ not set");
     ok = FaceContainsPointFromXYZ(id,xyzres,on_the_face);// check if the projected point lies on the face...
-    if (!ok) Msg::Error("GenericFace::ERROR from FaceContainsPointFromXYZ ! " );
+    if (!ok) Msg::Error("GenericFace::parFromPoint::ERROR from FaceContainsPointFromXYZ ! " );
     if (!on_the_face) Msg::Warning("GenericFace::parFromPoint::Warning !!!! The returned point does not lies on the face ! " );
   }
   return SPoint2(uvres[0],uvres[1]);
 }
 
-GEntity::GeomType GenericFace::geomType() const
-{
+//------------------------------------------------------------------------
+
+GEntity::GeomType GenericFace::geomType() const{
   string s;
   if (!FaceGeomType) Msg::Fatal("Genericface::ERROR: Callback FaceGeomType not set");
   bool ok = FaceGeomType(id,s);
@@ -184,24 +199,25 @@ GEntity::GeomType GenericFace::geomType() const
   return Unknown;
 }
 
-double GenericFace::curvatureMax(const SPoint2 &param) const
-{
+//------------------------------------------------------------------------
+
+double GenericFace::curvatureMax(const SPoint2 &param) const{
   vector<double> dirMax(3,0.);
   vector<double> dirMin(3,0.);
   double curvMax,curvMin;
   vector<double> par(2,0.);
-  for (int i=0;i<2;i++) par[i] = param[i];
+  for (int i=0;i<3;i++) par[i] = param[i];
   if (!FaceCurvatures) Msg::Fatal("Genericface::ERROR: Callback FaceCurvatures not set");
   bool ok = FaceCurvatures(id,par,dirMax,dirMin,curvMax,curvMin);
   if (!ok) Msg::Error("GenericFace::ERROR from FaceCurvatures ! " );
   return std::max(fabs(curvMax), fabs(curvMin));
 }
 
-double GenericFace::curvatures(const SPoint2 &_param,SVector3 *_dirMax,SVector3 *_dirMin,
-                               double *curvMax,double *curvMin) const
-{
+//------------------------------------------------------------------------
+
+double GenericFace::curvatures(const SPoint2 &_param,SVector3 *_dirMax,SVector3 *_dirMin,double *curvMax,double *curvMin) const{
   vector<double> param(2,0.);
-  for (int i=0;i<2;i++) param[i] = _param[i];
+  for (int i=0;i<3;i++) param[i] = _param[i];
   vector<double> dirMax(3,0.);
   vector<double> dirMin(3,0.);
 
@@ -214,14 +230,74 @@ double GenericFace::curvatures(const SPoint2 &_param,SVector3 *_dirMax,SVector3
   return *curvMax;
 }
 
-bool GenericFace::containsPoint(const SPoint3 &pt) const
-{
+//------------------------------------------------------------------------
+
+bool GenericFace::containsPoint(const SPoint3 &pt) const{
   bool res;
   vector<double> queryPoint(3,0.);
   for (int i=0;i<3;i++) queryPoint[i] = pt[i];
   if (!FaceContainsPointFromXYZ) Msg::Fatal("Genericface::ERROR: Callback FaceContainsPointFromXYZ not set");
   bool ok = FaceContainsPointFromXYZ(id,queryPoint,res);
-  if (!ok) Msg::Error("GenericFace::ERROR from FaceContainsPointFromXYZ ! " );
+  if (!ok) Msg::Error("GenericFace::containsPoint::ERROR from FaceContainsPointFromXYZ ! " );
   return res;
 }
 
+//------------------------------------------------------------------------
+
+void GenericFace::computePeriodicity(){
+  if (!FacePeriodicInfo) Msg::Fatal("Genericface::ERROR: Callback FacePeriodicInfo not set");
+  bool ok = FacePeriodicInfo(id,_periodic[0], _periodic[1],_period[0],_period[1]);
+  if (!ok) Msg::Error("GenericFace::ERROR from FacePeriodicInfo ! " );
+
+
+}
+
+//------------------------------------------------------------------------
+
+void GenericFace::createLoops(){
+  const bool debug = false;
+
+  if (debug){
+    cout << "Initial loops:";
+    for (std::list<GEdge *>::iterator it = l_edges.begin();it!=l_edges.end();it++){
+      cout << (*it)->tag() << " ";
+    }
+    cout << endl;
+  }
+
+
+  edgeLoops.clear();
+  l_edges.clear();
+  l_dirs.clear();
+  if (debug) cout << "Creating loops, face " << tag() << " (" << getNativeInt() << ")" << endl;
+
+  for (set<int>::iterator it_loop = loopsnumber.begin();it_loop!=loopsnumber.end();it_loop++){// for each loop
+    if (debug) cout << "Loop #" << *it_loop << " made of edges ";
+    pair<multimap<int, pair<GEdge*,int> >::iterator, multimap<int, pair<GEdge*,int> >::iterator> range = bnd.equal_range(*it_loop);
+    list<GEdge*> l_wire;
+    for (multimap<int, pair<GEdge*,int> >::iterator it = range.first;it!=range.second;it++){// for all edges
+      if (debug) cout << it->second.first->tag() << " ";
+      l_wire.push_back(it->second.first);      
+    }
+    if (debug) cout << endl << "-> GEdgeLoop made of edges (sign) ";
+    GEdgeLoop el(l_wire);
+    for(GEdgeLoop::citer it = el.begin(); it != el.end(); ++it){
+      l_edges.push_back(it->ge);
+      l_dirs.push_back(it->_sign);
+      if (debug) cout << it->ge->tag() << "(" << ((it->_sign<0.) ? "-" : "+") << ") " ;
+      if (el.count() == 2){
+        it->ge->meshAttributes.minimumMeshSegments = 
+          std::max(it->ge->meshAttributes.minimumMeshSegments,2);
+      }
+      if (el.count() == 1){
+        it->ge->meshAttributes.minimumMeshSegments = 
+          std::max(it->ge->meshAttributes.minimumMeshSegments,3);
+      }
+    }
+    if (debug) cout << endl;
+    edgeLoops.push_back(el);   
+  }
+}
+
+//------------------------------------------------------------------------
+
diff --git a/Geo/GenericFace.h b/Geo/GenericFace.h
index dbacfa021a7ccf8d54d56597642fb4c0b8ea95ef..395b982f208fcc6ca2c6c4019af69651d78f4dba 100644
--- a/Geo/GenericFace.h
+++ b/Geo/GenericFace.h
@@ -1,4 +1,4 @@
-// Gmsh - Copyright (C) 1997-2015 C. Geuzaine, J.-F. Remacle
+// Gmsh - Copyright (C) 1997-2014 C. Geuzaine, J.-F. Remacle
 //
 // See the LICENSE.txt file for license information. Please report all
 // bugs and problems to the public mailing list <gmsh@geuz.org>.
@@ -16,80 +16,88 @@
 
 using namespace std;
 
-/* The set of Generic Entities is a generic interface to any other modeler.
-   Callbacks (function pointers) are given, sending requests, enquiries, to the
-   native modeler. */
+/*The set of Generic Entities is a generic interface to any other modeler.
+  Callbacks (function pointers) are given, sending requests, enquiries, to the native modeler. */
 
 class GenericFace : public GFace {
- protected:
-  int id;
-  double umin, umax, vmin, vmax;// face uv bounds
-  bool _periodic[2];// is periodic in u, v
+  protected:
+    int id;
+    double umin, umax, vmin, vmax;// face uv bounds
+    bool _periodic[2];// is periodic in u, v
+    double _period[2];
 
   public:
-  // callbacks typedef
-  typedef bool (*ptrfunction_int_refstring)(int, string&);
-  typedef bool (*ptrfunction_int_refvector_refvector)(const int , const vector<double> &, vector<double> &);
-  typedef bool (*ptrfunction_int_refvector_refvector_refvector)
-               (const int , const vector<double> &, vector<double> &, vector<double> &);
-  typedef bool (*ptrfunction_int_refvector_refbool)(const int , const vector<double> &, bool &);
-  typedef bool (*ptrfunction_int_int_refdouble_refdouble)(const int, const int, double &, double &);
-  typedef bool (*ptrfunction_int_refvector_refvector_refvector_refdouble_refdouble)
-               (const int, const vector<double> &, vector<double> &, vector<double> &, double &, double &);
-  typedef bool (*ptrfunction_int_refvector_refvector_refvector_refvector)
-               (const int , const vector<double> &, vector<double> &, vector<double> &, vector<double> &);
-
-  GenericFace(GModel *m, int num, int _native_id);
-  virtual ~GenericFace();
-
-  Range<double> parBounds(int i) const;
-  virtual GPoint point(double par1, double par2) const;
-  virtual GPoint closestPoint(const SPoint3 & queryPoint,
-                              const double initialGuess[2]) const;
-  virtual bool containsPoint(const SPoint3 &pt) const;
-  virtual SVector3 normal(const SPoint2 &param) const;
-  virtual Pair<SVector3,SVector3> firstDer(const SPoint2 &param) const;
-  virtual void secondDer(const SPoint2 &, SVector3 *, SVector3 *, SVector3 *) const;
-  virtual GEntity::GeomType geomType() const;
-  virtual SPoint2 parFromPoint(const SPoint3 &, bool onSurface=true) const;
-  virtual double curvatureMax(const SPoint2 &param) const;
-  virtual double curvatures(const SPoint2 &param, SVector3 *dirMax, SVector3 *dirMin,
-                            double *curvMax, double *curvMin) const;
-
-  ModelType getNativeType() const { return GenericModel; }
-  virtual int getNativeInt()const{return id;};
-
-  void addBndInfo(int loop, GEdge *ptr, int sign){
-    bnd.insert(make_pair(loop, make_pair(ptr,sign)));
-    l_dirs.push_back(sign);
-    l_edges.push_back(ptr);
-    ptr->addFace(this);
-  };
-
-  // sets callbacks
-  static void setFaceGeomType(ptrfunction_int_refstring fct){FaceGeomType = fct;};
-  static void setFaceUVFromXYZ(ptrfunction_int_refvector_refvector fct){FaceUVFromXYZ = fct;};
-  static void setFaceClosestPoint(ptrfunction_int_refvector_refvector_refvector fct){FaceClosestPoint = fct;};
-  static void setFaceContainsPointFromXYZ(ptrfunction_int_refvector_refbool fct){FaceContainsPointFromXYZ = fct;};
-  static void setFaceContainsPointFromUV(ptrfunction_int_refvector_refbool fct){FaceContainsPointFromUV = fct;};
-  static void setFaceXYZFromUV(ptrfunction_int_refvector_refvector fct){FaceXYZFromUV = fct;};
-  static void setFaceParBounds(ptrfunction_int_int_refdouble_refdouble fct){FaceParBounds = fct;};
-  static void setFaceCurvatures(ptrfunction_int_refvector_refvector_refvector_refdouble_refdouble fct){FaceCurvatures = fct;};
-  static void setFaceEvalNormal(ptrfunction_int_refvector_refvector fct){FaceEvalNormal = fct;};
-  static void setFaceFirstDer(ptrfunction_int_refvector_refvector_refvector fct){FaceFirstDer = fct;};
-  static void setFaceSecondDer(ptrfunction_int_refvector_refvector_refvector_refvector fct){FaceSecondDer = fct;};
-
-private:
-  multimap<int, pair<GEdge*,int> > bnd;
-
-  // the callbacks:
-  static ptrfunction_int_refstring FaceGeomType;
-  static ptrfunction_int_refvector_refvector FaceUVFromXYZ,FaceXYZFromUV,FaceEvalNormal;
-  static ptrfunction_int_refvector_refvector_refvector FaceClosestPoint,FaceFirstDer;
-  static ptrfunction_int_refvector_refbool FaceContainsPointFromXYZ,FaceContainsPointFromUV;
-  static ptrfunction_int_int_refdouble_refdouble FaceParBounds;
-  static ptrfunction_int_refvector_refvector_refvector_refdouble_refdouble FaceCurvatures;
-  static ptrfunction_int_refvector_refvector_refvector_refvector FaceSecondDer;
+    // callbacks typedef
+    typedef bool (*ptrfunction_int_refstring)(int, string&);
+    typedef bool (*ptrfunction_int_refbool_refbool_refdouble_refdouble)(int, bool&, bool&, double&, double&);
+    typedef bool (*ptrfunction_int_refvector_refvector)(const int , const vector<double> &, vector<double> &);
+    typedef bool (*ptrfunction_int_refvector_refvector_refvector)(const int , const vector<double> &, vector<double> &, vector<double> &);
+    typedef bool (*ptrfunction_int_refvector_refbool)(const int , const vector<double> &, bool &);
+    typedef bool (*ptrfunction_int_int_refdouble_refdouble)(const int, const int, double &, double &);
+    typedef bool (*ptrfunction_int_refvector_refvector_refvector_refdouble_refdouble)(const int, const vector<double> &, vector<double> &, vector<double> &, double &, double &);
+    typedef bool (*ptrfunction_int_refvector_refvector_refvector_refvector)(const int , const vector<double> &, vector<double> &, vector<double> &, vector<double> &);
+
+    GenericFace(GModel *m, int num, int _native_id);
+    virtual ~GenericFace();
+
+    Range<double> parBounds(int i) const;
+    virtual GPoint point(double par1, double par2) const;
+    virtual GPoint closestPoint(const SPoint3 & queryPoint,
+        const double initialGuess[2]) const;
+    virtual bool containsPoint(const SPoint3 &pt) const;
+    virtual SVector3 normal(const SPoint2 &param) const;
+    virtual Pair<SVector3,SVector3> firstDer(const SPoint2 &param) const;
+    virtual void secondDer(const SPoint2 &, SVector3 *, SVector3 *, SVector3 *) const;
+    virtual GEntity::GeomType geomType() const;
+    virtual SPoint2 parFromPoint(const SPoint3 &, bool onSurface=true) const;
+    virtual double curvatureMax(const SPoint2 &param) const;
+    virtual double curvatures(const SPoint2 &param, SVector3 *dirMax, SVector3 *dirMin,
+        double *curvMax, double *curvMin) const;
+    virtual bool periodic(int dir) const { return false;}//_periodic[dir]; }// TODO ?
+    virtual double period(int dir) const{ return _period[dir]; };
+
+    ModelType getNativeType() const { return GenericModel; }
+    virtual int getNativeInt()const{return id;};
+
+    void addBndInfo(int loop, GEdge *ptr, int sign){
+      bnd.insert(make_pair(loop, make_pair(ptr,sign)));
+      l_dirs.push_back(sign);
+      l_edges.push_back(ptr);
+      loopsnumber.insert(loop);
+      ptr->addFace(this);
+    };
+    void createLoops();
+
+    void computePeriodicity();
+
+    // sets callbacks
+    static void setFaceGeomType(ptrfunction_int_refstring fct){FaceGeomType = fct;};
+    static void setFaceUVFromXYZ(ptrfunction_int_refvector_refvector fct){FaceUVFromXYZ = fct;};
+    static void setFaceClosestPoint(ptrfunction_int_refvector_refvector_refvector fct){FaceClosestPoint = fct;};
+    static void setFaceContainsPointFromXYZ(ptrfunction_int_refvector_refbool fct){FaceContainsPointFromXYZ = fct;};
+    static void setFaceContainsPointFromUV(ptrfunction_int_refvector_refbool fct){FaceContainsPointFromUV = fct;};
+    static void setFaceXYZFromUV(ptrfunction_int_refvector_refvector fct){FaceXYZFromUV = fct;};
+    static void setFaceParBounds(ptrfunction_int_int_refdouble_refdouble fct){FaceParBounds = fct;};
+    static void setFaceCurvatures(ptrfunction_int_refvector_refvector_refvector_refdouble_refdouble fct){FaceCurvatures = fct;};
+    static void setFaceEvalNormal(ptrfunction_int_refvector_refvector fct){FaceEvalNormal = fct;};
+    static void setFaceFirstDer(ptrfunction_int_refvector_refvector_refvector fct){FaceFirstDer = fct;};
+    static void setFaceSecondDer(ptrfunction_int_refvector_refvector_refvector_refvector fct){FaceSecondDer = fct;};
+    static void setFacePeriodicInfo(ptrfunction_int_refbool_refbool_refdouble_refdouble fct){FacePeriodicInfo = fct;};
+
+  private:
+    multimap<int, pair<GEdge*,int> > bnd;
+    set<int> loopsnumber;
+
+    // the callbacks:
+    static ptrfunction_int_refstring FaceGeomType;
+    static ptrfunction_int_refvector_refvector FaceUVFromXYZ,FaceXYZFromUV,FaceEvalNormal;
+    static ptrfunction_int_refvector_refvector_refvector FaceClosestPoint,FaceFirstDer;
+    static ptrfunction_int_refvector_refbool FaceContainsPointFromXYZ,FaceContainsPointFromUV;
+    static ptrfunction_int_int_refdouble_refdouble FaceParBounds;
+    static ptrfunction_int_refvector_refvector_refvector_refdouble_refdouble FaceCurvatures;
+    static ptrfunction_int_refvector_refvector_refvector_refvector FaceSecondDer;
+    static ptrfunction_int_refbool_refbool_refdouble_refdouble FacePeriodicInfo;
+
 };
 
 #endif
diff --git a/Geo/GenericRegion.cpp b/Geo/GenericRegion.cpp
index 8f965649300581747d298d63ff3b0a97f0e2f17f..a6b8f1a408777c553d70eddad8ebeb541a48d391 100644
--- a/Geo/GenericRegion.cpp
+++ b/Geo/GenericRegion.cpp
@@ -1,4 +1,4 @@
-// Gmsh - Copyright (C) 1997-2015 C. Geuzaine, J.-F. Remacle
+// Gmsh - Copyright (C) 1997-2014 C. Geuzaine, J.-F. Remacle
 //
 // See the LICENSE.txt file for license information. Please report all
 // bugs and problems to the public mailing list <gmsh@geuz.org>.
@@ -11,15 +11,24 @@
 #include "GenericFace.h"
 #include "GenericRegion.h"
 
+//------------------------------------------------------------------------
+
 GenericRegion::GenericRegion(GModel *m, int num, int _native_id):GRegion(m, num), id(_native_id)
 {
 }
 
+//------------------------------------------------------------------------
+
 GenericRegion::~GenericRegion()
 {
 }
 
+//------------------------------------------------------------------------
+
 GEntity::GeomType GenericRegion::geomType() const
 {
   return Unknown;
 }
+
+//------------------------------------------------------------------------
+
diff --git a/Geo/GenericRegion.h b/Geo/GenericRegion.h
index 8c00eb862c15beefcda1ac3cc5876c753f4a0bf2..74850cc8b2e7cb3535f48a167b283391adc97e48 100644
--- a/Geo/GenericRegion.h
+++ b/Geo/GenericRegion.h
@@ -1,4 +1,4 @@
-// Gmsh - Copyright (C) 1997-2015 C. Geuzaine, J.-F. Remacle
+// Gmsh - Copyright (C) 1997-2014 C. Geuzaine, J.-F. Remacle
 //
 // See the LICENSE.txt file for license information. Please report all
 // bugs and problems to the public mailing list <gmsh@geuz.org>.
@@ -10,30 +10,28 @@
 #include "GRegion.h"
 
 /*The set of Generic Entities is a generic interface to any other modeler.
-  Callbacks (function pointers) are given, sending requests, enquiries, to the
-  native modeler. */
+  Callbacks (function pointers) are given, sending requests, enquiries, to the native modeler. */
 
 class GenericRegion : public GRegion {
- public:
-  GenericRegion(GModel *m, int num, int _native_id);
-  virtual ~GenericRegion();
-
-  virtual GeomType geomType() const;
-
-  ModelType getNativeType() const { return GenericModel; }
-  virtual int getNativeInt()const{return id;};
-
-  // TODO: When using GRegion->l_dirs and l_faces, what is the convention for
-  // l_dirs ? For now, assuming positive value for normals pointing inside the
-  // region.
-  void addFace(GenericFace *ptr, int sign){
-    l_dirs.push_back(sign);
-    l_faces.push_back(ptr);
-    ptr->addRegion(this);
-  };
-
- private:
-  int id;
+  public:
+    GenericRegion(GModel *m, int num, int _native_id);
+    virtual ~GenericRegion();
+
+    virtual GeomType geomType() const;
+
+    ModelType getNativeType() const { return GenericModel; }
+    virtual int getNativeInt()const{return id;};
+
+    // TODO: When using GRegion->l_dirs and l_faces, what is the convention for l_dirs ? For now, assuming positive value for normals pointing inside the region.
+    void addFace(GenericFace *ptr, int sign){
+      l_dirs.push_back(sign);
+      l_faces.push_back(ptr);
+      ptr->addRegion(this);
+    };
+
+
+  private:
+    int id;
 };
 
 #endif
diff --git a/Geo/GenericVertex.cpp b/Geo/GenericVertex.cpp
index 019875437e9324301d0308205494b9959f794905..819f26a98a8a4ae096aa48badf9fe4201b658866 100644
--- a/Geo/GenericVertex.cpp
+++ b/Geo/GenericVertex.cpp
@@ -1,4 +1,4 @@
-// Gmsh - Copyright (C) 1997-2015 C. Geuzaine, J.-F. Remacle
+// Gmsh - Copyright (C) 1997-2014 C. Geuzaine, J.-F. Remacle
 //
 // See the LICENSE.txt file for license information. Please report all
 // bugs and problems to the public mailing list <gmsh@geuz.org>.
@@ -12,10 +12,14 @@
 #include "GenericFace.h"
 #include<algorithm>
 
+//------------------------------------------------------------------------
+
 GenericVertex::ptrfunction_int_vector GenericVertex::VertexXYZ = NULL;
+GenericVertex::ptrfunction_int_doubleptr_voidptr GenericVertex::VertexMeshSize = NULL;
 
-GenericVertex::GenericVertex(GModel *m, int num, int _native_id):GVertex(m, num), id(_native_id)
-{
+//------------------------------------------------------------------------
+
+GenericVertex::GenericVertex(GModel *m, int num, int _native_id):GVertex(m, num), id(_native_id){
   if (!VertexXYZ)
     Msg::Fatal("GenericVertex::ERROR: Callback not set");
 
@@ -27,48 +31,34 @@ GenericVertex::GenericVertex(GModel *m, int num, int _native_id):GVertex(m, num)
   _z=vec[2];
 }
 
-GenericVertex::~GenericVertex()
-{
+//------------------------------------------------------------------------
+
+GenericVertex::GenericVertex(GModel *m, int num, int _native_id, const vector<double> &vec):GVertex(m, num), id(_native_id){
+  if (!VertexXYZ)
+    Msg::Fatal("GenericVertex::ERROR: Callback not set");
+
+  _x=vec[0];
+  _y=vec[1];
+  _z=vec[2];
 }
 
+//------------------------------------------------------------------------
+
+GenericVertex::~GenericVertex(){
+}
+
+//------------------------------------------------------------------------
+
 SPoint2 GenericVertex::reparamOnFace(const GFace *gf, int dir) const
 {
-  std::list<GEdge*>::const_iterator it = l_edges.begin();
-  //      // TODO: isSeam not ready yet !!!
-  //  while(it != l_edges.end()){
-  //    std::list<GEdge*> l = gf->edges();
-  //    if(std::find(l.begin(), l.end(), *it) != l.end()){
-  //      if((*it)->isSeam(gf)){
-  //        const TopoDS_Face *s = (TopoDS_Face*)gf->getNativePtr();
-  //        const TopoDS_Edge *c = (TopoDS_Edge*)(*it)->getNativePtr();
-  //        double s1,s0;
-  //        Handle(Geom2d_Curve) curve2d = BRep_Tool::CurveOnSurface(*c, *s, s0, s1);
-  //        if((*it)->getBeginVertex() == this)
-  //          return (*it)->reparamOnFace(gf, s0, dir);
-  //        else if((*it)->getEndVertex() == this)
-  //          return (*it)->reparamOnFace(gf, s1, dir);
-  //      }
-  //    }
-  //    ++it;
-  //  }
-  it = l_edges.begin();
-  while(it != l_edges.end()){
-    std::list<GEdge*> l = gf->edges();
-    if(std::find(l.begin(), l.end(), *it) != l.end()){
-      // the vertex belongs to an edge, which belongs to a face -> just find if begin or end vertex and edge->reparamOnFace(begin or end param)
-      Range<double> vec = (*it)->parBounds(0);
-      if((*it)->getBeginVertex() == this)
-        return (*it)->reparamOnFace(gf, vec.low(), dir);
-      else if((*it)->getEndVertex() == this)
-        return (*it)->reparamOnFace(gf, vec.high(), dir);
-    }
-    ++it;
-  }
 
-  // normally never here
-  return GVertex::reparamOnFace(gf, dir);
+  SPoint3 pt(_x,_y,_z);
+  return gf->parFromPoint(pt,true);
+
 }
 
+//------------------------------------------------------------------------
+
 void GenericVertex::setPosition(GPoint &p)
 {
   _x = p.x();
@@ -80,3 +70,5 @@ void GenericVertex::setPosition(GPoint &p)
     mesh_vertices[0]->z() = p.z();
   }
 }
+
+//------------------------------------------------------------------------
diff --git a/Geo/GenericVertex.h b/Geo/GenericVertex.h
index 39eb13c44ae0e9f29bdb5a67f32fddceaa8c108e..46f5c60f4343d1fd840311f4acdeaf9446f153e3 100644
--- a/Geo/GenericVertex.h
+++ b/Geo/GenericVertex.h
@@ -1,4 +1,4 @@
-// Gmsh - Copyright (C) 1997-2015 C. Geuzaine, J.-F. Remacle
+// Gmsh - Copyright (C) 1997-2014 C. Geuzaine, J.-F. Remacle
 //
 // See the LICENSE.txt file for license information. Please report all
 // bugs and problems to the public mailing list <gmsh@geuz.org>.
@@ -9,43 +9,66 @@
 #include "GmshConfig.h"
 #include "GModel.h"
 #include "GVertex.h"
+#include "Context.h"
 #include <vector>
+#include <iostream>
 
 using namespace std;
 
-/* The set of Generic Entities is a generic interface to any other modeler.
-  Callbacks (function pointers) are given, sending requests, enquiries, to the
-  native modeler. */
+/*The set of Generic Entities is a generic interface to any other modeler.
+  Callbacks (function pointers) are given, sending requests, enquiries, to the native modeler. */
 
 class GenericVertex : public GVertex {
- protected:
-  int id;
-  double _x, _y, _z;
- public:
-  // callbacks typedef
-  typedef bool (*ptrfunction_int_vector)(int, vector<double>&);
-
-  GenericVertex(GModel *m, int num, int _native_id);
-  virtual ~GenericVertex();
-
-  virtual GPoint point() const { return GPoint(x(), y(), z()); }
-  virtual double x() const { return _x; }
-  virtual double y() const { return _y; }
-  virtual double z() const { return _z; }
-
-  virtual void setPosition(GPoint &p);
-  virtual SPoint2 reparamOnFace(const GFace *gf, int) const;
-
-  ModelType getNativeType() const { return GenericModel; }
-  virtual int getNativeInt()const{return id;};
-
-  // sets the callbacks
-  static void setVertexXYZ(ptrfunction_int_vector fct){VertexXYZ = fct;};
-
-private:
-  // the callbacks:
-  // fills vector xyz for vertex of int id
-  static ptrfunction_int_vector VertexXYZ;
+  protected:
+    int id;
+    double _x, _y, _z;
+  public:
+    // callbacks typedef
+    typedef bool (*ptrfunction_int_vector)(int, vector<double>&);
+    typedef bool (*ptrfunction_int_doubleptr_voidptr)(int, double*, void*);
+
+    GenericVertex(GModel *m, int num, int _native_id);
+    GenericVertex(GModel *m, int num, int _native_id, const vector<double> &vec);
+    virtual ~GenericVertex();
+
+    virtual GPoint point() const { return GPoint(x(), y(), z()); }
+    virtual double x() const { return _x; }
+    virtual double y() const { return _y; }
+    virtual double z() const { return _z; }
+
+    virtual void setPosition(GPoint &p);
+    virtual SPoint2 reparamOnFace(const GFace *gf, int) const;
+
+    ModelType getNativeType() const { return GenericModel; }
+    virtual int getNativeInt()const{return id;};
+
+    // sets the callbacks
+    static void setVertexXYZ(ptrfunction_int_vector fct){VertexXYZ = fct;};
+    static void setVertexMeshSize(ptrfunction_int_doubleptr_voidptr fct){VertexMeshSize = fct;};
+
+    // meshing-related methods:
+    virtual void setPrescribedMeshSizeAtVertex(double l) {
+      cout << "GenericVertex:: setPrescribedMeshSizeAtVertex !!!, aborting" << endl;
+      throw;
+    }
+    virtual inline double prescribedMeshSizeAtVertex() const {
+      double size;
+      void *chose = NULL;
+      if (!VertexMeshSize(id,&size,chose)){
+        cout<< "GenericVertex::ERROR from callback VertexMeshSize "<< endl;
+        Msg::Warning("GenericVertex::ERROR from callback VertexMeshSize ");
+        cout << "Check for ring edges !!!" << endl;
+        return CTX::instance()->lc;
+      }
+      return size;
+    }
+    
+  private:
+    // the callbacks:
+    // --------------
+    // fills vector xyz for vertex of int id
+    static ptrfunction_int_vector VertexXYZ;
+    static ptrfunction_int_doubleptr_voidptr VertexMeshSize;
 };
 
 #endif
diff --git a/Geo/MElement.h b/Geo/MElement.h
index ee814f22e567be51ee4160301cbf78e729efc514..c19d8e4a4add67796c46af3d1cd04c1cfa04365e 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -130,7 +130,7 @@ class MElement
   }
 
   // get the edges
-  virtual int getNumEdges() = 0;
+  virtual int getNumEdges() const = 0;
   virtual MEdge getEdge(int num) const= 0;
 
   // give an MEdge as input and get its local number and sign
diff --git a/Geo/MElementCut.h b/Geo/MElementCut.h
index df265f9e5f7e06491939b195d0e8685bcdd6357c..d96fd3379ef04fbe3114fb708f8537a20493dbda 100644
--- a/Geo/MElementCut.h
+++ b/Geo/MElementCut.h
@@ -70,7 +70,7 @@ class MPolyhedron : public MElement {
   {
     return (num < (int)_vertices.size()) ? _vertices[num] : _innerVertices[num - _vertices.size()];
   }
-  virtual int getNumEdges() { return _edges.size(); }
+  virtual int getNumEdges()const { return _edges.size(); }
   virtual MEdge getEdge(int num) const{ return _edges[num]; }
   virtual int getNumEdgesRep(bool curved) { return _edges.size(); }
   virtual void getEdgeRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n)
@@ -227,7 +227,7 @@ class MPolygon : public MElement {
   {
     return (num < (int)_vertices.size()) ? _vertices[num] : _innerVertices[num - _vertices.size()];
   }
-  virtual int getNumEdges() { return _edges.size(); }
+  virtual int getNumEdges()const { return _edges.size(); }
   virtual MEdge getEdge(int num) const{ return _edges[num]; }
   virtual int getNumEdgesRep(bool curved) { return getNumEdges(); }
   virtual void getEdgeRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n)
diff --git a/Geo/MElementOctree.cpp b/Geo/MElementOctree.cpp
index ffbe1e89970cd7167c5a3842fdbf7e28d506bab5..c49fc04336326811b32ce15c8c8c872ce48ee867 100644
--- a/Geo/MElementOctree.cpp
+++ b/Geo/MElementOctree.cpp
@@ -5,6 +5,7 @@
 
 #include "GModel.h"
 #include "MElement.h"
+#include "MTetrahedron.h"
 #include "MElementOctree.h"
 #include "Octree.h"
 #include "Context.h"
@@ -117,6 +118,42 @@ MElementOctree::MElementOctree(std::vector<MElement*> &v) : _gm(0), _elems(v)
   Octree_Arrange(_octree);
 }
 
+
+MElementOctree::MElementOctree(std::vector<MTetrahedron*> &v) : _gm(0)
+{
+  // exactly the same as MElementOctree::MElementOctree(vector<MElement*>) !!!
+  for (std::vector<MTetrahedron*>::iterator it = v.begin();it!=v.end();it++){
+    _elems.push_back(*it);
+  }
+  // exactly the same as MElementOctree::MElementOctree(vector<MElement*>) !!!
+
+
+  SBoundingBox3d bb;
+  for (unsigned int i = 0; i < v.size(); i++){
+    for(int j = 0; j < v[i]->getNumVertices(); j++){
+      //if (!_gm) _gm = v[i]->getVertex(j)->onWhat()->model();
+      bb += SPoint3(v[i]->getVertex(j)->x(),
+                    v[i]->getVertex(j)->y(),
+                    v[i]->getVertex(j)->z());
+    }
+  }
+  // make bounding box larger up to (absolute) geometrical tolerance
+  double eps = CTX::instance()->geom.tolerance;
+  SPoint3 bbmin = bb.min(), bbmax = bb.max(), bbeps(eps, eps, eps);
+  bbmin -= bbeps;
+  bbmax += bbeps;
+  double min[3] = {bbmin.x(), bbmin.y(), bbmin.z()};
+  double size[3] = {bbmax.x() - bbmin.x(),
+                    bbmax.y() - bbmin.y(),
+                    bbmax.z() - bbmin.z()};
+  const int maxElePerBucket = 100; // memory vs. speed trade-off
+  _octree = Octree_Create(maxElePerBucket, min, size,
+                          MElementBB, MElementCentroid, MElementInEle);
+  for (unsigned int i = 0; i < v.size(); i++)
+    Octree_Insert(v[i], _octree);
+  Octree_Arrange(_octree);
+}
+
 MElementOctree::~MElementOctree()
 {
   Octree_Delete(_octree);
diff --git a/Geo/MElementOctree.h b/Geo/MElementOctree.h
index a517f8b68d211f750489a2e3052a2784caf4374c..61beb1e6a4f6ed11f03bbef1fff64f58dc7d2c54 100644
--- a/Geo/MElementOctree.h
+++ b/Geo/MElementOctree.h
@@ -11,6 +11,7 @@
 class Octree;
 class GModel;
 class MElement;
+class MTetrahedron;
 
 class MElementOctree{
  private:
@@ -20,6 +21,7 @@ class MElementOctree{
  public:
   MElementOctree(GModel *);
   MElementOctree(std::vector<MElement*> &);
+  MElementOctree(std::vector<MTetrahedron*> &);
   ~MElementOctree();
   MElement *find(double x, double y, double z, int dim = -1, bool strict = false) const;
   Octree *getInternalOctree(){ return _octree; }
diff --git a/Geo/MHexahedron.h b/Geo/MHexahedron.h
index 7fa2f57ac7fca046202fbea5c5a796e3d9a114e7..28c5c5b2b0b4d9c2f7edbb8faa5fc679d4adedbb 100644
--- a/Geo/MHexahedron.h
+++ b/Geo/MHexahedron.h
@@ -65,7 +65,7 @@ class MHexahedron : public MElement {
     static const int map[8] = {2, 3, 7, 6, 0, 1, 5, 4};
     return getVertex(map[num]);
   }
-  virtual int getNumEdges(){ return 12; }
+  virtual int getNumEdges()const{ return 12; }
   virtual MEdge getEdge(int num) const
   {
     return MEdge(_v[edges_hexa(num, 0)], _v[edges_hexa(num, 1)]);
diff --git a/Geo/MLine.h b/Geo/MLine.h
index 05f33059ec4c1cc93f295d43043412b036de8350..cb4251b3058da2adeba1c8dd86331ea9881f8605 100644
--- a/Geo/MLine.h
+++ b/Geo/MLine.h
@@ -47,7 +47,7 @@ class MLine : public MElement {
   {
     ithVertex = _v[0] == vertex ? 0 : 1;
   }
-  virtual int getNumEdges(){ return 1; }
+  virtual int getNumEdges()const{ return 1; }
   virtual MEdge getEdge(int num) const{ return MEdge(_v[0], _v[1]); }
   virtual int getNumEdgesRep(bool curved){ return 1; }
   virtual void getEdgeRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n)
diff --git a/Geo/MPoint.h b/Geo/MPoint.h
index a589707eae909ffb25e3397200780d636d431c2a..301f7a2d2d315306064eae84a5ee408b0ee6ed5b 100644
--- a/Geo/MPoint.h
+++ b/Geo/MPoint.h
@@ -33,7 +33,7 @@ class MPoint : public MElement {
   virtual MVertex *getVertex(int num){ return _v[0]; }
   virtual const MVertex *getVertex(int num) const { return _v[0]; }
   virtual void setVertex(int num,  MVertex *v){ _v[0] = v; }
-  virtual int getNumEdges(){ return 0; }
+  virtual int getNumEdges()const{ return 0; }
   virtual MEdge getEdge(int num) const{ return MEdge(); }
   virtual int getNumEdgesRep(bool curved){ return 0; }
   virtual void getEdgeRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n){}
diff --git a/Geo/MPrism.h b/Geo/MPrism.h
index 3df36f57c749517c7b524cd93344d7ddc5cdf730..a827966e1b6e048a65d7db8cb819e71e17900881 100644
--- a/Geo/MPrism.h
+++ b/Geo/MPrism.h
@@ -66,7 +66,7 @@ class MPrism : public MElement {
   virtual MVertex *getVertex(int num){ return _v[num]; }
   virtual const MVertex *getVertex(int num)const{ return _v[num]; }
   virtual void setVertex(int num,   MVertex *v){ _v[num] = v; }
-  virtual int getNumEdges(){ return 9; }
+  virtual int getNumEdges()const{ return 9; }
   virtual MEdge getEdge(int num) const
   {
     return MEdge(_v[edges_prism(num, 0)], _v[edges_prism(num, 1)]);
diff --git a/Geo/MPyramid.h b/Geo/MPyramid.h
index 31a98b0590eb35542e2326a472c5ecc713d068d4..6b26d9ff4a05ac2d2b0ce0c8e092c6c278ac0d05 100644
--- a/Geo/MPyramid.h
+++ b/Geo/MPyramid.h
@@ -70,7 +70,7 @@ class MPyramid : public MElement {
   virtual const MVertex *getVertex(int num) const{ return _v[num]; }
   virtual void setVertex(int num,  MVertex *v){ _v[num] = v; }
   virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
-  virtual int getNumEdges(){ return 8; }
+  virtual int getNumEdges()const{ return 8; }
   virtual MEdge getEdge(int num) const
   {
     return MEdge(_v[edges_pyramid(num, 0)], _v[edges_pyramid(num, 1)]);
diff --git a/Geo/MQuadrangle.h b/Geo/MQuadrangle.h
index c0c27727d4e46960cdb4d15577a270092eb6c4bc..1f04ca4351fa92fd98c853f375e4dcd4cb0c427c 100644
--- a/Geo/MQuadrangle.h
+++ b/Geo/MQuadrangle.h
@@ -63,7 +63,7 @@ class MQuadrangle : public MElement {
     static const int map[4] = {0, 1, 3, 2};
     return getVertex(map[num]);
   }
-  virtual int getNumEdges(){ return 4; }
+  virtual int getNumEdges()const{ return 4; }
   virtual MEdge getEdge(int num) const
   {
     return MEdge(_v[edges_quad(num, 0)], _v[edges_quad(num, 1)]);
diff --git a/Geo/MTetrahedron.h b/Geo/MTetrahedron.h
index f5884bae2b87b473073711222e9af0d0b5b78666..ebed2a7e256591e6723017b0df5ae0f3ae48f808 100644
--- a/Geo/MTetrahedron.h
+++ b/Geo/MTetrahedron.h
@@ -62,7 +62,7 @@ class MTetrahedron : public MElement {
   virtual MVertex *getVertex(int num){ return _v[num]; }
   virtual const MVertex *getVertex(int num) const { return _v[num]; }
   virtual void setVertex(int num,  MVertex *v){ _v[num] = v; }
-  virtual int getNumEdges(){ return 6; }
+  virtual int getNumEdges()const{ return 6; }
   virtual MEdge getEdge(int num) const
   {
     return MEdge(_v[edges_tetra(num, 0)], _v[edges_tetra(num, 1)]);
diff --git a/Geo/MTriangle.h b/Geo/MTriangle.h
index 5f7874839aca8582ce276493a8946bb969791890..f68a3e0a70b0e3cebb4d59ec060e5f73c37b6253 100644
--- a/Geo/MTriangle.h
+++ b/Geo/MTriangle.h
@@ -67,7 +67,7 @@ class MTriangle : public MElement {
     if(_v[2] != v1 && _v[2] != v2) return _v[2];
     return 0;
   }
-  virtual int getNumEdges(){ return 3; }
+  virtual int getNumEdges()const{ return 3; }
   virtual MEdge getEdge(int num) const
   {
     return MEdge(_v[edges_tri(num, 0)], _v[edges_tri(num, 1)]);
diff --git a/Geo/SPoint2.h b/Geo/SPoint2.h
index 06e91dceecaf090dd0297e9f4d60fb0ac608f6d4..f7dab5b8680b446677319a4d4a3c57b116b06661 100644
--- a/Geo/SPoint2.h
+++ b/Geo/SPoint2.h
@@ -29,7 +29,7 @@ class SPoint2 {
   void operator+=(const SPoint2 &p);
   void operator-=(const SPoint2 &p);
   void operator*=(double mult);
-  SPoint2 operator*(double mult);
+  SPoint2 operator*(double mult)const;
   operator double *() { return P; }
   bool operator < (const SPoint2  &other) const
   {
@@ -79,7 +79,7 @@ inline void SPoint2::operator-=(const SPoint2 &p)
 inline void SPoint2::operator*=(double mult)
 { P[0] *= mult; P[1] *= mult; }
 
-inline SPoint2 SPoint2::operator*(double mult)
+inline SPoint2 SPoint2::operator*(double mult)const
 { return SPoint2(P[0]*mult, P[1]*mult); }
 
 inline double SPoint2::distance(const SPoint2 &p)const
diff --git a/Geo/STensor3.h b/Geo/STensor3.h
index ba44057fed2bd25d0ccf598b211d1fcb53778a11..06b0ea5890e09d71210bc3f43aae7995c1c661f2 100644
--- a/Geo/STensor3.h
+++ b/Geo/STensor3.h
@@ -281,6 +281,10 @@ class STensor3 {
     return ithis;
   }
 
+  void operator = (const STensor3 &other){
+    for (int i = 0; i < 9; i++) _val[i] = other._val[i];
+  }
+
   STensor3 operator + (const STensor3 &other) const
   {
     STensor3 res(*this);