diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 226a7b1273f2db7a8372128b58107f0ef186c558..109e4a56f99eda5406401ffece4a597a2111222e 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -14,7 +14,7 @@ CellComplex::CellComplex(GModel* model,
                          bool saveOriginalComplex) :
   _model(model), _dim(0), _simplicial(true), _saveorig(saveOriginalComplex)
 {
-
+  _deleteCount = 0;
   _insertCells(subdomainElements, 1);
   _insertCells(domainElements, 0);
 
@@ -76,22 +76,31 @@ CellComplex::~CellComplex()
       for(citer cit = _ocells[i].begin(); cit != _ocells[i].end(); cit++){
         Cell* cell = *cit;
         delete cell;
+        _deleteCount++;
       }
     }
     else {
       for(citer cit = _cells[i].begin(); cit != _cells[i].end(); cit++){
         Cell* cell = *cit;
         delete cell;
+        _deleteCount++;
       }
     }
   }
-  for(unsigned int i = 0; i < _newcells.size(); i++) delete _newcells.at(i);
-  for(unsigned int i = 0; i < _removedcells.size(); i++) delete _removedcells.at(i);
+  for(unsigned int i = 0; i < _newcells.size(); i++) {
+    delete _newcells.at(i);
+    _deleteCount++;
+  }
+  for(unsigned int i = 0; i < _removedcells.size(); i++) {
+    delete _removedcells.at(i);
+    _deleteCount++;
+  }
+  Msg::Debug("Total number of cells created: %d", _deleteCount);
 }
 
 void CellComplex::insertCell(Cell* cell)
 {
-  _newcells.push_back(cell);
+  if(_saveorig) _newcells.push_back(cell);
   std::pair<citer, bool> insertInfo = _cells[cell->getDim()].insert(cell);
   if(!insertInfo.second){
     Msg::Debug("Cell not inserted");
@@ -103,35 +112,26 @@ void CellComplex::insertCell(Cell* cell)
 
 void CellComplex::removeCell(Cell* cell, bool other)
 {
-  //if(!hasCell(cell)) return;
   std::map<Cell*, short int, Less_Cell > coboundary;
   cell->getCoboundary(coboundary);
   std::map<Cell*, short int, Less_Cell > boundary;
   cell->getBoundary(boundary);
 
-  for( std::map<Cell*, short int, Less_Cell>::iterator it = coboundary.begin(); it != coboundary.end(); it++){
+  for(std::map<Cell*, short int, Less_Cell>::iterator it =
+        coboundary.begin(); it != coboundary.end(); it++){
     Cell* cbdCell = (*it).first;
     cbdCell->removeBoundaryCell(cell, other);
   }
 
-  for( std::map<Cell*, short int, Less_Cell>::iterator it = boundary.begin(); it != boundary.end(); it++){
+  for(std::map<Cell*, short int, Less_Cell>::iterator it =
+        boundary.begin(); it != boundary.end(); it++){
     Cell* bdCell = (*it).first;
     bdCell->removeCoboundaryCell(cell, other);
   }
 
-  /* for(Cell::biter it = cell->firstCoboundary(); it != cell->lastCoboundary(); it++){
-    Cell* cbdCell = (*it).first;
-    cbdCell->removeBoundaryCell(cell, false);
-  }
-
-  for(Cell::biter it = cell->firstBoundary(); it != cell->lastBoundary(); it++){
-    Cell* bdCell = (*it).first;
-    bdCell->removeCoboundaryCell(cell, false);
-    }*/
-
   int erased = _cells[cell->getDim()].erase(cell);
   if(!erased) Msg::Debug("Tried to remove a cell from the cell complex \n");
-  if(!_saveorig && !cell->isCombined()) _removedcells.push_back(cell);
+  if(!_saveorig) _removedcells.push_back(cell);
 }
 
 void CellComplex::enqueueCells(std::map<Cell*, short int, Less_Cell>& cells,
@@ -362,7 +362,6 @@ int CellComplex::coreduceComplex(bool docombine, bool omit)
   if(docombine) cocombine(2);
   coreduction(3, false, empty);
   coherent();
-
   Msg::Debug(" %d volumes, %d faces, %d edges and %d vertices",
             getSize(3), getSize(2), getSize(1), getSize(0));
 
@@ -573,11 +572,15 @@ bool CellComplex::restoreComplex()
     for(unsigned int i = 0; i < _newcells.size(); i++){
       Cell* cell = _newcells.at(i);
       delete cell;
+      _deleteCount++;
     }
     _newcells.clear();
     return true;
   }
-  else return false;
+  else {
+    Msg::Error("Cannot restore cell complex");
+    return false;
+  }
 }
 
 void CellComplex::printComplex(int dim)
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index c32124083edff2fcad8574e0e9b0344d70eaf17e..71b362050a8b9e9c14c7cf9a3f1235ca537bf490 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -42,6 +42,8 @@ class CellComplex
   int _dim;
   bool _saveorig;
 
+  int _deleteCount;
+
   // for constructor
   bool _insertCells(std::vector<MElement*>& elements,  int domain);
 
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index f158d5f4bcc9d8d1f003797ccbc35c140026a8c0..af2342ce6c4b6c9d95eb0d2ee5865f223c4ec59f 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -2764,7 +2764,8 @@ void GModel::computeHomology()
     std::pair<std::multimap<dpair, std::string>::iterator,
               std::multimap<dpair, std::string>::iterator> itp =
       _homologyRequests.equal_range(*it);
-    bool prepareToRestore = (itp.first != itp.second);
+    bool prepareToRestore = (itp.first != --itp.second);
+    itp.second++;
     Homology* homology = new Homology(this, itp.first->first.first,
                                       itp.first->first.second, false, prepareToRestore);
     CellComplex *cellcomplex = homology->createCellComplex();