diff --git a/Geo/GEdge.h b/Geo/GEdge.h
index d9cf68151692dab28178c81d296cb5a9b6141cbf..834283d28b9dc98bc4d057110e5221fac7a84005 100644
--- a/Geo/GEdge.h
+++ b/Geo/GEdge.h
@@ -103,8 +103,9 @@ class GEdge : public GEntity {
   virtual double prescribedMeshSizeAtVertex() const { return meshAttributes.meshSize; }
 
   // true if start == end and no more than 2 segments
-  void setTooSmall ( bool b) {_tooSmall = b;}
-  bool isMeshDegenerated() const{ return _tooSmall || (v0 == v1 && mesh_vertices.size() < 2); }
+  void setTooSmall(bool b) { _tooSmall = b; }
+  bool isMeshDegenerated() const 
+  { return _tooSmall || (v0 == v1 && mesh_vertices.size() < 2); }
 
   // number of types of elements
   int getNumElementTypes() const { return 1; }
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index a2a4f68e9dbe9f5f0749fc8ae6aedfca80a3de06..3cef450ce890093d3ff26af602e4d571710c9585 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -2426,13 +2426,20 @@ int GModel::writeDIFF(const std::string &name, bool binary, bool saveAll,
   std::vector<GEntity*> entities;
   getEntities(entities);
 
+  // find max dimension of mesh elements we need to save
+  int dim = 0;
+  for(unsigned int i = 0; i < entities.size(); i++)
+    if(entities[i]->physicals.size() || saveAll)
+      for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++)
+        dim = std::max(dim, entities[i]->getMeshElement(j)->getDim());
+  
   // loop over all elements we need to save
   int numElements = 0, maxNumNodesPerElement = 0;
   for(unsigned int i = 0; i < entities.size(); i++){
     if(entities[i]->physicals.size() || saveAll){
       for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
         MElement *e = entities[i]->getMeshElement(j);
-        if(e->getStringForDIFF()){
+        if(e->getStringForDIFF() && e->getDim() == dim){
           numElements++;
           maxNumNodesPerElement = std::max(maxNumNodesPerElement, e->getNumVertices());
         }
@@ -2494,7 +2501,7 @@ int GModel::writeDIFF(const std::string &name, bool binary, bool saveAll,
     if(entities[i]->physicals.size() || saveAll){
       for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
         MElement *e = entities[i]->getMeshElement(j);
-        if(e->getStringForDIFF())
+        if(e->getStringForDIFF() && e->getDim() == dim)
           e->writeDIFF(fp, ++num, binary, entities[i]->tag());
       }
     }
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 3a55d93127bee90dd9e64e3bbbc2ffeda68225f5..1bed85769b36d0b089713e1c1cd151ba9fc47641 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -127,15 +127,20 @@ static double F_Lc(GEdge *ge, double t)
 
 static double F_Transfinite(GEdge *ge, double t)
 {
-  double val, r;
+  double length = ge->length();
+  if(length == 0.0){
+    Msg::Error("Zero-length curve in transfinite mesh");
+    return 1.;
+  }
 
   SVector3 der = ge->firstDer(t) ;
   double d = norm(der);
-
   double coef = ge->meshAttributes.coeffTransfinite;
   int type = ge->meshAttributes.typeTransfinite;
   int nbpt = ge->meshAttributes.nbPointsTransfinite;
 
+  double val;
+
   if(coef <= 0.0 || coef == 1.0) {
     // coef < 0 should never happen
     val = d * coef / ge->length();
@@ -145,12 +150,9 @@ static double F_Transfinite(GEdge *ge, double t)
 
     case 1: // Geometric progression ar^i; Sum of n terms = length = a (r^n-1)/(r-1)
       {
-        if(sign(type) >= 0)
-          r = coef;
-        else
-          r = 1. / coef;
-        double a = ge->length() * (r - 1.) / (pow(r, nbpt - 1.) - 1.);
-        int i = (int)(log(t * ge->length() / a * (r - 1.) + 1.) / log(r));
+        double r = (sign(type) >= 0) ? coef : 1. / coef;
+        double a = length * (r - 1.) / (pow(r, nbpt - 1.) - 1.);
+        int i = (int)(log(t * length / a * (r - 1.) + 1.) / log(r));
         val = d / (a * pow(r, (double)i));
       }
       break;
@@ -159,18 +161,16 @@ static double F_Transfinite(GEdge *ge, double t)
       {
         double a;
         if(coef > 1.0) {
-          a = -4. * sqrt(coef - 1.) *
-            atan2(1., sqrt(coef - 1.)) /
-            ((double)nbpt *  ge->length());
+          a = -4. * sqrt(coef - 1.) * atan2(1., sqrt(coef - 1.)) /
+            ((double)nbpt *  length);
         }
         else {
-          a = 2. * sqrt(1. - coef) *
-            log(fabs((1. + 1. / sqrt(1. - coef))
-                     / (1. - 1. / sqrt(1. - coef))))
-            / ((double)nbpt * ge->length());
+          a = 2. * sqrt(1. - coef) * log(fabs((1. + 1. / sqrt(1. - coef)) /
+                                              (1. - 1. / sqrt(1. - coef))))
+            / ((double)nbpt * length);
         }
-        double b = -a * ge->length() * ge->length() / (4. * (coef - 1.));
-        val = d / (-a * SQU(t * ge->length() - (ge->length()) * 0.5) + b);
+        double b = -a * length * length / (4. * (coef - 1.));
+        val = d / (-a * SQU(t * length - (length) * 0.5) + b);
       }
       break;
       
@@ -288,8 +288,7 @@ void meshGEdge::operator() (GEdge *ge)
     Msg::Debug("Curve %d has a zero length", ge->tag());
 
   // TEST
-  if (length < CTX.mesh.tolerance_edge_length)ge->setTooSmall(true);
-
+  if (length < CTX.mesh.tolerance_edge_length) ge->setTooSmall(true);
 
   // Integrate detJ/lc du 
   double a;