diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 26a3f83b6363efea75779d11ea81b4a483e9687a..b47dfc3c4d1588655005869369bb01fec39e7b82 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -1382,7 +1382,7 @@ bool GFace::fillPointCloud(double maxDist,
 static void meshCompound (GFace* gf, bool verbose) {
 
   discreteFace *df = new discreteFace (gf->model(), gf->tag() + 100000);
-  
+
   std::set<int> ec;
   for (unsigned int i=0;i<gf->_compound.size();i++){
     GFace *c = (GFace*)gf->_compound[i];
@@ -1410,15 +1410,14 @@ static void meshCompound (GFace* gf, bool verbose) {
 void GFace::mesh(bool verbose)
 {
 #if defined(HAVE_MESH)
-
   meshGFace mesher;
   mesher(this, verbose);
-  if (!_compound.empty()){ // Some faces are meshed together 
+  if (!_compound.empty()){ // Some faces are meshed together
     if (_compound[0] == this){ //  I'm the one that makes the compound job
       bool ok = true;
       for (unsigned int i=0;i<_compound.size();i++){
 	GFace *gf = (GFace*)_compound[i];
-	ok &= (gf->meshStatistics.status == GFace::DONE);            
+	ok &= (gf->meshStatistics.status == GFace::DONE);
       }
       if (!ok)meshStatistics.status = GFace::PENDING;
       else {
@@ -1427,7 +1426,6 @@ void GFace::mesh(bool verbose)
       }
     }
   }
-    
 #endif
 }
 
@@ -1575,7 +1573,7 @@ void GFace::setMeshMaster(GFace* master, const std::vector<double>& tfo)
   for (mvIter=m_vertices.begin();mvIter!=m_vertices.end();++mvIter) {
 
     GVertex* m_vertex = *mvIter;
-    
+
     SPoint3 xyzOri((*mvIter)->x(),(*mvIter)->y(),(*mvIter)->z());
     SPoint3 xyzTfo(0,0,0);
 
@@ -1584,7 +1582,7 @@ void GFace::setMeshMaster(GFace* master, const std::vector<double>& tfo)
       for (int j=0;j<3;j++) xyzTfo[i] += xyzOri[j] * tfo[idx++];
       xyzTfo[i] += tfo[idx++];
     }
-    
+
     GVertex* l_vertex = NULL;
 
     std::set<GVertex*>::iterator lvIter = l_vertices.begin();
@@ -1645,7 +1643,7 @@ void GFace::setMeshMaster(GFace* master, const std::vector<double>& tfo)
       return;
     }
     GEdge* masterEdge = mv2eIter->second;
-    
+
     if (masterEdge->meshMaster() != localEdge) {
       localEdge->setMeshMaster(masterEdge,tfo);
       Msg::Info("Setting edge master %d - %d",
diff --git a/Geo/discreteRegion.cpp b/Geo/discreteRegion.cpp
index 5a27ec58274f8e16153699916aa5a0879eaf641b..34b995889bfc50d663d20e8ec1dfd71d23fc2226 100644
--- a/Geo/discreteRegion.cpp
+++ b/Geo/discreteRegion.cpp
@@ -7,6 +7,11 @@
 #include "discreteRegion.h"
 #include "MVertex.h"
 #include "Geo.h"
+#include "Context.h"
+
+#if defined(HAVE_MESH)
+#include "meshGRegionDelaunayInsertion.h"
+#endif
 
 discreteRegion::discreteRegion(GModel *model, int num) : GRegion(model, num)
 {
@@ -53,3 +58,16 @@ void discreteRegion::findFaces(std::map<MFace, std::vector<int>, Less_Face> &map
   }
 }
 
+void discreteRegion::remesh()
+{
+#if defined(HAVE_MESH)
+  if(CTX::instance()->mesh.oldRefinement){
+    insertVerticesInRegion(this, 2000000000, true);
+  }
+  else{
+    insertVerticesInRegion(this, 0);
+    void edgeBasedRefinement(const int numThreads, const int nptsatonce, GRegion *gr);
+    edgeBasedRefinement(1, 1, this);
+  }
+#endif
+}
diff --git a/Geo/discreteRegion.h b/Geo/discreteRegion.h
index 7d8c58bbae2bba15990f8ed844b07c47bbe5d327..684f170a4b133818606685b2c9c17ce0e66d40ef 100644
--- a/Geo/discreteRegion.h
+++ b/Geo/discreteRegion.h
@@ -17,6 +17,7 @@ class discreteRegion : public GRegion {
   virtual GeomType geomType() const { return DiscreteVolume; }
   void setBoundFaces(std::set<int> tagFaces);
   void findFaces(std::map<MFace, std::vector<int>, Less_Face> &map_faces);
+  virtual void remesh();
 };
 
 #endif
diff --git a/utils/api_demos/CMakeLists.txt b/utils/api_demos/CMakeLists.txt
index 82b2f9683a1902e67ebf7ca2af2db9091690a1a0..021bf4c47576b940b19ac6d2aa2a2046fcb1caff 100644
--- a/utils/api_demos/CMakeLists.txt
+++ b/utils/api_demos/CMakeLists.txt
@@ -56,3 +56,6 @@ target_link_libraries(mainSimple shared)
 add_executable(mainGeoFactory mainGeoFactory.cpp)
 target_link_libraries(mainGeoFactory shared)
 
+add_executable(mainRemesh mainRemesh.cpp)
+target_link_libraries(mainRemesh shared)
+
diff --git a/utils/api_demos/mainRemesh.cpp b/utils/api_demos/mainRemesh.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f6219494a9ac741e5c23900d2c6910076c3654b1
--- /dev/null
+++ b/utils/api_demos/mainRemesh.cpp
@@ -0,0 +1,29 @@
+#include <stdio.h>
+#include "Gmsh.h"
+#include "GModel.h"
+#include "MElement.h"
+#include "discreteRegion.h"
+
+int main(int argc, char **argv)
+{
+  GmshInitialize(argc, argv);
+  GmshSetOption("Mesh", "CharacteristicLengthExtendFromBoundary", 1.);
+  GmshSetOption("Mesh", "OldRefinement", 1.);
+  GmshSetOption("Mesh", "CharacteristicLengthMin", 0.1);
+  GmshSetOption("Mesh", "CharacteristicLengthMax", 0.1);
+
+  GModel *m = new GModel();
+  m->readMSH("cube.msh");
+
+  for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it){
+    discreteRegion *r = dynamic_cast<discreteRegion*>(*it);
+    if(r){
+      printf("found a discrete region (%d) with %d tetrahedra\n", r->tag(), r->tetrahedra.size());
+      r->remesh();
+      printf("after remeshing: %d tetrahedra\n", r->tetrahedra.size());
+    }
+  }
+  m->writeMSH("new.msh");
+  delete m;
+  GmshFinalize();
+}