diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index 6e7b80af3f075dad0fd47affcb63d4ceb8c91590..8ea69893b401621c2b41107e8cfd152529c7921f 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -305,7 +305,8 @@ const double GOLDEN2 = 2 - GOLDEN;
 // x1 and x3 are the current bounds; the minimum is between them.
 // x2 is the center point, which is closer to x1 than to x3
 
-double goldenSectionSearch(const GEdge *ge, const SPoint3 &q, double x1, double x2, double x3, double tau)
+double goldenSectionSearch(const GEdge *ge, const SPoint3 &q, double x1, 
+                           double x2, double x3, double tau)
 { 
   // Create a new possible center in the area between x2 and x3, closer to x2
   double x4 = x2 + GOLDEN2 * (x3 - x2);
@@ -328,7 +329,6 @@ double goldenSectionSearch(const GEdge *ge, const SPoint3 &q, double x1, double
 
 GPoint GEdge::closestPoint(const SPoint3 &q, double &t) const
 {
-
   //  printf("looking for closest point in curve %d to point %g %g\n",tag(),q.x(),q.y());
   
   const int nbSamples = 100;
@@ -434,9 +434,9 @@ bool GEdge::XYZToU(const double X, const double Y, const double Z,
   return false;
 }
 
-void GEdge::replaceEndingPoints (GVertex *replOfv0, GVertex *replOfv1)
+void GEdge::replaceEndingPoints(GVertex *replOfv0, GVertex *replOfv1)
 {
-  replaceEndingPointsInternals (replOfv0, replOfv1);
+  replaceEndingPointsInternals(replOfv0, replOfv1);
   if (replOfv0 != v0){
     v0->delEdge(this);
     replOfv0->addEdge(this);
@@ -448,4 +448,3 @@ void GEdge::replaceEndingPoints (GVertex *replOfv0, GVertex *replOfv1)
     v1 = replOfv1;
   }  
 }
-
diff --git a/Geo/GEdge.h b/Geo/GEdge.h
index 77caffbccf50dea3fe4f300328730726373d152e..f967c2ffdf818abc44e637a09861b840f8168730 100644
--- a/Geo/GEdge.h
+++ b/Geo/GEdge.h
@@ -77,7 +77,7 @@ class GEdge : public GEntity {
 
   // get the point for the given parameter location
   virtual GPoint point(double p) const = 0;
-  
+
   // true if the edge contains the given parameter
   virtual bool containsParam(double pt) const;
 
@@ -181,16 +181,6 @@ class GEdge : public GEntity {
   // gluing
   void replaceEndingPoints(GVertex *, GVertex *);
 
-  //get bounds
-  inline double getLowBound() const {  
-    Range<double> bounds = parBounds(0);
-    return bounds.low();
-  }
-  inline double getHighBound() const {  
-    Range<double> bounds = parBounds(0);
-    return bounds.high();
-  }
-
   struct {
     char Method;
     double coeffTransfinite;
diff --git a/Geo/GEntity.h b/Geo/GEntity.h
index 004c1bdea4f7266cac25a3f2d3d8725235ae5966..cc1fdd0410c95bef6ce5f76142070e9651875d8d 100644
--- a/Geo/GEntity.h
+++ b/Geo/GEntity.h
@@ -222,6 +222,9 @@ class GEntity {
   // true if there are parametric degeneracies in the "dim" direction.
   virtual bool degenerate(int dim) const { return false; }
 
+  // does the entity have a parametrization?
+  virtual bool haveParametrization(){ return true; }
+
   // parametric bounds of the entity in the "i" direction.
   virtual Range<double> parBounds(int i) const { return Range<double>(0., 0.); }
 
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 671a0f82784349c0d70b8c7e33d5247a7dfc38d2..2653288b72d246670404a5a79c35a7cc902ed67f 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -518,8 +518,8 @@ int GModel::mesh(int dimension)
 #endif
 }
 
-
-int GModel::adaptMesh(int technique, simpleFunction<double> *f, std::vector<double> parameters, bool meshAll)
+int GModel::adaptMesh(int technique, simpleFunction<double> *f, 
+                      std::vector<double> parameters, bool meshAll)
 {
 #if defined(HAVE_MESH)
 
@@ -1472,7 +1472,6 @@ void GModel::createTopologyFromMesh(int ignoreHoles)
   double t1 = Cpu();
 
   removeDuplicateMeshVertices(CTX::instance()->geom.tolerance);
-
   makeDiscreteRegionsSimplyConnected();
   makeDiscreteFacesSimplyConnected();
 
@@ -1484,13 +1483,13 @@ void GModel::createTopologyFromMesh(int ignoreHoles)
   createTopologyFromRegions(discRegions);
  
   // create topology for all discrete faces
-   std::vector<discreteFace*> discFaces;
-   for(fiter it = firstFace(); it != lastFace(); it++)
-     if((*it)->geomType() == GEntity::DiscreteSurface)
-       discFaces.push_back((discreteFace*) *it);
-   createTopologyFromFaces(discFaces, ignoreHoles);
-
-  //create old format (necessary for boundary layers)
+  std::vector<discreteFace*> discFaces;
+  for(fiter it = firstFace(); it != lastFace(); it++)
+    if((*it)->geomType() == GEntity::DiscreteSurface)
+      discFaces.push_back((discreteFace*) *it);
+  createTopologyFromFaces(discFaces, ignoreHoles);
+  
+  //create old format (necessary e.g. for old-style extruded boundary layers)
   exportDiscreteGEOInternals();
  
   double t2 = Cpu();
@@ -1723,7 +1722,7 @@ void GModel::createTopologyFromFaces(std::vector<discreteFace*> &discFaces, int
         for (unsigned int i = 0; i < (*itE)->getNumMeshElements(); i++){
           MEdge me = (*itE)->getMeshElement(i)->getEdge(0);
           std::set<MEdge, Less_Edge >::iterator itset = myEdges.find(me);
-          myEdges.erase(itset);
+          if (itset != myEdges.end()) myEdges.erase(itset);
         }
         for (std::vector<int>::iterator itFace = tagFaces.begin(); 
              itFace != tagFaces.end(); itFace++) {
@@ -1754,7 +1753,7 @@ void GModel::createTopologyFromFaces(std::vector<discreteFace*> &discFaces, int
       }
       std::vector<std::vector<MEdge> > new_boundaries;
       new_boundaries.push_back(boundaries[index]);
-      boundaries =  new_boundaries;
+      boundaries = new_boundaries;
     }
 
     // create new discrete edges
@@ -1835,9 +1834,6 @@ void GModel::createTopologyFromFaces(std::vector<discreteFace*> &discFaces, int
   //    edgeLoops.push_back(el);
   //  }
 
-  Msg::Debug("Done creating topology for edges...");
-
-
   // we need to recreate all mesh elements because some mesh vertices
   // might have been changed during the parametrization process
   // (MVertices became MEdgeVertices)
@@ -1911,7 +1907,7 @@ void GModel::createTopologyFromFaces(std::vector<discreteFace*> &discFaces, int
     gr->pyramids = newPyramids;
   }
 
-  Msg::Debug("Done creating topology from edges");
+  Msg::Debug("Done creating topology for edges");
 }
 
 GModel *GModel::buildCutGModel(gLevelset *ls, bool cutElem, bool saveTri)
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index b59e0e4422456faea67090c587141b70c89721e0..3b8c67ffb0f0584640294b5a5f2c3d2fcecf59e9 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -295,6 +295,7 @@ static void getAllParameters(MVertex *v, GFace *gf, std::vector<SPoint2> &params
   }
   else if(v->onWhat()->dim() == 1){
     GEdge *ge = (GEdge*)v->onWhat();
+    if(!ge->haveParametrization()) return;
     double UU;
     v->getParameter(0, UU);
     if (UU == 0.0)
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index fb515254a0de658d6424284feb32c59e69592cac..f8a2f58be3a38858139cbdc4b6c2e8189e0aff1b 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -47,7 +47,7 @@ void discreteEdge::orderMLines()
   std::vector<MLine*> _m;
   std::list<MLine*> segments;
 
-  //store all lines in a list : segments
+  // store all lines in a list : segments
   for (unsigned int i = 0; i < lines.size(); i++){
     segments.push_back(lines[i]);
   }
@@ -65,7 +65,7 @@ void discreteEdge::orderMLines()
     else boundv.erase(it2);
   }
 
-  //find the first MLine and erase it from the list segments
+  // find the first MLine and erase it from the list segments
   MLine *firstLine;
   if (boundv.size() == 2){   // non periodic
     firstLine = (boundv.begin())->second;
@@ -86,7 +86,7 @@ void discreteEdge::orderMLines()
                tag(), boundv.size());
   }
 
-  //loop over all segments to order segments and store it in the list _m
+  // loop over all segments to order segments and store it in the list _m
   _m.push_back(firstLine);
   _orientation.push_back(1);
   MVertex *first = _m[0]->getVertex(0);
@@ -257,7 +257,6 @@ void discreteEdge::setBoundVertices()
 
   v0->addEdge(this);
   v1->addEdge(this);
-
 }
 
 /*
@@ -312,7 +311,7 @@ void discreteEdge::parametrize(std::map<GFace*, std::map<MVertex*, MVertex*,
    for(std::list<GFace*>::iterator iFace = l_faces.begin(); 
        iFace != l_faces.end(); ++iFace){
 
-     //for each face find correspondane face2Vertex
+     // for each face find correspondane face2Vertex
      std::map<GFace*, std::map<MVertex*, MVertex*, 
        std::less<MVertex*> > >::iterator itmap = face2Vert.find(*iFace);
      if (itmap == face2Vert.end()) {
@@ -326,7 +325,7 @@ void discreteEdge::parametrize(std::map<GFace*, std::map<MVertex*, MVertex*,
        itmap->second = mapVert;
      }
 
-     //do the same for regions
+     // do the same for regions
      for ( int j = 0; j < (*iFace)->numRegions(); j++){
        GRegion *r = (*iFace)->getRegion(j);
        std::map<GRegion*,  std::map<MVertex*, MVertex*, 
@@ -399,7 +398,7 @@ void discreteEdge::computeNormals () const
 bool discreteEdge::getLocalParameter(const double &t, int &iLine,
                                      double &tLoc) const
 {
-
+  if(_pars.empty()) return false;
   for (iLine = 0; iLine < (int)lines.size(); iLine++){
     double tmin = _pars[iLine];
     double tmax = _pars[iLine+1];
@@ -416,7 +415,7 @@ GPoint discreteEdge::point(double par) const
 {
   double tLoc;
   int iEdge;
-  if(!getLocalParameter(par,iEdge,tLoc)) return GPoint();
+  if(!getLocalParameter(par, iEdge, tLoc)) return GPoint();
 
   double x, y, z;
   MVertex *vB = lines[iEdge]->getVertex(0);
@@ -475,9 +474,9 @@ GPoint discreteEdge::point(double par) const
 
 SVector3 discreteEdge::firstDer(double par) const
 {
-   double tLoc;
+  double tLoc;
   int iEdge;
-  getLocalParameter(par, iEdge, tLoc);
+  if(!getLocalParameter(par, iEdge, tLoc)) return SVector3();
 
   MVertex *vB = lines[iEdge]->getVertex(0);
   MVertex *vE = lines[iEdge]->getVertex(1);
@@ -492,10 +491,11 @@ SVector3 discreteEdge::firstDer(double par) const
   return der;
 }
 
-double discreteEdge::curvature(double par) const{
+double discreteEdge::curvature(double par) const
+{
   double tLoc;
   int iEdge;
-  if(!getLocalParameter(par,iEdge,tLoc)) return MAX_LC;
+  if(!getLocalParameter(par, iEdge, tLoc)) return MAX_LC;
 
   double c0, c1;
   Curvature& curvature  = Curvature::getInstance();
@@ -521,5 +521,6 @@ Range<double> discreteEdge::parBounds(int i) const
 void discreteEdge::writeGEO(FILE *fp)
 {
   if (getBeginVertex() && getEndVertex())
-    fprintf(fp, "Discrete Line(%d) = {%d,%d};\n",tag(), getBeginVertex()->tag(), getEndVertex()->tag());  
+    fprintf(fp, "Discrete Line(%d) = {%d,%d};\n", tag(), 
+            getBeginVertex()->tag(), getEndVertex()->tag());  
 }
diff --git a/Geo/discreteEdge.h b/Geo/discreteEdge.h
index ec54b77993005826c759893649dbe5021be2244a..0193c31bc7467829754e8ab18d19625701ce84e9 100644
--- a/Geo/discreteEdge.h
+++ b/Geo/discreteEdge.h
@@ -19,12 +19,14 @@ class discreteEdge : public GEdge {
  public:
   discreteEdge(GModel *model, int num, GVertex *_v0, GVertex *_v1);
   virtual ~discreteEdge() {}
-  bool getLocalParameter(const double &t, int &iEdge, double &tLoc) const;
   virtual GeomType geomType() const { return DiscreteCurve; }
   virtual GPoint point(double p) const;
   virtual SVector3 firstDer(double par) const;
   virtual double curvature(double par) const;
+  virtual bool haveParametrization(){ return !_pars.empty(); }
   virtual Range<double> parBounds(int) const;
+
+  bool getLocalParameter(const double &t, int &iEdge, double &tLoc) const;
   void parametrize(std::map<GFace*, std::map<MVertex*, MVertex*, 
                    std::less<MVertex*> > > &face2Verts, 
                    std::map<GRegion*, std::map<MVertex*, MVertex*,
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index ef5f80f3888b923fecf62234e7b793efb237b753..fd099233c11e9195442cd2377612910810a06850 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -144,8 +144,8 @@ void discreteFace::writeGEO(FILE *fp)
   int count = 0;
   for (std::list<GEdge*>::iterator it = l_edges.begin();
        it != l_edges.end() ;++it){
-    if (count == 0)fprintf(fp, "%d",(*it)->tag());    
-    else fprintf(fp, ",%d",(*it)->tag());    
+    if (count == 0) fprintf(fp, "%d", (*it)->tag());    
+    else fprintf(fp, ",%d", (*it)->tag());    
     count ++;
   }
   fprintf(fp, "};\n");