diff --git a/Geo/GModelIO_GEO.cpp b/Geo/GModelIO_GEO.cpp
index 00b9dab94064864aacf26151640264fc9da593ac..cfc6d4e27382fbe3cfc62f4c1dd990a65f701946 100644
--- a/Geo/GModelIO_GEO.cpp
+++ b/Geo/GModelIO_GEO.cpp
@@ -41,6 +41,7 @@ void GEO_Internals::_allocateAll()
   Volumes = Tree_Create(sizeof(Volume *), CompareVolume);
 
   PhysicalGroups = List_Create(5, 5, sizeof(PhysicalGroup *));
+  DelPhysicalGroups = List_Create(5, 5, sizeof(PhysicalGroup *));
 
   DelPoints = Tree_Create(sizeof(Vertex *), CompareVertex);
   DelCurves = Tree_Create(sizeof(Curve *), CompareCurve);
@@ -68,6 +69,7 @@ void GEO_Internals::_freeAll()
   Tree_Action(DelVolumes, FreeVolume); Tree_Delete(DelVolumes);
 
   List_Action(PhysicalGroups, FreePhysicalGroup); List_Delete(PhysicalGroups);
+  List_Action(DelPhysicalGroups, FreePhysicalGroup); List_Delete(DelPhysicalGroups);
 
   _changed = true;
 }
@@ -468,6 +470,7 @@ void GEO_Internals::_transform(std::vector<int> tags[4], int mode,
   case 2: DilatShapes(x, y, z, a, b, c, list); break;
   case 3: SymmetryShapes(a, b, c, d, list); break;
   }
+  _changed = true;
 }
 
 void GEO_Internals::translate(std::vector<int> tags[4],
@@ -498,6 +501,7 @@ void GEO_Internals::symmetry(std::vector<int> tags[4],
 
 int GEO_Internals::copy(int dim, int tag)
 {
+  _changed = true;
   if(dim == 0){
     Vertex *v = FindPoint(tag);
     if(!v){
@@ -549,6 +553,7 @@ void GEO_Internals::remove(int dim, int tag)
 void GEO_Internals::resetPhysicalGroups()
 {
   List_Action(PhysicalGroups, FreePhysicalGroup);
+  List_Action(DelPhysicalGroups, FreePhysicalGroup);
   List_Reset(PhysicalGroups);
   _changed = true;
 }
@@ -600,6 +605,7 @@ void GEO_Internals::modifyPhysicalGroup(int dim, int num, int op, std::vector<in
   else{
     Msg::Error("Unsupported operation on physical %s %d", str.c_str(), num);
   }
+  _changed = true;
 }
 
 void GEO_Internals::removeAllDuplicates()
@@ -683,6 +689,7 @@ void GEO_Internals::setTransfiniteLine(int tag, int nPoints, int type, double co
       c->coeffTransfinite = coef;
     }
   }
+  _changed = true;
 }
 
 void GEO_Internals::setTransfiniteSurface(int tag, int arrangement,
@@ -719,6 +726,7 @@ void GEO_Internals::setTransfiniteSurface(int tag, int arrangement,
       }
     }
   }
+  _changed = true;
 }
 
 void GEO_Internals::setTransfiniteVolume(int tag, std::vector<int> cornerTags)
@@ -749,6 +757,7 @@ void GEO_Internals::setTransfiniteVolume(int tag, std::vector<int> cornerTags)
       }
     }
   }
+  _changed = true;
 }
 
 void GEO_Internals::setTransfiniteVolumeQuadTri(int tag)
@@ -767,6 +776,7 @@ void GEO_Internals::setTransfiniteVolumeQuadTri(int tag)
     if(v)
       v->QuadTri = TRANSFINITE_QUADTRI_1;
   }
+  _changed = true;
 }
 
 void GEO_Internals::setRecombine(int dim, int tag, double angle)
@@ -807,6 +817,7 @@ void GEO_Internals::setRecombine(int dim, int tag, double angle)
       }
     }
   }
+  _changed = true;
 }
 
 void GEO_Internals::setSmoothing(int tag, int val)
@@ -824,6 +835,7 @@ void GEO_Internals::setSmoothing(int tag, int val)
     Surface *s = FindSurface(tag);
     if(s) s->TransfiniteSmoothing = val;
   }
+  _changed = true;
 }
 
 void GEO_Internals::setReverseMesh(int dim, int tag)
@@ -858,12 +870,18 @@ void GEO_Internals::setReverseMesh(int dim, int tag)
       if(s) s->ReverseMesh = 1;
     }
   }
+  _changed = true;
 }
 
 void GEO_Internals::synchronize(GModel *model)
 {
   Msg::Debug("Syncing GEO_Internals with GModel");
 
+  // if the entities do not exist, we create them; if they exist, we update the
+  // pointer (and the underlying dependencies, e.g. surface boundaries): this is
+  // necessary because a GEO entity can change (while keeping the same tag!!),
+  // due e.g. to ReplaceDuplicates.
+
   if(Tree_Nbr(Points)) {
     List_T *points = Tree2List(Points);
     for(int i = 0; i < List_Nbr(points); i++){
@@ -875,6 +893,8 @@ void GEO_Internals::synchronize(GModel *model)
         model->add(v);
       }
       else{
+        if(v->getNativeType() == GEntity::GmshModel)
+          ((gmshVertex*)v)->resetNativePtr(p);
         v->resetMeshAttributes();
       }
     }
@@ -894,8 +914,7 @@ void GEO_Internals::synchronize(GModel *model)
                      "have been created", c->Num);
         }
         else if(!e && c->beg && c->end){
-          e = new gmshEdge(model, c,
-                           model->getVertexByTag(c->beg->Num),
+          e = new gmshEdge(model, c, model->getVertexByTag(c->beg->Num),
                            model->getVertexByTag(c->end->Num));
           model->add(e);
         }
@@ -904,6 +923,10 @@ void GEO_Internals::synchronize(GModel *model)
           model->add(e);
         }
         else{
+          if(e->getNativeType() == GEntity::GmshModel &&
+             c->Typ != MSH_SEGM_COMPOUND)
+            ((gmshEdge*)e)->resetNativePtr(c, model->getVertexByTag(c->beg->Num),
+                                           model->getVertexByTag(c->end->Num));
           e->resetMeshAttributes();
         }
         if(c->degenerated) e->setTooSmall(true);
@@ -958,13 +981,13 @@ void GEO_Internals::synchronize(GModel *model)
         }
         int param = CTX::instance()->mesh.remeshParam;
         GFaceCompound::typeOfCompound typ = GFaceCompound::HARMONIC_CIRCLE;
-        if (param == 1) typ =  GFaceCompound::CONFORMAL_SPECTRAL;
-        if (param == 2) typ =  GFaceCompound::RADIAL_BASIS;
-        if (param == 3) typ =  GFaceCompound::HARMONIC_PLANE;
-        if (param == 4) typ =  GFaceCompound::CONVEX_CIRCLE;
-        if (param == 5) typ =  GFaceCompound::CONVEX_PLANE;
-        if (param == 6) typ =  GFaceCompound::HARMONIC_SQUARE;
-        if (param == 7) typ =  GFaceCompound::CONFORMAL_FE;
+        if (param == 1) typ = GFaceCompound::CONFORMAL_SPECTRAL;
+        if (param == 2) typ = GFaceCompound::RADIAL_BASIS;
+        if (param == 3) typ = GFaceCompound::HARMONIC_PLANE;
+        if (param == 4) typ = GFaceCompound::CONVEX_CIRCLE;
+        if (param == 5) typ = GFaceCompound::CONVEX_PLANE;
+        if (param == 6) typ = GFaceCompound::HARMONIC_SQUARE;
+        if (param == 7) typ = GFaceCompound::CONFORMAL_FE;
         int algo = CTX::instance()->mesh.remeshAlgo;
         f = new GFaceCompound(model, s->Num, comp, b[0], b[1], b[2], b[3], typ, algo);
         f->meshAttributes.recombine = s->Recombine;
@@ -989,8 +1012,9 @@ void GEO_Internals::synchronize(GModel *model)
         model->add(f);
       }
       else{
-        // recompute in case geom has changed
-        if(s->Typ == MSH_SURF_PLAN) f->computeMeanPlane();
+        if(f->getNativeType() == GEntity::GmshModel &&
+           s->Typ != MSH_SURF_COMPOUND)
+          ((gmshFace*)f)->resetNativePtr(s);
         f->resetMeshAttributes();
       }
     }
@@ -1017,6 +1041,9 @@ void GEO_Internals::synchronize(GModel *model)
         model->add(r);
       }
       else{
+        if(r->getNativeType() == GEntity::GmshModel &&
+           v->Typ != MSH_VOLUME_COMPOUND)
+          ((gmshRegion*)r)->resetNativePtr(v);
         r->resetMeshAttributes();
       }
     }
diff --git a/Geo/GModelIO_GEO.h b/Geo/GModelIO_GEO.h
index 3eca3918725be5d219b565ca5a7910401bcdf523..fc4e0be5bed309709e0dda0f86e403ccc62b645c 100644
--- a/Geo/GModelIO_GEO.h
+++ b/Geo/GModelIO_GEO.h
@@ -18,7 +18,7 @@ class GEO_Internals{
   // is complete
   Tree_T *Points, *Curves, *EdgeLoops, *Surfaces, *SurfaceLoops, *Volumes;
   Tree_T *DelPoints, *DelCurves, *DelSurfaces, *DelVolumes;
-  List_T *PhysicalGroups;
+  List_T *PhysicalGroups, *DelPhysicalGroups;
   int MaxPointNum, MaxLineNum, MaxLineLoopNum, MaxSurfaceNum;
   int MaxSurfaceLoopNum, MaxVolumeNum, MaxPhysicalNum;
   std::multimap<int, std::vector<int> > meshCompounds;
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index 399decbb29a51eadea2b3fa86c4cbf751b8ab781..ece4062682bd2ce15f9e0f61b5ba5fbf6b5fda3e 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -1060,7 +1060,7 @@ void DeletePhysicalPoint(int num)
   if(p){
     List_Suppress(GModel::current()->getGEOInternals()->PhysicalGroups, &p,
                   ComparePhysicalGroup);
-    FreePhysicalGroup(&p, NULL);
+    List_Add(GModel::current()->getGEOInternals()->DelPhysicalGroups, &p);
   }
   GModel::current()->deletePhysicalGroup(0, num);
 }
@@ -1071,7 +1071,7 @@ void DeletePhysicalLine(int num)
   if(p){
     List_Suppress(GModel::current()->getGEOInternals()->PhysicalGroups, &p,
                   ComparePhysicalGroup);
-    FreePhysicalGroup(&p, NULL);
+    List_Add(GModel::current()->getGEOInternals()->DelPhysicalGroups, &p);
   }
   GModel::current()->deletePhysicalGroup(1, num);
 }
@@ -1082,7 +1082,7 @@ void DeletePhysicalSurface(int num)
   if(p){
     List_Suppress(GModel::current()->getGEOInternals()->PhysicalGroups, &p,
                   ComparePhysicalGroup);
-    FreePhysicalGroup(&p, NULL);
+    List_Add(GModel::current()->getGEOInternals()->DelPhysicalGroups, &p);
   }
   GModel::current()->deletePhysicalGroup(2, num);
 }
@@ -1093,7 +1093,7 @@ void DeletePhysicalVolume(int num)
   if(p){
     List_Suppress(GModel::current()->getGEOInternals()->PhysicalGroups, &p,
                   ComparePhysicalGroup);
-    FreePhysicalGroup(&p, NULL);
+    List_Add(GModel::current()->getGEOInternals()->DelPhysicalGroups, &p);
   }
   GModel::current()->deletePhysicalGroup(3, num);
 }
@@ -2007,7 +2007,7 @@ static void ReplaceDuplicatePointsNew(double tol = -1.)
   for(unsigned int i = 0; i < unused.size(); i++){
     Vertex *V = v2V[unused[i]];
     Tree_Suppress(GModel::current()->getGEOInternals()->Points, &V);
-    FreeVertex(&V, NULL);
+    Tree_Add(GModel::current()->getGEOInternals()->DelPoints, &V);
     delete unused[i];
   }
   for(unsigned int i = 0; i < used.size(); i++){
@@ -2182,7 +2182,10 @@ static void ReplaceDuplicatePoints(std::map<int, int> * v_report = 0)
     }
   }
 
-  Tree_Action(points2delete, FreeVertex);
+  List_T *tmp = Tree2List(points2delete);
+  for(int i = 0; i < List_Nbr(tmp); i++)
+    Tree_Add(GModel::current()->getGEOInternals()->DelPoints, List_Pointer(tmp, i));
+  List_Delete(tmp);
   Tree_Delete(points2delete);
   Tree_Delete(allNonDuplicatedPoints);
 
@@ -2208,7 +2211,11 @@ static void ReplaceDuplicateCurves(std::map<int, int> * c_report = 0)
         if(!(c2 = FindCurve(-c->Num))) {
           Msg::Error("Unknown curve %d", -c->Num);
           List_Delete(All);
-          Tree_Action(curves2delete, FreeCurve);
+          List_T *tmp = Tree2List(curves2delete);
+          for(int i = 0; i < List_Nbr(tmp); i++)
+            Tree_Add(GModel::current()->getGEOInternals()->DelCurves,
+                     List_Pointer(tmp, i));
+          List_Delete(tmp);
           Tree_Delete(curves2delete);
           Tree_Delete(allNonDuplicatedCurves);
           return;
@@ -2247,7 +2254,11 @@ static void ReplaceDuplicateCurves(std::map<int, int> * c_report = 0)
   int end = Tree_Nbr(GModel::current()->getGEOInternals()->Curves);
 
   if(start == end) {
-    Tree_Action(curves2delete, FreeCurve);
+    List_T *tmp = Tree2List(curves2delete);
+    for(int i = 0; i < List_Nbr(tmp); i++)
+      Tree_Add(GModel::current()->getGEOInternals()->DelCurves,
+               List_Pointer(tmp, i));
+    List_Delete(tmp);
     Tree_Delete(curves2delete);
     Tree_Delete(allNonDuplicatedCurves);
     return;
@@ -2324,7 +2335,11 @@ static void ReplaceDuplicateCurves(std::map<int, int> * c_report = 0)
     }
   }
 
-  Tree_Action(curves2delete, FreeCurve);
+  List_T *tmp = Tree2List(curves2delete);
+  for(int i = 0; i < List_Nbr(tmp); i++)
+    Tree_Add(GModel::current()->getGEOInternals()->DelCurves,
+             List_Pointer(tmp, i));
+  List_Delete(tmp);
   Tree_Delete(curves2delete);
   Tree_Delete(allNonDuplicatedCurves);
 }
@@ -2525,7 +2540,11 @@ static void ReplaceDuplicateSurfaces(std::map<int, int> *s_report = 0)
   int end = Tree_Nbr(GModel::current()->getGEOInternals()->Surfaces);
 
   if(start == end) {
-    Tree_Action(surfaces2delete, FreeSurface);
+    List_T *tmp = Tree2List(surfaces2delete);
+    for(int i = 0; i < List_Nbr(tmp); i++)
+      Tree_Add(GModel::current()->getGEOInternals()->DelSurfaces,
+               List_Pointer(tmp, i));
+    List_Delete(tmp);
     Tree_Delete(surfaces2delete);
     Tree_Delete(allNonDuplicatedSurfaces);
     return;
@@ -2599,7 +2618,11 @@ static void ReplaceDuplicateSurfaces(std::map<int, int> *s_report = 0)
     }
   }
 
-  Tree_Action(surfaces2delete, FreeSurface);
+  List_T *tmp = Tree2List(surfaces2delete);
+  for(int i = 0; i < List_Nbr(tmp); i++)
+    Tree_Add(GModel::current()->getGEOInternals()->DelSurfaces,
+             List_Pointer(tmp, i));
+  List_Delete(tmp);
   Tree_Delete(surfaces2delete);
   Tree_Delete(allNonDuplicatedSurfaces);
 }
@@ -3091,14 +3114,13 @@ int ExtrudeSurface(int type, int is,
       c->Extrude->mesh = e->mesh;
   }
 
-  // FIXME: this is a really ugly hack for backward compatibility, so
-  // that we don't screw up the old .geo files too much. (Before
-  // version 1.54, we didn't always create new volumes during "Extrude
-  // Surface". Now we do, but with "CTX::instance()->geom.oldNewreg==1", this
-  // bumps the NEWREG() counter, and thus changes the whole automatic
-  // numbering sequence.) So we locally force oldNewreg to 0: in most
-  // cases, since we define points, curves, etc., before defining
-  // volumes, the NEWVOLUME() call below will return a fairly low
+  // FIXME: this is a really ugly hack for backward compatibility, so that we
+  // don't screw up the old .geo files too much. (Before version 1.54, we didn't
+  // always create new volumes during "Extrude Surface". Now we do, but with
+  // "CTX::instance()->geom.oldNewreg==1", this bumps the NEWREG() counter, and
+  // thus changes the whole automatic numbering sequence.) So we locally force
+  // oldNewreg to 0: in most cases, since we define points, curves, etc., before
+  // defining volumes, the NEWVOLUME() call below will return a fairly low
   // number, that will not interfere with the other numbers...
   int tmp = CTX::instance()->geom.oldNewreg;
   CTX::instance()->geom.oldNewreg = 0;
diff --git a/Geo/gmshEdge.cpp b/Geo/gmshEdge.cpp
index 8796caea924720cdfb9bb4a80422fe9194dd1b81..1283ce10c341a60cc4f96354103fe22d3722e265 100644
--- a/Geo/gmshEdge.cpp
+++ b/Geo/gmshEdge.cpp
@@ -20,6 +20,14 @@ gmshEdge::gmshEdge(GModel *m, Curve *edge, GVertex *v1, GVertex *v2)
   resetMeshAttributes();
 }
 
+void gmshEdge::resetNativePtr(Curve *edge, GVertex *_v1, GVertex *_v2)
+{
+  c = edge;
+  v0 = _v1; v1 = _v2;
+  if(v0) v0->addEdge(this);
+  if(v1 && v1 != v0) v1->addEdge(this);
+}
+
 bool gmshEdge::degenerate(int dim) const
 {
   if (c->beg == c->end &&
diff --git a/Geo/gmshEdge.h b/Geo/gmshEdge.h
index 0b01c8b27e22551a1b1f145e83f99068ef89b3fa..26e2c52bafadda0c9a940310f87df0f5478113ab 100644
--- a/Geo/gmshEdge.h
+++ b/Geo/gmshEdge.h
@@ -23,7 +23,7 @@ class gmshEdge : public GEdge {
   virtual SVector3 firstDer(double par) const;
   virtual SVector3 secondDer(double par) const;
   ModelType getNativeType() const { return GmshModel; }
-  void * getNativePtr() const { return c; }
+  void *getNativePtr() const { return c; }
   virtual std::string getAdditionalInfoString();
   virtual int minimumMeshSegments() const;
   virtual int minimumDrawSegments() const;
@@ -32,6 +32,7 @@ class gmshEdge : public GEdge {
   virtual void writeGEO(FILE *fp);
   void discretize(double tol, std::vector<SPoint3> &dpts, std::vector<double> &ts);
   virtual bool degenerate(int dim) const;
+  void resetNativePtr(Curve *edge, GVertex *v1, GVertex *v2);
 };
 
 #endif
diff --git a/Geo/gmshFace.cpp b/Geo/gmshFace.cpp
index ae84d9bc13f1cc1c4ef5d84eacb5db46e31b7652..334e7816e6f4685043904acc29ef502d97c65387 100644
--- a/Geo/gmshFace.cpp
+++ b/Geo/gmshFace.cpp
@@ -21,19 +21,23 @@
 #endif
 
 gmshFace::gmshFace(GModel *m, Surface *face)
-  : GFace(m, face->Num), s(face), isSphere(false), radius(0.)
+  : GFace(m, face->Num), isSphere(false), radius(0.)
 {
+  resetNativePtr(face);
   resetMeshAttributes();
+}
 
-  // edgeCounterparts = s->edgeCounterparts;
-  // affineTransform = s->affineTransform;
-
+void gmshFace::resetNativePtr(Surface *face)
+{
+  s = face;
+  l_edges.clear();
+  l_dirs.clear();
   std::vector<GEdge*> eds;
   std::vector<int> nums;
   for(int i = 0; i < List_Nbr(s->Generatrices); i++){
     Curve *c;
     List_Read(s->Generatrices, i, &c);
-    GEdge *e = m->getEdgeByTag(abs(c->Num));
+    GEdge *e = model()->getEdgeByTag(abs(c->Num));
     if(e){
       eds.push_back(e);
       nums.push_back(c->Num);
@@ -44,7 +48,7 @@ gmshFace::gmshFace(GModel *m, Surface *face)
   for(int i = 0; i < List_Nbr(s->GeneratricesByTag); i++){
     int j;
     List_Read(s->GeneratricesByTag, i, &j);
-    GEdge *e = m->getEdgeByTag(abs(j));
+    GEdge *e = model()->getEdgeByTag(abs(j));
     if(e){
       eds.push_back(e);
       nums.push_back(j);
diff --git a/Geo/gmshFace.h b/Geo/gmshFace.h
index 9c118ba7a349b1efc9bb83512f01c3ad3eb93b5d..466da57068cb190fba4916c837df67cb6b095be2 100644
--- a/Geo/gmshFace.h
+++ b/Geo/gmshFace.h
@@ -35,6 +35,7 @@ class gmshFace : public GFace {
   void *getNativePtr() const { return s; }
   virtual SPoint2 parFromPoint(const SPoint3 &, bool onSurface=true) const;
   virtual void resetMeshAttributes();
+  void resetNativePtr(Surface *_s);
 };
 
 #endif
diff --git a/Geo/gmshRegion.cpp b/Geo/gmshRegion.cpp
index 9b4983da5b6a93553199d429153116f290984417..890654c441c09e381e079f6917b01da50884cf27 100644
--- a/Geo/gmshRegion.cpp
+++ b/Geo/gmshRegion.cpp
@@ -10,14 +10,21 @@
 #include "GmshMessage.h"
 
 gmshRegion::gmshRegion(GModel *m, ::Volume *volume)
-  : GRegion(m, volume->Num), v(volume)
+  : GRegion(m, volume->Num)
 {
+  resetNativePtr(volume);
+  resetMeshAttributes();
+}
+
+void gmshRegion::resetNativePtr(::Volume *volume)
+{
+  v = volume;
   for(int i = 0; i < List_Nbr(v->Surfaces); i++){
     Surface *s;
     List_Read(v->Surfaces, i, &s);
     int ori;
     List_Read(v->SurfacesOrientations, i, &ori);
-    GFace *f = m->getFaceByTag(abs(s->Num));
+    GFace *f = model()->getFaceByTag(abs(s->Num));
     if(f){
       l_faces.push_back(f);
       l_dirs.push_back(ori);
@@ -29,7 +36,7 @@ gmshRegion::gmshRegion(GModel *m, ::Volume *volume)
   for(int i = 0; i < List_Nbr(v->SurfacesByTag); i++){
     int is;
     List_Read(v->SurfacesByTag, i, &is);
-    GFace *f = m->getFaceByTag(abs(is));
+    GFace *f = model()->getFaceByTag(abs(is));
     if(f){
       l_faces.push_back(f);
       l_dirs.push_back(gmsh_sign(is));
@@ -38,7 +45,6 @@ gmshRegion::gmshRegion(GModel *m, ::Volume *volume)
     else
       Msg::Error("Unknown surface %d", is);
   }
-  resetMeshAttributes();
 }
 
 void gmshRegion::resetMeshAttributes()
diff --git a/Geo/gmshRegion.h b/Geo/gmshRegion.h
index 1ea9cf10fa0cbb547dd8fdab8c50027ebcec2892..43d6b3b33ac71a79f7afb005b8886efd97e71fd1 100644
--- a/Geo/gmshRegion.h
+++ b/Geo/gmshRegion.h
@@ -18,8 +18,9 @@ class gmshRegion : public GRegion {
   virtual ~gmshRegion() {}
   virtual GeomType geomType() const;
   ModelType getNativeType() const { return GmshModel; }
-  void * getNativePtr() const { return v; }
+  void *getNativePtr() const { return v; }
   virtual void resetMeshAttributes();
+  void resetNativePtr(::Volume *_v);
 };
 
 #endif
diff --git a/Geo/gmshVertex.cpp b/Geo/gmshVertex.cpp
index 2bef70159d5ba0851b49ef88b880446f3f68e905..f6c27d2708d0b42d00dd576356122866c954fb6e 100644
--- a/Geo/gmshVertex.cpp
+++ b/Geo/gmshVertex.cpp
@@ -18,6 +18,11 @@ gmshVertex::gmshVertex(GModel *m, Vertex *_v)
   resetMeshAttributes();
 }
 
+void gmshVertex::resetNativePtr(Vertex *_v)
+{
+  v = _v;
+}
+
 void gmshVertex::resetMeshAttributes()
 {
   meshSize = v->lc;
diff --git a/Geo/gmshVertex.h b/Geo/gmshVertex.h
index cc947775b0c60695bb526a2bb83fa84c25968314..84f8bb878af449c33cd4d5775589306edf7c2f29 100644
--- a/Geo/gmshVertex.h
+++ b/Geo/gmshVertex.h
@@ -29,6 +29,7 @@ class gmshVertex : public GVertex {
   virtual void setPrescribedMeshSizeAtVertex(double l);
   virtual SPoint2 reparamOnFace(const GFace *gf, int) const;
   virtual void writeGEO(FILE *fp, const std::string &meshSizeParameter="");
+  void resetNativePtr(Vertex *_v);
 };
 
 #endif
diff --git a/Parser/Gmsh.tab.cpp b/Parser/Gmsh.tab.cpp
index 6ba76f55fbbc9dee0518b780e70a759204f78b11..5d3d3384a29b0447f760d706bf75be4504d80aaa 100644
--- a/Parser/Gmsh.tab.cpp
+++ b/Parser/Gmsh.tab.cpp
@@ -14212,12 +14212,8 @@ void setColor(std::vector<int> tags[4], unsigned int val, bool recursive)
   if(GModel::current()->getOCCInternals() &&
      GModel::current()->getOCCInternals()->getChanged())
     GModel::current()->getOCCInternals()->synchronize(GModel::current());
-
-  // FIXME: sync in the middle of some geo files leads to crashes: disabling
-  // until I figure it out
-
-  //if(GModel::current()->getGEOInternals()->getChanged())
-  //    GModel::current()->getGEOInternals()->synchronize(GModel::current());
+  if(GModel::current()->getGEOInternals()->getChanged())
+    GModel::current()->getGEOInternals()->synchronize(GModel::current());
 
   for(int dim = 0; dim < 4; dim++){
     for(unsigned int i = 0; i < tags[dim].size(); i++){
diff --git a/Parser/Gmsh.y b/Parser/Gmsh.y
index 68150ba2ff4bf21e756fb6bd9c5eca73f7e0b3f0..a55b82c2ce7829df7a8c37ff85d9e507fb8d4daf 100644
--- a/Parser/Gmsh.y
+++ b/Parser/Gmsh.y
@@ -6739,12 +6739,8 @@ void setColor(std::vector<int> tags[4], unsigned int val, bool recursive)
   if(GModel::current()->getOCCInternals() &&
      GModel::current()->getOCCInternals()->getChanged())
     GModel::current()->getOCCInternals()->synchronize(GModel::current());
-
-  // FIXME: sync in the middle of some geo files leads to crashes: disabling
-  // until I figure it out
-
-  //if(GModel::current()->getGEOInternals()->getChanged())
-  //    GModel::current()->getGEOInternals()->synchronize(GModel::current());
+  if(GModel::current()->getGEOInternals()->getChanged())
+    GModel::current()->getGEOInternals()->synchronize(GModel::current());
 
   for(int dim = 0; dim < 4; dim++){
     for(unsigned int i = 0; i < tags[dim].size(); i++){