diff --git a/Geo/Cell.cpp b/Geo/Cell.cpp
index a7a63fdbd3c3e2a14f640bbb3cfcefafb362b558..86dd29c8bfaccf12c178d6910b25f8b827a2cbd0 100755
--- a/Geo/Cell.cpp
+++ b/Geo/Cell.cpp
@@ -12,24 +12,22 @@
 
 bool Less_Cell::operator()(const Cell* c1, const Cell* c2) const 
 {  
-  //cells with fever vertices first
+  if(c1->isCombined() && !c2->isCombined()) return false;
+  if(!c1->isCombined() && c2->isCombined()) return true;
+  if(c1->isCombined() && c2->isCombined()){
+    return (c1->getNum() < c2->getNum());
+  }
   
   if(c1->getNumSortedVertices() != c2->getNumSortedVertices()){
     return (c1->getNumSortedVertices() < c2->getNumSortedVertices());
-  }
-
-  // "natural number" -like ordering 
-  // (the number of a vertice corresponds a digit)
-  
+  }    
   for(int i=0; i < c1->getNumSortedVertices();i++){
     if(c1->getSortedVertex(i) < c2->getSortedVertex(i)) return true;
     else if (c1->getSortedVertex(i) > c2->getSortedVertex(i)) return false;
   }
-  
   return false;
 }
 
-
 Cell::Cell(MElement* image) :  
   _combined(false), _index(0), _immune(false), _image(NULL), 
   _delimage(false), _subdomain(false) 
@@ -273,43 +271,31 @@ void Cell::printCoboundary(bool orig)
   }
 }
 
+int CombinedCell::_globalNum = 0;
 CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell() 
 {  
   // use "smaller" cell as c2
-  if(c1->getNumSortedVertices() < c2->getNumSortedVertices()){
+  if(c1->getNumCells() < c2->getNumCells()){
     Cell* temp = c1;
     c1 = c2;
     c2 = temp;
   }
   
+  _num = ++_globalNum;
+
   _index = c1->getIndex();
   _dim = c1->getDim();
   _subdomain = c1->inSubdomain();
   _combined = true;
   _image = NULL;
 
-  // vertices
-  _vs.reserve(c1->getNumSortedVertices() + c2->getNumSortedVertices());
-  for(int i = 0; i < c1->getNumSortedVertices(); i++){
-    _vs.push_back(c1->getSortedVertex(i));
-  }
-  std::vector<int> v;
-  for(int i = 0; i < c2->getNumSortedVertices(); i++){
-    if(!this->hasVertex(c2->getSortedVertex(i))){
-      v.push_back(c2->getSortedVertex(i)); 
-    }
-  }
-  for(unsigned int i = 0; i < v.size(); i++) _vs.push_back(v.at(i));
-  std::sort(_vs.begin(), _vs.end());
- 
   // cells
   c1->getCells(_cells);
-  std::list< std::pair<int, Cell*> > c2Cells;
+  std::map< Cell*, int, Less_Cell > c2Cells;
   c2->getCells(c2Cells);
-  for(std::list< std::pair<int, Cell*> >::iterator it = c2Cells.begin();
-      it != c2Cells.end(); it++){
-    if(!orMatch) (*it).first = -1*(*it).first;
-    _cells.push_back(*it);
+  for(citer cit  = c2Cells.begin(); cit != c2Cells.end(); cit++){
+    if(!orMatch) (*cit).second = -1*(*cit).second;
+    _cells.insert(*cit);
   }
 
   // boundary cells
diff --git a/Geo/Cell.h b/Geo/Cell.h
index 3750a5ae4c5123551d1cdbed918ee97cfae32312..7bafb9ddcfcf6d582d670289b0c1dd87300cb55c 100644
--- a/Geo/Cell.h
+++ b/Geo/Cell.h
@@ -104,18 +104,19 @@ class Cell
   // get boundary cell orientation
   int getFacetOri(Cell* cell);
 
-  virtual int getDim() const { return _dim; };
-  virtual int getIndex() const { return _index; };
-  virtual void setIndex(int index) { _index = index; };
-  virtual void setImmune(bool immune) { _immune = immune; };
-  virtual bool getImmune() const { return _immune; };
-  virtual bool inSubdomain() const { return _subdomain; }
-  virtual void setInSubdomain(bool subdomain)  { _subdomain = subdomain; }
+  int getDim() const { return _dim; };
+  int getIndex() const { return _index; };
+  void setIndex(int index) { _index = index; };
+  void setImmune(bool immune) { _immune = immune; };
+  bool getImmune() const { return _immune; };
+  bool inSubdomain() const { return _subdomain; }
+  void setInSubdomain(bool subdomain)  { _subdomain = subdomain; }
   virtual int getNumSortedVertices() const { return _vs.size(); }
   virtual int getSortedVertex(int vertex) const { return _vs.at(vertex); }
-  
+  virtual int getNum() const { return 0; }
+
   // restores the cell information to its original state before reduction
-  virtual void restoreCell();
+  void restoreCell();
   
   // true if this cell has given vertex
   virtual bool hasVertex(int vertex) const;
@@ -131,33 +132,33 @@ class Cell
   biter lastCoboundary(bool orig=false){ 
     return orig ? _ocbd.end() : _cbd.end(); }
 
-  virtual int getBoundarySize() { return _bd.size(); }
-  virtual int getCoboundarySize() { return _cbd.size(); }
+  int getBoundarySize() { return _bd.size(); }
+  int getCoboundarySize() { return _cbd.size(); }
    
   // get the cell boundary
-  virtual void getBoundary(std::map<Cell*, int, Less_Cell >& boundary, 
-                           bool orig=false){
+  void getBoundary(std::map<Cell*, int, Less_Cell >& boundary, 
+		   bool orig=false){
     orig ? boundary = _obd : boundary =  _bd; }
-  virtual void getCoboundary(std::map<Cell*, int, Less_Cell >& coboundary,
-                             bool orig=false){
+  void getCoboundary(std::map<Cell*, int, Less_Cell >& coboundary,
+		     bool orig=false){
     orig ? coboundary = _ocbd : coboundary = _cbd; }
   
   // add (co)boundary cell
-  virtual void addBoundaryCell(int orientation, Cell* cell, 
-                               bool orig=false, bool other=true); 
-  virtual void addCoboundaryCell(int orientation, Cell* cell, 
-                                 bool orig=false, bool other=true);
+  void addBoundaryCell(int orientation, Cell* cell, 
+		       bool orig=false, bool other=true); 
+  void addCoboundaryCell(int orientation, Cell* cell, 
+			 bool orig=false, bool other=true);
   
   // remove (co)boundary cell
-  virtual void removeBoundaryCell(Cell* cell, bool other=true);
-  virtual void removeCoboundaryCell(Cell* cell, bool other=true);
+  void removeBoundaryCell(Cell* cell, bool other=true);
+  void removeCoboundaryCell(Cell* cell, bool other=true);
   
   // true if has given cell on (original) (co)boundary
-  virtual bool hasBoundary(Cell* cell, bool orig=false);
-  virtual bool hasCoboundary(Cell* cell, bool orig=false);
+  bool hasBoundary(Cell* cell, bool orig=false);
+  bool hasCoboundary(Cell* cell, bool orig=false);
   
-  virtual void clearBoundary() { _bd.clear(); }
-  virtual void clearCoboundary() { _cbd.clear(); }
+  void clearBoundary() { _bd.clear(); }
+  void clearCoboundary() { _cbd.clear(); }
   
   // print cell debug info
   virtual void printCell();
@@ -165,16 +166,19 @@ class Cell
   virtual void printCoboundary(bool orig=false);   
   
   // tools for combined cells
-  virtual bool isCombined() { return _combined; }
-  virtual void getCells( std::list< std::pair<int, Cell*> >& cells) {  
+  bool isCombined() const { return _combined; }
+  virtual void getCells( std::map< Cell*, int, Less_Cell >& cells) {  
     cells.clear();
-    cells.push_back( std::make_pair(1, this)); 
+    cells[this] = 1; 
     return;
   }
-  virtual int getNumCells() {return 1;}
-  
+  virtual int getNumCells() const {return 1;}
+  typedef std::map<Cell*, int, Less_Cell>::iterator citer;
+
   // equivalence
-  bool operator==(const Cell& c2) const {  
+  virtual bool operator==(const Cell& c2) const {  
+    if(this->isCombined() != c2.isCombined()) return false;
+    
     if(this->getNumSortedVertices() != c2.getNumSortedVertices()){
       return false;
     }
@@ -192,18 +196,27 @@ class CombinedCell : public Cell{
   
  private:
   // list of cells this cell is a combination of
-  std::list< std::pair<int, Cell*> > _cells;
-  
+  std::map< Cell*, int, Less_Cell > _cells;
+
+  // combined cell number, used for ordering
+  int _num;
+  static int _globalNum;
+
  public:
   
   CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co=false);
   ~CombinedCell() {}
   
-  void getCells(std::list< std::pair<int, Cell*> >& cells) { cells = _cells; }
-  int getNumCells() {return _cells.size();}
-  
-};
+  void getCells(std::map< Cell*, int, Less_Cell >& cells) { cells = _cells; }
+  int getNumCells() const {return _cells.size();}  
+  int getNum() const { return _num; }
 
+  bool operator==(const Cell& c2) const {
+    if(this->getNum() != c2.getNum()) return false;
+    else return true;
+  }
+};
+  
 #endif
 
 #endif
diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index c85982e77fb5fb3bfa5926eb9473675e0ea3dabf..4de9f9904924a84f6f044c6b5d5f83e01b79fb1b 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -109,7 +109,12 @@ void CellComplex::insertCell(Cell* cell)
 {
   _newcells.push_back(cell);      
   std::pair<citer, bool> insertInfo = _cells[cell->getDim()].insert(cell);
-  if(!insertInfo.second) printf("Warning: Cell not inserted! \n");
+  if(!insertInfo.second){
+    printf("Warning: Cell not inserted! \n");
+    Cell* oldCell = (*insertInfo.first);
+    cell->printCell();
+    oldCell->printCell();
+  }
 }
 
 void CellComplex::removeCell(Cell* cell, bool other)
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index c65b079c99056fbb7eedc26a461af33b755cdd4f..3bffbd8d68504da11a652c86d860641e1dd7e314 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -458,12 +458,11 @@ void ChainComplex::getChain(std::map<Cell*, int, Less_Cell>& chain,
       continue;
     }
 
-    std::list< std::pair<int, Cell*> > subCells;
+    std::map<Cell*, int, Less_Cell > subCells;
     cell->getCells(subCells);
-    for(std::list< std::pair<int, Cell*> >::iterator it = 
-          subCells.begin(); it != subCells.end(); it++){
-      Cell* subCell = (*it).second;
-      int coeff = (*it).first;
+    for(Cell::citer it = subCells.begin(); it != subCells.end(); it++){
+      Cell* subCell = (*it).first;
+      int coeff = (*it).second;
       chain[subCell] = coeff*elemi;
     }
   }
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index 213594166d795c6d6a63a751d4d95cc39585b0b4..c0e571a8f8252882fcdfe86c7d80531bf746340b 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -446,12 +446,11 @@ Chain::Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs,
     _dim = cell->getDim();
     if((int)coeffs.size() > i){
       if(coeffs.at(i) != 0){
-        std::list< std::pair<int, Cell*> > subCells;
+        std::map<Cell*, int, Less_Cell > subCells;
         cell->getCells(subCells);
-        for(std::list< std::pair<int, Cell*> >::iterator it = 
-              subCells.begin(); it != subCells.end(); it++){
-          Cell* subCell = (*it).second;
-          int coeff = (*it).first;
+        for(Cell::citer it = subCells.begin(); it != subCells.end(); it++){
+          Cell* subCell = (*it).first;
+          int coeff = (*it).second;
           _cells.insert( std::make_pair(subCell, coeffs.at(i)*coeff));
         }
       }
@@ -474,11 +473,10 @@ void Homology::storeCells(CellComplex* cellComplex, int dim)
       cit != cellComplex->lastCell(dim); cit++){
     Cell* cell = *cit;
     
-    std::list< std::pair<int, Cell*> > cells;
+    std::map<Cell*, int, Less_Cell > cells;
     cell->getCells(cells);
-    for(std::list< std::pair<int, Cell*> >::iterator it = cells.begin();
-        it != cells.end(); it++){
-      Cell* subCell = it->second;
+    for(Cell::citer it = cells.begin(); it != cells.end(); it++){
+      Cell* subCell = it->first;
       
       MElement* e = subCell->getImageMElement();
       subCell->setDeleteImage(false);