Some more.

Geo/Cell.cpp

@@ -283,27 +283,21 @@ CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell() {
   std::map< Cell*, int, Less_Cell > c2Boundary; 
   c2->getBoundary(c2Boundary);
   
-  for(std::map<Cell*, int, Less_Cell>::iterator it = c1Boundary.begin();
-      it != c1Boundary.end(); it++){
+  for(biter it = c1Boundary.begin(); it != c1Boundary.end(); it++){
     Cell* cell = (*it).first;
     int ori = (*it).second;
     cell->removeCoboundaryCell(c1); 
     if(this->addBoundaryCell(ori, cell)) cell->addCoboundaryCell(ori, this);
   }
-  for(std::map<Cell*, int, Less_Cell >::iterator it = c2Boundary.begin(); 
-      it != c2Boundary.end(); it++){
+  for(biter it = c2Boundary.begin(); it != c2Boundary.end(); it++){
     Cell* cell = (*it).first;
     if(!orMatch) (*it).second = -1*(*it).second;
     int ori = (*it).second;
     cell->removeCoboundaryCell(c2);    
     if(co){
-      std::map<Cell*, int, Less_Cell >::iterator it2 = c1Boundary.find(cell);
-      bool old = false;
-      if(it2 != c1Boundary.end()) old = true;
-      if(!old){  
-	if(this->addBoundaryCell(ori, cell)) { 
-	  cell->addCoboundaryCell(ori, this);
-	}
+      biter it2 = c1Boundary.find(cell);
+      if(it2 == c1Boundary.end() && this->addBoundaryCell(ori, cell)) { 
+	cell->addCoboundaryCell(ori, this);
       }
     }
     else{
@@ -317,27 +311,21 @@ CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell() {
   std::map<Cell*, int, Less_Cell > c2Coboundary;
   c2->getCoboundary(c2Coboundary);
   
-  for(std::map<Cell*, int, Less_Cell>::iterator it = c1Coboundary.begin(); 
-      it != c1Coboundary.end(); it++){
+  for(biter it = c1Coboundary.begin(); it != c1Coboundary.end(); it++){
     Cell* cell = (*it).first;
     int ori = (*it).second;
     cell->removeBoundaryCell(c1); 
     if(this->addCoboundaryCell(ori, cell)) cell->addBoundaryCell(ori, this);
   }
-  for(std::map<Cell*, int, Less_Cell>::iterator it = c2Coboundary.begin();
-      it != c2Coboundary.end(); it++){
+  for(biter it = c2Coboundary.begin(); it != c2Coboundary.end(); it++){
     Cell* cell = (*it).first;
     if(!orMatch) (*it).second = -1*(*it).second;
     int ori = (*it).second;
     cell->removeBoundaryCell(c2);    
     if(!co){
-      std::map<Cell*, int, Less_Cell >::iterator it2 = c1Coboundary.find(cell);
-      bool old = false;
-      if(it2 != c1Coboundary.end()) old = true;
-      if(!old) { 
-	if(this->addCoboundaryCell(ori, cell)){
-	  cell->addBoundaryCell(ori, this); 
-	}
+      biter it2 = c1Coboundary.find(cell);
+      if(it2 == c1Coboundary.end() && this->addCoboundaryCell(ori, cell)){
+	cell->addBoundaryCell(ori, this); 
       }
     }
     else {

Geo/Cell.h

@@ -116,14 +116,18 @@ class Cell
   
   // (co)boundary cell iterator
   typedef std::map<Cell*, int, Less_Cell>::iterator biter;
-  
+  biter firstBoundary(){ return _boundary.begin(); }
+  biter lastBoundary(){ return _boundary.end(); }
+  biter firstCoboundary(){ return _coboundary.begin(); }
+  biter lastCoboundary(){ return _coboundary.end(); }
+
   virtual int getBoundarySize() { return _boundary.size(); }
   virtual int getCoboundarySize() { return _coboundary.size(); }
    
   // get the cell boundary
   virtual void getBoundary(std::map<Cell*, int, Less_Cell >& boundary){
     boundary =  _boundary; }
-  virtual void getCoboundary(  std::map<Cell*, int, Less_Cell >& coboundary ){
+  virtual void getCoboundary(std::map<Cell*, int, Less_Cell >& coboundary){
     coboundary = _coboundary; }
   
   // add (co)boundary cell

Geo/CellComplex.cpp

@@ -277,14 +277,12 @@ void CellComplex::removeCell(Cell* cell, bool other){
   std::map<Cell*, int, Less_Cell > boundary; 
   cell->getBoundary(boundary);
   
-  for(std::map<Cell*, int, Less_Cell>::iterator it = coboundary.begin();
-      it != coboundary.end(); it++){
+  for(Cell::biter it = coboundary.begin(); it != coboundary.end(); it++){
     Cell* cbdCell = (*it).first;
     cbdCell->removeBoundaryCell(cell, other);
   } 
   
-  for(std::map<Cell*, int, Less_Cell>::iterator it = boundary.begin();
-      it != boundary.end(); it++){
+  for(Cell::biter it = boundary.begin(); it != boundary.end(); it++){
     Cell* bdCell = (*it).first;
     bdCell->removeCoboundaryCell(cell, other);
   }
@@ -293,7 +291,7 @@ void CellComplex::removeCell(Cell* cell, bool other){
   
 }
 
-void CellComplex::removeCellQset(Cell*& cell, 
+void CellComplex::removeCellQset(Cell* cell, 
 				 std::set<Cell*, Less_Cell>& Qset){
   Qset.erase(cell);
 }
@@ -325,8 +323,7 @@ int CellComplex::coreduction(Cell* generator, int omitted){
   
   std::map<Cell*, int, Less_Cell > bd_s;
   std::map<Cell*, int, Less_Cell > cbd_c;
-  //std::list<Cell*> bd_s;
-  //std::list<Cell*> cbd_c;
+
   Cell* s;
   int round = 0;
   while( !Q.empty() ){
@@ -336,9 +333,9 @@ int CellComplex::coreduction(Cell* generator, int omitted){
     s = Q.front();
     Q.pop();
     removeCellQset(s, Qset);
-    s->getBoundary(bd_s);
-    
-    if( bd_s.size() == 1 && inSameDomain(s, bd_s.begin()->first) ){
+    if(s->getBoundarySize() == 1 
+       && inSameDomain(s, s->firstBoundary()->first) ){
+      s->getBoundary(bd_s);
       removeCell(s);
       bd_s.begin()->first->getCoboundary(cbd_c);
       enqueueCells(cbd_c, Q, Qset);
@@ -347,10 +344,8 @@ int CellComplex::coreduction(Cell* generator, int omitted){
 	_store.at(omitted-1).insert(bd_s.begin()->first);
       }
       coreductions++;
-      
-      
     }
-    else if(bd_s.empty()){
+    else if(s->getBoundarySize() == 0){
       s->getCoboundary(cbd_c);
       enqueueCells(cbd_c, Q, Qset);
     }
@@ -363,7 +358,6 @@ int CellComplex::coreduction(Cell* generator, int omitted){
 
 int CellComplex::reduction(int dim, int omitted){
   if(dim < 1 || dim > 3) return 0;
-  //std::list<Cell*> cbd_c;
   std::map<Cell*, int, Less_Cell > cbd_c;
   int count = 0;
   
@@ -377,28 +371,25 @@ int CellComplex::reduction(int dim, int omitted){
       
       
       Cell* cell = *cit;
-      cell->getCoboundary(cbd_c);
-      if( cbd_c.size() == 1 && inSameDomain(cell, cbd_c.begin()->first)){
-          //&& ( (!cell->getImmune() && !cbd_c.front()->getImmune() ) )){
-        ++cit;
-        removeCell(cbd_c.begin()->first);
-        removeCell(cell);
-        if(dim == getDim() && omitted > 0){
-	  _store.at(omitted-1).insert(cbd_c.begin()->first);    
+      if( cell->getCoboundarySize() == 1 
+	  && inSameDomain(cell, cell->firstCoboundary()->first)){
+	++cit;
+	if(dim == getDim() && omitted > 0){
+	  _store.at(omitted-1).insert(cell->firstCoboundary()->first);    
 	}
-        
-        count++;
-        reduced = true;
-        
+	removeCell(cell->firstCoboundary()->first);
+	removeCell(cell);
+	count++;
+	reduced = true;
+	
       }
+    
       if(getSize(dim) == 0 || getSize(dim-1) == 0) break;
       cit++;
-      
     }
-    //if(!reduced && ignoreCells) { ignoreCells = false; reduced = true;}
     
   }
-
+  
   return count;
 }
 /*
@@ -641,7 +632,6 @@ int CellComplex::cocombine(int dim){
   std::queue<Cell*> Q;
   std::set<Cell*, Less_Cell> Qset;
   std::map<Cell*, int, Less_Cell> cbd_c;
-  std::map<Cell*, int, Less_Cell > bd_c;
   int count = 0;
   
   for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
@@ -654,10 +644,8 @@ int CellComplex::cocombine(int dim){
       Cell* s = Q.front();
       Q.pop();
       
-      s->getBoundary(bd_c);
       if(s->getBoundarySize() == 2){
-      
-        std::map<Cell*, int, Less_Cell>::iterator it = bd_c.begin();
+	Cell::biter it = s->firstBoundary();
         int or1 = (*it).second;
         Cell* c1 = (*it).first;
         it++;
@@ -706,7 +694,7 @@ int CellComplex::combine(int dim){
   
   std::queue<Cell*> Q;
   std::set<Cell*, Less_Cell> Qset;
-  std::map<Cell*, int, Less_Cell> cbd_c;
+  //std::map<Cell*, int, Less_Cell> cbd_c;
   std::map<Cell*, int, Less_Cell> bd_c;
   int count = 0;
 
@@ -720,11 +708,9 @@ int CellComplex::combine(int dim){
       
       Cell* s = Q.front();
       Q.pop(); 
-      s->getCoboundary(cbd_c);
 
       if(s->getCoboundarySize() == 2){
-
-        std::map<Cell*, int, Less_Cell>::iterator it = cbd_c.begin();
+	Cell::biter it = s->firstCoboundary();
         int or1 = (*it).second;
         Cell* c1 = (*it).first;
         it++;

Geo/CellComplex.h

@@ -69,7 +69,7 @@ class CellComplex
   void enqueueCells(std::map<Cell*, int, Less_Cell>& cells, 
 		    std::queue<Cell*>& Q, std::set<Cell*, Less_Cell>& Qset);
   // remove cell from the queue set
-  void removeCellQset(Cell*& cell, std::set<Cell*, Less_Cell>& Qset);
+  void removeCellQset(Cell* cell, std::set<Cell*, Less_Cell>& Qset);
   
   // for constructor 
   void insert_cells(bool subdomain, bool boundary);

Geo/ChainComplex.cpp

@@ -76,10 +76,8 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){
 	   cit != cellComplex->lastCell(dim); cit++){
         Cell* cell = *cit;
         if(!cell->inSubdomain()){
-          std::map<Cell*, int, Less_Cell> bdCell;
-	  cell->getBoundary(bdCell);
-          for(std::map<Cell*, int, Less_Cell>::iterator it = bdCell.begin();
-	      it != bdCell.end(); it++){
+          for(Cell::biter it = cell->firstBoundary();
+	      it != cell->lastBoundary(); it++){
             Cell* bdCell = (*it).first;
             if(!bdCell->inSubdomain()){
               int old_elem = 0;