diff --git a/Geo/Cell.cpp b/Geo/Cell.cpp
index 379b5bdff590c8d9a590b0785f71a56b35ca7603..c0a27460b719eba5f6d3f918d13c99722f2db0ac 100755
--- a/Geo/Cell.cpp
+++ b/Geo/Cell.cpp
@@ -36,8 +36,8 @@ Cell::Cell(MElement* image) :
 {
   _dim = image->getDim();
   _image = image;
-  for(int i = 0; i < image->getNumVertices(); i++) 
-    _vs.push_back(image->getVertex(i)->getNum()); 
+  for(int i = 0; i < getNumVertices(); i++) 
+    _vs.push_back(getVertex(i)->getNum()); 
   std::sort(_vs.begin(), _vs.end());
 }
 
@@ -167,65 +167,58 @@ void Cell::restoreCell(){
   _immune = false;   
 }
 
-bool Cell::addBoundaryCell(int orientation, Cell* cell, 
+void Cell::addBoundaryCell(int orientation, Cell* cell, 
 			   bool orig, bool other) 
 {
   if(orig) _obd.insert( std::make_pair(cell, orientation ) );
-  if(other) cell->addCoboundaryCell(orientation, this);
   biter it = _bd.find(cell);
   if(it != _bd.end()){
-    (*it).second = (*it).second + orientation;
-    if((*it).second == 0) {
+    int newOrientation = (*it).second + orientation;
+    if(newOrientation != 0) (*it).second = newOrientation;
+    else {
       _bd.erase(it);
       (*it).first->removeCoboundaryCell(this,false);
-      return false;
+      return;
     }
-    return true;
   }
-  _bd.insert( std::make_pair(cell, orientation ) );
-  return true;
+  else _bd.insert( std::make_pair(cell, orientation ) );
+  if(other) cell->addCoboundaryCell(orientation, this, orig, false);
 }
 
-bool Cell::addCoboundaryCell(int orientation, Cell* cell, 
+void Cell::addCoboundaryCell(int orientation, Cell* cell, 
 			     bool orig, bool other) 
 {
   if(orig) _ocbd.insert( std::make_pair(cell, orientation ) );
-  if(other) cell->addBoundaryCell(orientation, this);
   biter it = _cbd.find(cell);
   if(it != _cbd.end()){
-    (*it).second = (*it).second + orientation;
-    if((*it).second == 0) {
+    int newOrientation = (*it).second + orientation;
+    if(newOrientation != 0) (*it).second = newOrientation;
+    else {
       _cbd.erase(it);
       (*it).first->removeBoundaryCell(this,false);
-      return false;
+      return;
     }
-    return true;
   }
-  _cbd.insert( std::make_pair(cell, orientation ) );
-  return true;
+  else _cbd.insert( std::make_pair(cell, orientation ) );
+  if(other) cell->addBoundaryCell(orientation, this, orig, false);
 }
 
-int Cell::removeBoundaryCell(Cell* cell, bool other) 
+void Cell::removeBoundaryCell(Cell* cell, bool other) 
 {
   biter it = _bd.find(cell);
   if(it != _bd.end()){
     _bd.erase(it);
     if(other) (*it).first->removeCoboundaryCell(this, false);
-    return (*it).second;
   }
-  
-  return 0;
 }
  
-int Cell::removeCoboundaryCell(Cell* cell, bool other) 
+void Cell::removeCoboundaryCell(Cell* cell, bool other) 
 {
   biter it = _cbd.find(cell);
   if(it != _cbd.end()){
     _cbd.erase(it);
     if(other) (*it).first->removeBoundaryCell(this, false);
-    return (*it).second;
   }
-  return 0;
 }
    
 bool Cell::hasBoundary(Cell* cell, bool orig)
@@ -329,7 +322,7 @@ CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell()
     Cell* cell = (*it).first;
     int ori = (*it).second;
     cell->removeCoboundaryCell(c1); 
-    if(this->addBoundaryCell(ori, cell)) cell->addCoboundaryCell(ori, this);
+    this->addBoundaryCell(ori, cell, false, true);
   }
   for(biter it = c2Boundary.begin(); it != c2Boundary.end(); it++){
     Cell* cell = (*it).first;
@@ -338,13 +331,11 @@ CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell()
     cell->removeCoboundaryCell(c2);    
     if(co){
       biter it2 = c1Boundary.find(cell);
-      if(it2 == c1Boundary.end() && this->addBoundaryCell(ori, cell)) { 
-	cell->addCoboundaryCell(ori, this);
+      if(it2 == c1Boundary.end()){
+	this->addBoundaryCell(ori, cell, false, true);
       }
     }
-    else{
-      if(this->addBoundaryCell(ori, cell)) cell->addCoboundaryCell(ori, this);
-    }
+    else this->addBoundaryCell(ori, cell, false, true);
   }
 
   // coboundary cells
@@ -357,7 +348,7 @@ CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell()
     Cell* cell = (*it).first;
     int ori = (*it).second;
     cell->removeBoundaryCell(c1); 
-    if(this->addCoboundaryCell(ori, cell)) cell->addBoundaryCell(ori, this);
+    this->addCoboundaryCell(ori, cell, false, true);
   }
   for(biter it = c2Coboundary.begin(); it != c2Coboundary.end(); it++){
     Cell* cell = (*it).first;
@@ -366,13 +357,11 @@ CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell()
     cell->removeBoundaryCell(c2);    
     if(!co){
       biter it2 = c1Coboundary.find(cell);
-      if(it2 == c1Coboundary.end() && this->addCoboundaryCell(ori, cell)){
-	cell->addBoundaryCell(ori, this); 
+      if(it2 == c1Coboundary.end()){
+	this->addCoboundaryCell(ori, cell, false, true);
       }
     }
-    else {
-      if(this->addCoboundaryCell(ori, cell)) cell->addBoundaryCell(ori, this);
-    }
+    else this->addCoboundaryCell(ori, cell, false, true);
   }
 
 }
diff --git a/Geo/Cell.h b/Geo/Cell.h
index d79a7e10a5870a5e0666a99c62496eda047ed2f5..27d2002834b35fe4f88b6664f70602042fb397c1 100644
--- a/Geo/Cell.h
+++ b/Geo/Cell.h
@@ -75,7 +75,7 @@ class Cell
    // get the number of vertices this cell has
   int getNumVertices() const { 
     if(_image == NULL) printf("ERROR: No image mesh element for cell. \n");
-    return _image->getNumVertices(); }
+    return _image->getNumPrimaryVertices(); }
   // get the number of facets of this cell
   int getNumFacets() const;
   // get the vertices on a facet of this cell
@@ -143,14 +143,14 @@ class Cell
     orig ? coboundary = _ocbd : coboundary = _cbd; }
   
   // add (co)boundary cell
-  virtual bool addBoundaryCell(int orientation, Cell* cell, 
-			       bool orig=false, bool other=false); 
-  virtual bool addCoboundaryCell(int orientation, Cell* cell, 
-				 bool orig=false, bool other=false);
+  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);
   
   // remove (co)boundary cell
-  virtual int removeBoundaryCell(Cell* cell, bool other=true);
-  virtual int removeCoboundaryCell(Cell* cell, bool other=true);
+  virtual void removeBoundaryCell(Cell* cell, bool other=true);
+  virtual void removeCoboundaryCell(Cell* cell, bool other=true);
   
   // true if has given cell on (original) (co)boundary
   virtual bool hasBoundary(Cell* cell, bool orig=false);
diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 2698f839f1edba8fbceffd4584b768115bfa404c..143032694b0b54a9049152d47bb2f9257605e3f5 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -81,8 +81,7 @@ bool CellComplex::insert_cells(std::vector<MElement*>& elements,
 	}
 	if(!subdomain) {
 	  int ori = cell->getFacetOri(newCell);
-	  cell->addBoundaryCell( ori, newCell, true);
-	  newCell->addCoboundaryCell( ori, cell, true);
+	  cell->addBoundaryCell( ori, newCell, true, true);
 	}
       }
     }
@@ -267,9 +266,6 @@ int CellComplex::reduceComplex(bool omit)
   int count = 0;
   for(int i = 3; i > 0; i--) count = count + reduction(i);
 
-  printf(" %d volumes, %d faces, %d edges and %d vertices. \n",
-         getSize(3), getSize(2), getSize(1), getSize(0));
-
   if(omit){
     int omitted = 0;
     _store.clear();
@@ -290,12 +286,15 @@ int CellComplex::reduceComplex(bool omit)
     }
   }
   
+  printf(" %d volumes, %d faces, %d edges and %d vertices. \n",
+       getSize(3), getSize(2), getSize(1), getSize(0));
+
   combine(3);
   reduction(2);
   combine(2);
   reduction(1);
   combine(1);
-
+  
   printf(" %d volumes, %d faces, %d edges and %d vertices. \n",
 	 getSize(3), getSize(2), getSize(1), getSize(0));
   
@@ -332,10 +331,6 @@ int CellComplex::coreduceComplex(bool omit)
     }
   } 
   
-  printf(" %d volumes, %d faces, %d edges and %d vertices. \n",
-	 getSize(3), getSize(2), getSize(1), getSize(0));
-
-  
   int omitted = 0;
   if(omit){
     _store.clear();
@@ -352,6 +347,9 @@ int CellComplex::coreduceComplex(bool omit)
       coreduction(cell, omitted);
     }
   }
+
+  printf(" %d volumes, %d faces, %d edges and %d vertices. \n",
+	 getSize(3), getSize(2), getSize(1), getSize(0));
   
   cocombine(0);
   coreduction(1);
@@ -508,7 +506,7 @@ bool CellComplex::coherent()
         }
         if(!bdCell->hasCoboundary(cell)){
           printf("Warning! Incoherent boundary/coboundary pair! Fixed. \n");
-	  bdCell->addCoboundaryCell(ori, cell);
+	  bdCell->addCoboundaryCell(ori, cell, false, false);
           coherent = false;
         }
         
@@ -527,7 +525,7 @@ bool CellComplex::coherent()
         }
         if(!cbdCell->hasBoundary(cell)){
           printf("Warning! Incoherent coboundary/boundary pair! Fixed. \n");
-	  cbdCell->addBoundaryCell(ori, cell);
+	  cbdCell->addBoundaryCell(ori, cell, false, false);
           coherent = false;
         }