diff --git a/Geo/discreteRegion.cpp b/Geo/discreteRegion.cpp
index 34b995889bfc50d663d20e8ec1dfd71d23fc2226..6c8bdfcba5b4234b37de1573b330a8f83db8e9e6 100644
--- a/Geo/discreteRegion.cpp
+++ b/Geo/discreteRegion.cpp
@@ -10,6 +10,7 @@
 #include "Context.h"
 
 #if defined(HAVE_MESH)
+#include "meshGRegion.h"
 #include "meshGRegionDelaunayInsertion.h"
 #endif
 
@@ -61,13 +62,29 @@ void discreteRegion::findFaces(std::map<MFace, std::vector<int>, Less_Face> &map
 void discreteRegion::remesh()
 {
 #if defined(HAVE_MESH)
+
+  bool classify = false;
   if(CTX::instance()->mesh.oldRefinement){
-    insertVerticesInRegion(this, 2000000000, true);
+    insertVerticesInRegion(this, 2000000000, classify);
   }
   else{
-    insertVerticesInRegion(this, 0);
+    insertVerticesInRegion(this, 0, classify);
     void edgeBasedRefinement(const int numThreads, const int nptsatonce, GRegion *gr);
     edgeBasedRefinement(1, 1, this);
   }
+
+  // not functional yet: need boundaries
+  for(int i = 0; i < std::max(CTX::instance()->mesh.optimize,
+                              CTX::instance()->mesh.optimizeNetgen); i++){
+    if(CTX::instance()->mesh.optimize >= i){
+      printf("optimizing!\n");
+      optimizeMeshGRegionGmsh opt; opt(this, true);
+    }
+    if(CTX::instance()->mesh.optimizeNetgen >= i){
+      printf("optimizing netgen!\n");
+      optimizeMeshGRegionNetgen opt; opt(this, true);
+    }
+  }
+
 #endif
 }
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 9a1d3dd2d9089e4c0eb796caa611b2264d27edc9..7e913150f8f95a7448ce7cf4090d9727bd789566 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -997,11 +997,11 @@ void meshGRegion::operator() (GRegion *gr)
 
 }
 
-void optimizeMeshGRegionNetgen::operator() (GRegion *gr)
+void optimizeMeshGRegionNetgen::operator() (GRegion *gr, bool always)
 {
   gr->model()->setCurrentMeshEntity(gr);
 
-  if(gr->geomType() == GEntity::DiscreteVolume) return;
+  if(!always && gr->geomType() == GEntity::DiscreteVolume) return;
 
   // don't optimize transfinite or extruded meshes
   if(gr->meshAttributes.method == MESH_TRANSFINITE) return;
@@ -1026,11 +1026,11 @@ void optimizeMeshGRegionNetgen::operator() (GRegion *gr)
 #endif
 }
 
-void optimizeMeshGRegionGmsh::operator() (GRegion *gr)
+void optimizeMeshGRegionGmsh::operator() (GRegion *gr, bool always)
 {
   gr->model()->setCurrentMeshEntity(gr);
 
-  if(gr->geomType() == GEntity::DiscreteVolume) return;
+  if(!always && gr->geomType() == GEntity::DiscreteVolume) return;
 
   // don't optimize extruded meshes
   if(gr->meshAttributes.method == MESH_TRANSFINITE) return;
diff --git a/Mesh/meshGRegion.h b/Mesh/meshGRegion.h
index 25f12794bf80024d06a1a804e77db1550a2152ce..1d50f5c5e40807f1319a541645bae96cd4ba6586 100644
--- a/Mesh/meshGRegion.h
+++ b/Mesh/meshGRegion.h
@@ -35,13 +35,13 @@ class meshGRegionExtruded {
 // Optimize the mesh of the region using gmsh's algo
 class optimizeMeshGRegionGmsh {
  public :
-  void operator () (GRegion *);
+  void operator () (GRegion *, bool always=false);
 };
 
 // Optimize the mesh of the region using netgen's algo
 class optimizeMeshGRegionNetgen {
  public :
-  void operator () (GRegion *);
+  void operator () (GRegion *, bool always=false);
 };
 
 // destroy the mesh of the region
diff --git a/utils/api_demos/mainRemesh.cpp b/utils/api_demos/mainRemesh.cpp
index f6219494a9ac741e5c23900d2c6910076c3654b1..69713b774e20ef7f1524d3af011d46a70bed4ccc 100644
--- a/utils/api_demos/mainRemesh.cpp
+++ b/utils/api_demos/mainRemesh.cpp
@@ -7,11 +7,15 @@
 int main(int argc, char **argv)
 {
   GmshInitialize(argc, argv);
+  GmshSetOption("General", "Terminal", 1.);
+  GmshSetOption("General", "Verbosity", 4.);
   GmshSetOption("Mesh", "CharacteristicLengthExtendFromBoundary", 1.);
   GmshSetOption("Mesh", "OldRefinement", 1.);
   GmshSetOption("Mesh", "CharacteristicLengthMin", 0.1);
   GmshSetOption("Mesh", "CharacteristicLengthMax", 0.1);
 
+  GmshSetOption("Mesh", "Optimize", 0.); // not yet: need boundary!
+
   GModel *m = new GModel();
   m->readMSH("cube.msh");