diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 5b57a935abe4a602070e5091e4ea01ca56a966e9..87a8515e670f9008227eceee1283d44dae7c867e 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -497,21 +497,20 @@ int CellComplex::reduceComplex(){
 }
 
 
-void CellComplex::removeSubdomain(){
-  
+void CellComplex::removeSubdomain()
+{
+  std::vector<Cell*> toRemove;
   for(int i = 0; i < 4; i++){
-    for(citer cit = firstCell(i); cit != lastCell(i); ){
-      Cell* cell = *cit;
-      cit++;
-      if(cell->inSubdomain()) removeCell(cell);
+    for(citer cit = firstCell(i); cit != lastCell(i); ++cit){
+      Cell *cell = *cit;
+      if(cell->inSubdomain()) toRemove.push_back(cell);
     }
   }
-
-  return;
+  for(unsigned int i = 0; i < toRemove.size(); i++) removeCell(toRemove[i]);
 }
 
-
-int CellComplex::coreduceComplex(){
+int CellComplex::coreduceComplex()
+{
 
   Msg::Debug("Cell complex before coreduction: %d volumes, %d faces, %d edges and %d vertices.\n",
          getSize(3), getSize(2), getSize(1), getSize(0));