diff --git a/Geo/Cell.cpp b/Geo/Cell.cpp
index 86dd29c8bfaccf12c178d6910b25f8b827a2cbd0..f59cc141332e9154eba9df04d172ffa75fbd0eba 100755
--- a/Geo/Cell.cpp
+++ b/Geo/Cell.cpp
@@ -299,56 +299,45 @@ CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell()
   }
 
   // boundary cells
-  std::map< Cell*, int, Less_Cell > c1Boundary;
-  c1->getBoundary(c1Boundary);
-  std::map< Cell*, int, Less_Cell > c2Boundary; 
-  c2->getBoundary(c2Boundary);
-  
-  for(biter it = c1Boundary.begin(); it != c1Boundary.end(); it++){
+  for(biter it = c1->firstBoundary(); it != c1->lastBoundary(); it++){
     Cell* cell = (*it).first;
     int ori = (*it).second;
-    cell->removeCoboundaryCell(c1); 
+    cell->removeCoboundaryCell(c1, false); 
     this->addBoundaryCell(ori, cell, false, true);
   }
-  for(biter it = c2Boundary.begin(); it != c2Boundary.end(); it++){
+  for(biter it = c2->firstBoundary(); it != c2->lastBoundary(); it++){
     Cell* cell = (*it).first;
     if(!orMatch) (*it).second = -1*(*it).second;
     int ori = (*it).second;
-    cell->removeCoboundaryCell(c2);    
-    if(co){
-      biter it2 = c1Boundary.find(cell);
-      if(it2 == c1Boundary.end()){
-        this->addBoundaryCell(ori, cell, false, true);
-      }
+    cell->removeCoboundaryCell(c2, false);    
+    if(co && !c1->hasBoundary(cell)){
+      this->addBoundaryCell(ori, cell, false, true);
     }
-    else this->addBoundaryCell(ori, cell, false, true);
+    else if(!co) this->addBoundaryCell(ori, cell, false, true);
   }
+  c1->clearBoundary();
+  c2->clearBoundary();
 
   // coboundary cells
-  std::map<Cell*, int, Less_Cell > c1Coboundary;
-  c1->getCoboundary(c1Coboundary);
-  std::map<Cell*, int, Less_Cell > c2Coboundary;
-  c2->getCoboundary(c2Coboundary);
-  
-  for(biter it = c1Coboundary.begin(); it != c1Coboundary.end(); it++){
+  for(biter it = c1->firstCoboundary(); it != c1->lastCoboundary(); it++){
     Cell* cell = (*it).first;
     int ori = (*it).second;
-    cell->removeBoundaryCell(c1); 
+    cell->removeBoundaryCell(c1, false); 
     this->addCoboundaryCell(ori, cell, false, true);
   }
-  for(biter it = c2Coboundary.begin(); it != c2Coboundary.end(); it++){
+  for(biter it = c2->firstCoboundary(); it != c2->lastCoboundary(); it++){
     Cell* cell = (*it).first;
     if(!orMatch) (*it).second = -1*(*it).second;
     int ori = (*it).second;
-    cell->removeBoundaryCell(c2);    
-    if(!co){
-      biter it2 = c1Coboundary.find(cell);
-      if(it2 == c1Coboundary.end()){
-        this->addCoboundaryCell(ori, cell, false, true);
-      }
+    cell->removeBoundaryCell(c2, false);    
+    if(!co && !c1->hasCoboundary(cell)){
+      this->addCoboundaryCell(ori, cell, false, true);
     }
-    else this->addCoboundaryCell(ori, cell, false, true);
+    else if(co) this->addCoboundaryCell(ori, cell, false, true);
   }
+  c1->clearCoboundary();
+  c2->clearCoboundary();
+  
 
 }
 
diff --git a/Geo/Cell.h b/Geo/Cell.h
index 7bafb9ddcfcf6d582d670289b0c1dd87300cb55c..46b516af3b0cafbcdbe3ab70880692b691ea9fd2 100644
--- a/Geo/Cell.h
+++ b/Geo/Cell.h
@@ -145,13 +145,13 @@ class Cell
   
   // add (co)boundary cell
   void addBoundaryCell(int orientation, Cell* cell, 
-		       bool orig=false, bool other=true); 
+		       bool orig, bool other); 
   void addCoboundaryCell(int orientation, Cell* cell, 
-			 bool orig=false, bool other=true);
+			 bool orig, bool other);
   
   // remove (co)boundary cell
-  void removeBoundaryCell(Cell* cell, bool other=true);
-  void removeCoboundaryCell(Cell* cell, bool other=true);
+  void removeBoundaryCell(Cell* cell, bool other);
+  void removeCoboundaryCell(Cell* cell, bool other);
   
   // true if has given cell on (original) (co)boundary
   bool hasBoundary(Cell* cell, bool orig=false);
diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 4de9f9904924a84f6f044c6b5d5f83e01b79fb1b..4a87e357193fd655be500ee0ff0bc6cf78cbac6e 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -398,11 +398,7 @@ int CellComplex::cocombine(int dim)
         Cell* c2 = (*it).first;
 
         if(!(*c1 == *c2) && abs(or1) == abs(or2)
-           && inSameDomain(s, c1) && inSameDomain(s, c2)
-           && c1->getNumSortedVertices() < getSize(dim) 
-           // heuristics for mammoth cell birth control
-           && c2->getNumSortedVertices() < getSize(dim)){
-
+           && inSameDomain(s, c1) && inSameDomain(s, c2)){
           removeCell(s);
           
           c1->getCoboundary(cbd_c);
@@ -460,11 +456,7 @@ int CellComplex::combine(int dim)
         Cell* c2 = (*it).first;
 
         if(!(*c1 == *c2) && abs(or1) == abs(or2)
-           && inSameDomain(s, c1) && inSameDomain(s, c2)
-           && c1->getNumSortedVertices() < getSize(dim) 
-           // heuristics for mammoth cell birth control
-           && c2->getNumSortedVertices() < getSize(dim)){
-
+           && inSameDomain(s, c1) && inSameDomain(s, c2)){
           removeCell(s);
           
           c1->getBoundary(bd_c);
@@ -506,7 +498,7 @@ bool CellComplex::coherent()
         citer cit = _cells[bdCell->getDim()].find(bdCell);
         if(cit == lastCell(bdCell->getDim())){ 
           printf("Warning! Boundary cell not in cell complex! Boundary removed. \n");
-          cell->removeBoundaryCell(bdCell);
+          cell->removeBoundaryCell(bdCell, false);
           coherent = false;
         }
         if(!bdCell->hasCoboundary(cell)){
@@ -525,7 +517,7 @@ bool CellComplex::coherent()
         citer cit = _cells[cbdCell->getDim()].find(cbdCell);
         if(cit == lastCell(cbdCell->getDim())){ 
           printf("Warning! Coboundary cell not in cell complex! Coboundary removed. \n");
-          cell->removeCoboundaryCell(cbdCell);
+          cell->removeCoboundaryCell(cbdCell, false);
           coherent = false;
         }
         if(!cbdCell->hasBoundary(cell)){