diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index d1b536d0c40ca9cc4e118b23bfcf171a4545eb3d..871418679f20d58516e4f56061efd38c808ff1a7 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -91,6 +91,7 @@ std::vector<std::pair<std::string, std::string> > GetUsage()
   s.push_back(mp("-smooth int",        "Set number of mesh smoothing steps"));
   s.push_back(mp("-order int",         "Set mesh order (1, ..., 5)"));
   s.push_back(mp("-optimize[_netgen]", "Optimize quality of tetrahedral elements"));
+  s.push_back(mp("-optimize_threshold", "Optimize tetrahedral elements that have a qulaity less than a threshold"));
   s.push_back(mp("-optimize_ho",       "Optimize high order meshes"));
   s.push_back(mp("-ho_[min,max,nlayers]", "High-order optimization parameters"));
   s.push_back(mp("-optimize_lloyd",    "Optimize 2D meshes using Lloyd algorithm"));
@@ -480,9 +481,19 @@ void GetOptions(int argc, char *argv[])
         i++;
       }
       else if(!strcmp(argv[i] + 1, "optimize")) {
+	Msg::Warning("The optimize option is now obsolete");
+	Msg::Warning("Gmsh optimizes tetrahedral meshes by default");
+	Msg::Warning("Use the \"-optimize_threshold threshold \" to control which elements are optimized");
+	Msg::Warning("Option \"-optimize_threshold 0 \" leads to no optimization");
         CTX::instance()->mesh.optimize = 1;
         i++;
       }
+      else if(!strcmp(argv[i] + 1, "optimize_threshold")) {
+        i++;
+        if(argv[i])
+          opt_mesh_optimize_threshold(0, GMSH_SET, atof(argv[i++]));
+        i++;
+      }
       else if(!strcmp(argv[i] + 1, "optimize_netgen")) {
         CTX::instance()->mesh.optimizeNetgen = 1;
         i++;
diff --git a/Common/Context.h b/Common/Context.h
index 54a57f6079139adc91975be965aca509c7d5cb35..cadd17160ee36fe0656430dd1db7cbe6229a1c82 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -24,8 +24,8 @@ struct contextMeshOptions {
   int pyramids, trihedra;
   int surfacesEdges, surfacesFaces, volumesEdges, volumesFaces, numSubEdges;
   int pointsNum, linesNum, surfacesNum, volumesNum, qualityType, labelType;
-  int optimize, optimizeNetgen, optimizeLloyd, smoothCrossField, refineSteps;
-  double normals, tangents, explode, angleSmoothNormals, allowSwapEdgeAngle;
+  int optimize,  optimizeNetgen, optimizeLloyd, smoothCrossField, refineSteps;
+  double optimizeThreshold,normals, tangents, explode, angleSmoothNormals, allowSwapEdgeAngle;
   double mshFileVersion, mshFilePartitioned, pointSize, lineWidth;
   double qualityInf, qualitySup, radiusInf, radiusSup;
   double scalingFactor, lcFactor, randFactor, lcIntegrationPrecision;
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index aae3baa8b15539ea8ea41b55c7c7fb078e2fea82..8b0c970a0aa9751e46ee7d08e67ef0045c8e048c 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -1151,9 +1151,11 @@ StringXNumber MeshOptions_Number[] = {
 
   { F|O, "OldRefinement" , opt_mesh_old_refinement , 1 ,
     "Use old 3D point insertion algorithm" },
-  { F|O, "Optimize" , opt_mesh_optimize , 0. ,
+  { F|O, "Optimize" , opt_mesh_optimize , 1. ,
     "Optimize the mesh to improve the quality of tetrahedral elements" },
-  { F|O, "OptimizeNetgen" , opt_mesh_optimize_netgen , 0. ,
+  { F|O, "OptimizeThreshold" , opt_mesh_optimize_threshold , 0.3 ,
+    "Optimize tetrahedra that have a quality below ... " },
+  { F|O, "OptimizeNetgen" , opt_mesh_optimize_netgen , 0 ,
     "Optimize the mesh using Netgen to improve the quality of tetrahedral "
     "elements" },
 
diff --git a/Common/Options.cpp b/Common/Options.cpp
index c14aad51efdfc87840a64510df30352188c1981b..5a66df744ceeb5b3267debaf4afe0b62e795a326 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -4888,6 +4888,17 @@ double opt_mesh_optimize(OPT_ARGS_NUM)
   return CTX::instance()->mesh.optimize;
 }
 
+double opt_mesh_optimize_threshold(OPT_ARGS_NUM)
+{
+  if(action & GMSH_SET){
+    if(!(action & GMSH_SET_DEFAULT) && val != CTX::instance()->mesh.optimizeThreshold)
+      Msg::SetOnelabChanged(2);
+    CTX::instance()->mesh.optimizeThreshold = val;
+  }
+  return CTX::instance()->mesh.optimizeThreshold;
+}
+
+
 double opt_mesh_optimize_netgen(OPT_ARGS_NUM)
 {
   if(action & GMSH_SET){
diff --git a/Common/Options.h b/Common/Options.h
index 0fcad56213ba5b6f04a801c78a232c8409d7e08e..7a4393e9967437a0917cd9c1a69f037168efb769 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -403,6 +403,7 @@ double opt_geometry_exact_extrusion(OPT_ARGS_NUM);
 double opt_geometry_match_geom_and_mesh(OPT_ARGS_NUM);
 double opt_mesh_label_sampling(OPT_ARGS_NUM);
 double opt_mesh_optimize(OPT_ARGS_NUM);
+double opt_mesh_optimize_threshold(OPT_ARGS_NUM);
 double opt_mesh_optimize_netgen(OPT_ARGS_NUM);
 double opt_mesh_old_refinement(OPT_ARGS_NUM);
 double opt_mesh_refine_steps(OPT_ARGS_NUM);
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index 1512aad5f4c5c6fff2e267012966a8102e92b51c..955dccd9e616793e53d676d7d1d424e05a6b1963 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -788,7 +788,11 @@ void adaptMeshGRegion::operator () (GRegion *gr)
 
 void optimizeMesh(GRegion *gr, const qmTetrahedron::Measures &qm)
 {
-  // well, this should not be true !!!
+  double qMin = CTX::instance()->mesh.optimizeThreshold;
+
+  if (qMin <= 0.0) return;
+  
+   // well, this should not be true !!!
   // if (gr->hexahedra.size() ||
   //      gr->prisms.size() ||
   //      gr->pyramids.size())return;
@@ -853,7 +857,6 @@ void optimizeMesh(GRegion *gr, const qmTetrahedron::Measures &qm)
     }
   }
 
-  double qMin = 0.3;
   double sliverLimit = 0.04;
 
   int nbESwap = 0, nbFSwap = 0, nbReloc = 0;
@@ -1146,7 +1149,6 @@ void insertVerticesInRegion (GRegion *gr, int maxVert, bool _classify)
   { // leave this in a block so the map gets deallocated directly
     std::map<MVertex*, double,MVertexLessThanNum> vSizesMap;
     std::set<MVertex*,MVertexLessThanNum> bndVertices;
-
     std::list<GEdge*> e = gr->embeddedEdges();
     for (std::list<GEdge*>::iterator it = e.begin() ; it != e.end(); ++it){
       for (unsigned int i = 0; i < (*it)->lines.size(); i++){