diff --git a/Geo/Cell.cpp b/Geo/Cell.cpp
index 012b6ea9a88a227d208f9af3984723360da66aba..9a64b308adfefaf2c8a7bfab2a88172dae60c658 100755
--- a/Geo/Cell.cpp
+++ b/Geo/Cell.cpp
@@ -1,23 +1,11 @@
-// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle
-//
-// See the LICENSE.txt file for license information. Please report all
-// bugs and problems to <gmsh@geuz.org>.
-//
-// Contributed by Matti Pellikka <matti.pellikka@tut.fi>.
-
-#include "CellComplex.h"
 #include "Cell.h"
 
-#if defined(HAVE_KBIPACK)
-
 bool Less_Cell::operator()(const Cell* c1, const Cell* c2) const 
 {  
-  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 cell complex is done use enumeration (same as vertex order)
+  if(c1->getNum() != 0) return (c1->getNum() < c2->getNum());
   
+  // Otherwise order by vertex numbering (good heuristic for reduction)
   if(c1->getNumSortedVertices() != c2->getNumSortedVertices()){
     return (c1->getNumSortedVertices() < c2->getNumSortedVertices());
   }    
@@ -28,120 +16,214 @@ bool Less_Cell::operator()(const Cell* c1, const Cell* c2) const
   return false;
 }
 
-Cell::Cell(MElement* image) :  
-  _subdomain(false), _combined(false), _index(0), _immune(false), _image(NULL), 
-  _delimage(false)
-{
-  _image = image;
-  for(int i = 0; i < getNumVertices(); i++) 
-    _vs.push_back(getVertex(i)->getNum()); 
-  std::sort(_vs.begin(), _vs.end());
+int Cell::_globalNum = 0;
+Cell::Cell(int dim, std::vector<MVertex*>& v) 
+{  
+  _pnum = 0;
+  _dim = dim;
+  _domain = 0;
+  _combined = false;
+  _immune = false;
+  _v = v;
+  std::sort(_v.begin(), _v.end(), MVertexLessThanNum());
+
+  //_num = ++_globalNum;
+  _num = 0;
+  _index = 0;
 }
+/*Cell::Cell(int num, int dim, int type, std::vector<int>& v, int domain) 
+{  
+  _num = num;
+  _dim = dim;
+  _domain = domain;
+  _combined = false;
+  _immune = false;
+  _type = type;
+  _v = v;
+  _index = 0;
+  }*/
 
-Cell::~Cell() 
-{
-  if(_delimage) delete _image; 
+Cell::Cell(MElement* element, int domain)
+{  
+  _pnum = element->getNum();
+  _dim = element->getDim();
+  _domain = domain;
+  _combined = false;
+  _immune = false;
+  int type = element->getType();
+  if(type != TYPE_PNT && type != TYPE_LIN && type != TYPE_TRI && type != TYPE_TET) _type = 2;
+  else _type = 1;
+
+  for(int i = 0; i < element->getNumPrimaryVertices(); i++){
+    _v.push_back(element->getVertex(i));
+  }
+  std::sort(_v.begin(), _v.end(), MVertexLessThanNum());
+
+  //_num = ++_globalNum;
+  _num = 0;
+  _index = 0;
+  
+}
+
+Cell::Cell(Cell* parent, int i)
+{  
+  _pnum = parent->getParentNum();
+  _dim = parent->getDim()-1;
+  _domain = parent->getDomain();
+  _combined = false;
+  _immune = false;
+
+  parent->findBdElement(i, _type, _v);
+  std::sort(_v.begin(), _v.end(), MVertexLessThanNum());
+
+  //_num = ++_globalNum;
+  _num = 0;
+  _index = 0;
+  
 }
 
-bool Cell::findBoundaryElements(std::vector<MElement*>& bdElements)
+void Cell::findBdElement(int i, int& type, std::vector<MVertex*>& vertices) const 
 {
-  if(_combined) return false;
-  bdElements.clear();
-  MElementFactory factory;
-  for(int i = 0; i < getNumFacets(); i++){
-    std::vector<MVertex*> vertices;
-    getFacetVertices(i, vertices);
-    int type = _image->getType();
-    int newtype = 0;
-    if(getDim() == 3){
-      if(type == TYPE_TET) newtype = MSH_TRI_3;
-      else if(type == TYPE_HEX) newtype = MSH_QUA_4;
-      else if(type == TYPE_PRI) {
-	if(vertices.size() == 3) newtype = MSH_TRI_3;
-	else if(vertices.size() == 4) newtype = MSH_QUA_4;
+  vertices.clear();
+  type = 0;
+  if(_type == 1){ 
+    if(_dim == 1){ // i = 0 -> ori = -1, i = 1 -> ori = 1
+      type = 1;
+      vertices.push_back(_v[i]);
+    }
+    else if(_dim == 2){
+      type = 1;
+      if(i == 0){ // ori = 1
+	vertices.push_back(_v[0]);
+	vertices.push_back(_v[1]);
+      }
+      else if(i == 1){ // ori = 1
+	vertices.push_back(_v[1]);
+	vertices.push_back(_v[2]);
+      }
+      else if(i == 2){ // ori = 1
+	vertices.push_back(_v[2]);
+	vertices.push_back(_v[0]);
       }
     }
-    else if(getDim() == 2) newtype = MSH_LIN_2;
-    else if(getDim() == 1) newtype = MSH_PNT;
-    if(newtype == 0){
-      printf("Error: mesh element %d not implemented yet! \n", type);
-      return false;
+    else if(_dim == 3){
+      type = 1;
+      if(i == 0){
+	vertices.push_back(_v[0]);
+	vertices.push_back(_v[2]);
+	vertices.push_back(_v[1]);
+      }
+      else if(i == 1){
+	vertices.push_back(_v[0]);
+	vertices.push_back(_v[1]);
+	vertices.push_back(_v[3]);
+      }
+      else if(i == 2){
+	vertices.push_back(_v[0]);
+	vertices.push_back(_v[3]);
+	vertices.push_back(_v[2]);
+      }
+      else if(i == 3){
+	vertices.push_back(_v[3]);
+	vertices.push_back(_v[1]);
+	vertices.push_back(_v[2]);
+      }
     }
-    MElement* element = factory.create(newtype, vertices, 0, 
-				       _image->getPartition());
-    bdElements.push_back(element);
   }
-  return true;
 }
-
-int Cell::getNumFacets() const 
-{
-  if(_image == NULL || _combined){ 
-    printf("ERROR: No image mesh element for cell. \n");
-    return 0;
-  }
-  if(getDim() == 0) return 0;
-  else if(getDim() == 1) return 2;
-  else if(getDim() == 2) return _image->getNumEdges();
-  else if(getDim() == 3) return _image->getNumFaces();
+int Cell::getNumBdElements() const 
+{ 
+  if(_dim == 0) return 0;
+  else if(_dim == 1) return 2;
+  else if(_dim == 2 && _type == 1) return 3;
+  else if(_dim == 3 && _type == 1) return 4;
   else return 0;
 }
+int Cell::findBdCellOrientation(Cell* cell) const 
+{ 
+  std::vector<MVertex*> v;
+  cell->getMeshVertices(v);
+  if(_dim == 0) return 0;
+  else if(_dim == 1){
+    if(v[0]->getNum() == _v[0]->getNum()) return -1;
+    else if(v[0]->getNum() == _v[1]->getNum()) return 1;
+  }
+  else if(_type == 1){
+    if(_dim == 2){
+      if     (v[0]->getNum() == _v[0]->getNum() && v[1]->getNum() == _v[1]->getNum() ) return 1;
+      else if(v[0]->getNum() == _v[1]->getNum() && v[1]->getNum() == _v[0]->getNum() ) return -1;
+      
+      else if(v[0]->getNum() == _v[1]->getNum() && v[1]->getNum() == _v[2]->getNum() ) return 1;
+      else if(v[0]->getNum() == _v[2]->getNum() && v[1]->getNum() == _v[1]->getNum() ) return -1;
+      
+      else if(v[0]->getNum() == _v[2]->getNum() && v[1]->getNum() == _v[0]->getNum() ) return 1;
+      else if(v[0]->getNum() == _v[0]->getNum() && v[1]->getNum() == _v[2]->getNum() ) return -1;
+    } 
+    else if(_dim == 3) {
+      if     (v[0]->getNum() == _v[0]->getNum() && v[1]->getNum() == _v[2]->getNum() && v[2]->getNum() == _v[1]->getNum() ) return 1;
+      else if(v[0]->getNum() == _v[0]->getNum() && v[1]->getNum() == _v[1]->getNum() && v[2]->getNum() == _v[2]->getNum()) return -1;
+      
+      else if(v[0]->getNum() == _v[2]->getNum() && v[1]->getNum() == _v[1]->getNum() && v[2]->getNum() == _v[0]->getNum()) return 1;    
+      else if(v[0]->getNum() == _v[2]->getNum() && v[1]->getNum() == _v[0]->getNum() && v[2]->getNum() == _v[1]->getNum()) return -1;
+      
+      else if(v[0]->getNum() == _v[1]->getNum() && v[1]->getNum() == _v[0]->getNum() && v[2]->getNum() == _v[2]->getNum()) return 1;    
+      else if(v[0]->getNum() == _v[1]->getNum() && v[1]->getNum() == _v[2]->getNum() && v[2]->getNum() == _v[0]->getNum()) return -1;
+      
+      
+      else if(v[0]->getNum() == _v[0]->getNum() && v[1]->getNum() == _v[1]->getNum() && v[2]->getNum() == _v[3]->getNum() ) return 1;
+      else if(v[0]->getNum() == _v[0]->getNum() && v[1]->getNum() == _v[3]->getNum() && v[2]->getNum() == _v[1]->getNum()) return -1;
+      
+      else if(v[0]->getNum() == _v[1]->getNum() && v[1]->getNum() == _v[3]->getNum() && v[2]->getNum() == _v[0]->getNum()) return 1;    
+      else if(v[0]->getNum() == _v[1]->getNum() && v[1]->getNum() == _v[0]->getNum() && v[2]->getNum() == _v[3]->getNum()) return -1;
+      
+      else if(v[0]->getNum() == _v[3]->getNum() && v[1]->getNum() == _v[0]->getNum() && v[2]->getNum() == _v[1]->getNum()) return 1;    
+      else if(v[0]->getNum() == _v[3]->getNum() && v[1]->getNum() == _v[1]->getNum() && v[2]->getNum() == _v[0]->getNum()) return -1;
+      
+      
+      else if(v[0]->getNum() == _v[0]->getNum() && v[1]->getNum() == _v[3]->getNum() && v[2]->getNum() == _v[2]->getNum() ) return 1;
+      else if(v[0]->getNum() == _v[0]->getNum() && v[1]->getNum() == _v[2]->getNum() && v[2]->getNum() == _v[3]->getNum()) return -1;
+      
+      else if(v[0]->getNum() == _v[3]->getNum() && v[1]->getNum() == _v[2]->getNum() && v[2]->getNum() == _v[3]->getNum()) return 1;    
+      else if(v[0]->getNum() == _v[3]->getNum() && v[1]->getNum() == _v[3]->getNum() && v[2]->getNum() == _v[2]->getNum()) return -1;
+      
+      else if(v[0]->getNum() == _v[2]->getNum() && v[1]->getNum() == _v[0]->getNum() && v[2]->getNum() == _v[3]->getNum()) return 1;    
+      else if(v[0]->getNum() == _v[2]->getNum() && v[1]->getNum() == _v[3]->getNum() && v[2]->getNum() == _v[0]->getNum()) return -1;
+      
 
-void Cell::getFacetVertices(const int num, std::vector<MVertex*> &v) const 
-{
-  if(_image == NULL || _combined){ 
-    printf("ERROR: No image mesh element for cell. \n");
-    return;
+      else if(v[0]->getNum() == _v[3]->getNum() && v[1]->getNum() == _v[1]->getNum() && v[2]->getNum() == _v[2]->getNum() ) return 1;
+      else if(v[0]->getNum() == _v[3]->getNum() && v[1]->getNum() == _v[2]->getNum() && v[2]->getNum() == _v[1]->getNum()) return -1;
+      
+      else if(v[0]->getNum() == _v[1]->getNum() && v[1]->getNum() == _v[2]->getNum() && v[2]->getNum() == _v[3]->getNum()) return 1;    
+      else if(v[0]->getNum() == _v[1]->getNum() && v[1]->getNum() == _v[3]->getNum() && v[2]->getNum() == _v[2]->getNum()) return -1;
+      
+      else if(v[0]->getNum() == _v[2]->getNum() && v[1]->getNum() == _v[3]->getNum() && v[2]->getNum() == _v[1]->getNum()) return 1;    
+      else if(v[0]->getNum() == _v[2]->getNum() && v[1]->getNum() == _v[1]->getNum() && v[2]->getNum() == _v[3]->getNum()) return -1;
+    }
   }
-  if(getDim() == 0) return;
-  else if(getDim() == 1) { v.resize(1); v[0] = getVertex(num); }
-  else if(getDim() == 2) _image->getEdgeVertices(num, v);
-  else if(getDim() == 3) _image->getFaceVertices(num, v);
-  return;
+  return 0;
 }
 
-int Cell::findBoundaryCellOrientation(Cell* cell) 
+int Cell::getTypeMSH() const 
 {
-  if(_image == NULL || _combined){ 
-    printf("ERROR: No image mesh element for cell. \n");
-    return 0;
-  }
-  std::vector<MVertex*> v; 
-  for(int i = 0; i < cell->getNumVertices(); i++) {
-    v.push_back(cell->getVertex(i));
-  }
-  if(getDim() == 0) return 0;
-  else if(getDim() == 1){
-    if(v.size() != 1) return 0;
-    else if(v[0] == getVertex(0)) return -1;
-    else if(v[0] == getVertex(1)) return 1;
-    else return 0;
-  }
-  else if(getDim() == 2){
-    if(v.size() != 2) return 0;
-    MEdge facet = MEdge(v[0], v[1]);
-    int ithFacet = 0;
-    int sign = 0;
-    _image->getEdgeInfo(facet, ithFacet, sign);
-    return sign;
-  }
-  else if(getDim() == 3){
-    if(v.size() != 3) return 0;
-    MFace facet = MFace(v);
-    int ithFacet = 0;
-    int sign = 0;
-    int rot = 0;
-    _image->getFaceInfo(facet, ithFacet, sign, rot);
-    return sign;
+  if(_type == 1){
+    if(_dim == 0) return 15;
+    else if(_dim == 1) return 1;
+    else if(_dim == 2) return 2;
+    else if(_dim == 3) return 4;
   }
   else return 0;
-}  
+}
 
 bool Cell::hasVertex(int vertex) const 
 {
-  std::vector<int>::const_iterator it = std::find(_vs.begin(), _vs.end(), 
+  //_vs
+  std::vector<int> v;
+  for(unsigned int i = 0; i < _v.size(); i++) {
+    v.push_back(_v[i]->getNum());
+  }
+  std::vector<int>::const_iterator it = std::find(v.begin(), v.end(), 
 						  vertex);
-  if (it != _vs.end()) return true;
+  if (it != v.end()) return true;
   else return false;
 }
 
@@ -156,65 +238,76 @@ bool CombinedCell::hasVertex(int vertex) const
 
 void Cell::printCell() 
 {
-  printf("%d-cell: \n" , getDim());
+  printf("%d-cell %d: \n" , getDim(), getNum());
   printf("Vertices: ");
-  for(int i = 0; i < this->getNumSortedVertices(); i++){
-    printf("%d ", this->getSortedVertex(i));
+  for(int i = 0; i < this->getNumVertices(); i++){
+    printf("%d ", this->getVertex(i));
   }
   printf(", in subdomain: %d, ", inSubdomain());
   printf("combined: %d. \n" , isCombined() );
 };
 
 void Cell::restoreCell(){
-  _bd = _obd;
-  _cbd = _ocbd;
+  std::vector<Cell*> toRemove;
+  for(biter it = firstCoboundary(true); it != lastCoboundary(); it++){
+    it->second.reset();
+    if(it->second.get() == 0) toRemove.push_back(it->first);
+  }
+  for(unsigned int i = 0; i < toRemove.size(); i++) _cbd.erase(toRemove[i]);
+  toRemove.clear();
+  for(biter it = firstBoundary(true); it != lastBoundary(); it++){
+    it->second.reset();
+    if(it->second.get() == 0) toRemove.push_back(it->first);
+  }
+  for(unsigned int i = 0; i < toRemove.size(); i++) _bd.erase(toRemove[i]);
   _combined = false;
   _index = 0;
   _immune = false;   
 }
 
-void Cell::addBoundaryCell(int orientation, Cell* cell, 
-			   bool orig, bool other) 
+void Cell::addBoundaryCell(int orientation, Cell* cell, bool other) 
 {
-  if(orig) _obd.insert( std::make_pair(cell, orientation ) );
   biter it = _bd.find(cell);
   if(it != _bd.end()){
-    int newOrientation = (*it).second + orientation;
-    if(newOrientation != 0) (*it).second = newOrientation;
-    else {
-      _bd.erase(it);
-      (*it).first->removeCoboundaryCell(this,false);
+    int newOrientation = it->second.get() + orientation;
+    it->second.set(newOrientation);
+    if(newOrientation == 0){
+      it->first->removeCoboundaryCell(this, false);
+      if(it->second.geto() == 0) {
+	_bd.erase(it);
+      }
       return;
     }
   }
-  else _bd.insert( std::make_pair(cell, orientation ) );
-  if(other) cell->addCoboundaryCell(orientation, this, orig, false);
+  else _bd.insert( std::make_pair(cell, BdInfo(orientation) ) );
+  if(other) cell->addCoboundaryCell(orientation, this, false);
 }
 
-void Cell::addCoboundaryCell(int orientation, Cell* cell, 
-			     bool orig, bool other) 
+void Cell::addCoboundaryCell(int orientation, Cell* cell, bool other) 
 {
-  if(orig) _ocbd.insert( std::make_pair(cell, orientation ) );
   biter it = _cbd.find(cell);
   if(it != _cbd.end()){
-    int newOrientation = (*it).second + orientation;
-    if(newOrientation != 0) (*it).second = newOrientation;
-    else {
-      _cbd.erase(it);
-      (*it).first->removeBoundaryCell(this,false);
+    int newOrientation = it->second.get() + orientation;
+    it->second.set(newOrientation);
+    if(newOrientation == 0) {
+      it->first->removeBoundaryCell(this, false);
+      if(it->second.geto() == 0) {
+	_cbd.erase(it);
+      }
       return;
     }
   }
-  else _cbd.insert( std::make_pair(cell, orientation ) );
-  if(other) cell->addBoundaryCell(orientation, this, orig, false);
+  else _cbd.insert( std::make_pair(cell, BdInfo(orientation) ) );
+  if(other) cell->addBoundaryCell(orientation, this, false);
 }
 
 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);
+    it->second.set(0);
+    if(it->second.geto() == 0) _bd.erase(it);
+    if(other) it->first->removeCoboundaryCell(this, false);
   }
 }
  
@@ -222,8 +315,9 @@ 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);
+    it->second.set(0);
+    if(it->second.geto() == 0) _cbd.erase(it);
+    if(other) it->first->removeBoundaryCell(this, false);
   }
 }
    
@@ -231,12 +325,12 @@ bool Cell::hasBoundary(Cell* cell, bool orig)
 {
   if(!orig){
     biter it = _bd.find(cell);
-    if(it != _bd.end()) return true;
+    if(it != _bd.end() && it->second.get() != 0) return true;
     return false;
   }
   else{
-    biter it = _obd.find(cell);
-    if(it != _obd.end()) return true;
+    biter it = _bd.find(cell);
+    if(it != _bd.end() && it->second.geto() != 0) return true;
     return false;
   }
 }
@@ -245,16 +339,16 @@ bool Cell::hasCoboundary(Cell* cell, bool orig)
 {
   if(!orig){
     biter it = _cbd.find(cell);
-    if(it != _cbd.end()) return true;
+    if(it != _cbd.end() && it->second.get() != 0) return true;
     return false;
   }
   else{
-    biter it = _ocbd.find(cell);
-    if(it != _ocbd.end()) return true;
+    biter it = _cbd.find(cell);
+    if(it != _cbd.end() && it->second.geto() != 0) return true;
     return false;
   } 
 }
-
+/*
 void Cell::printBoundary(bool orig) 
 {  
   for(biter it = firstBoundary(orig); it != lastBoundary(orig); it++){
@@ -277,9 +371,8 @@ void Cell::printCoboundary(bool orig)
       printf("Cell coboundary is empty. \n");
     }
   }
-}
+  }*/
 
-int CombinedCell::_globalNum = 0;
 CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell() 
 {  
   // use "smaller" cell as c2
@@ -291,8 +384,9 @@ CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell()
   
   _num = ++_globalNum;
   _index = c1->getIndex();
-  _subdomain = c1->inSubdomain();
+  _domain = c1->getDomain();
   _combined = true;
+  _immune = false;
 
   // cells
   c1->getCells(_cells);
@@ -305,40 +399,44 @@ CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell()
 
   // boundary cells
   for(biter it = c1->firstBoundary(); it != c1->lastBoundary(); it++){
-    Cell* cell = (*it).first;
-    int ori = (*it).second;
+    Cell* cell = it->first;
+    int ori = it->second.get();
+    if(ori == 0) continue;
     cell->removeCoboundaryCell(c1, false); 
-    this->addBoundaryCell(ori, cell, false, true);
+    this->addBoundaryCell(ori, cell, true);
   }
   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* cell = it->first;
+    if(!orMatch) it->second.set(-1*it->second.get());
+    int ori = it->second.get();
+    if(ori == 0) continue;
     cell->removeCoboundaryCell(c2, false);    
     if(co && !c1->hasBoundary(cell)){
-      this->addBoundaryCell(ori, cell, false, true);
+      this->addBoundaryCell(ori, cell, true);
     }
-    else if(!co) this->addBoundaryCell(ori, cell, false, true);
+    else if(!co) this->addBoundaryCell(ori, cell, true);
   }
   c1->clearBoundary();
   c2->clearBoundary();
 
   // coboundary cells
   for(biter it = c1->firstCoboundary(); it != c1->lastCoboundary(); it++){
-    Cell* cell = (*it).first;
-    int ori = (*it).second;
+    Cell* cell = it->first;
+    int ori = it->second.get();
+    if(ori == 0) continue;
     cell->removeBoundaryCell(c1, false); 
-    this->addCoboundaryCell(ori, cell, false, true);
+    this->addCoboundaryCell(ori, cell, true);
   }
   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* cell = it->first;
+    if(!orMatch) it->second.set(-1*it->second.get());
+    int ori = it->second.get();
+    if(ori == 0) continue;
     cell->removeBoundaryCell(c2, false);    
     if(!co && !c1->hasCoboundary(cell)){
-      this->addCoboundaryCell(ori, cell, false, true);
+      this->addCoboundaryCell(ori, cell, true);
     }
-    else if(co) this->addCoboundaryCell(ori, cell, false, true);
+    else if(co) this->addCoboundaryCell(ori, cell, true);
   }
   c1->clearCoboundary();
   c2->clearCoboundary();
@@ -350,7 +448,7 @@ CombinedCell::CombinedCell(std::vector<Cell*>& cells) : Cell()
 {  
   _num = ++_globalNum;
   _index = cells.at(0)->getIndex();
-  _subdomain = cells.at(0)->inSubdomain();
+  _domain = cells.at(0)->getDomain();
   _combined = true;
 
   // cells
@@ -359,6 +457,3 @@ CombinedCell::CombinedCell(std::vector<Cell*>& cells) : Cell()
     _cells[c] = 1;
   }
 }
-
-
-#endif
diff --git a/Geo/Cell.h b/Geo/Cell.h
index b4d70205a16214a6de36c95145256936dd3e1632..b9b9e22c5521baab41948baa375634befaab641a 100644
--- a/Geo/Cell.h
+++ b/Geo/Cell.h
@@ -1,118 +1,100 @@
-// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle
-//
-// See the LICENSE.txt file for license information. Please report all
-// bugs and problems to <gmsh@geuz.org>.
-//
-// Contributed by Matti Pellikka <matti.pellikka@tut.fi>.
-
 #ifndef _CELL_H_
 #define _CELL_H_
 
-#include "GmshConfig.h"
-
-#if defined(HAVE_KBIPACK)
-
-#include <stdio.h>
-#include <string>
-#include <algorithm>
-#include <set>
-#include <list>
 #include <map>
-#include "CellComplex.h"
+#include <vector>
 #include "MElement.h"
-#include "MVertex.h"
 
 class Cell;
 
-// Ordering of the cells
-class Less_Cell{
-  public:
-   bool operator()(const Cell* c1, const Cell* c2) const;
+class Less_Cell {
+public:
+  bool operator()(const Cell* c1, const Cell* c2) const;
+};
+
+// Class to save cell boundary orientation information
+class BdInfo {
+ private:
+  short int _ori;
+  short int _origOri;
+
+ public:
+  BdInfo(int ori) { _ori = ori; _origOri = 0; }
+  
+  int get() const { return _ori; }
+  void reset() { _ori = _origOri; }
+  void init() { _origOri = _ori; }
+  void set(int ori) { _ori = ori; } 
+  int geto() const { return _origOri; }
+
 };
 
+
 // Class representing an elementary cell of a cell complex.
-class Cell
-{  
+class Cell {
+
  protected:
+  static int _globalNum;
   
-  // whether this cell belongs to a subdomain
-  // used in relative homology computation
-  bool _subdomain;
-  
+  int _num;
+  // mutable index for each cell (used to create boundary operator matrices)
+  int _index;
+
+  int _domain;
+
   // whether this cell a combinded cell of elemetary cells
   bool _combined; 
-  
-  // mutable index for each cell (used to create boundary operator matrices)
-  int _index; 
-  
   // for some algorithms to omit this cell
   bool _immune;
-  
-  // mutable list of cells on the boundary and on the coboundary of this cell
-  std::map< Cell*, int, Less_Cell > _bd;
-  std::map< Cell*, int, Less_Cell > _cbd;
-  // immutable original boundary and coboundary before the reduction of the
-  // cell complex
-  std::map<Cell*, int, Less_Cell > _obd;
-  std::map<Cell*, int, Less_Cell > _ocbd;
 
- private:
+  // list of cells on the boundary and on the coboundary of this cell
+  std::map<Cell*, BdInfo, Less_Cell> _bd;
+  std::map<Cell*, BdInfo, Less_Cell> _cbd;
   
+ private:
+
+  int _dim;
+  int _pnum;
+  int _type;
   // sorted vertices of this cell (used for ordering of the cells)
-  std::vector<int> _vs;
-  
-  // the mesh element that is the image of this cell
-  MElement* _image;
-  // whether to delete the mesh element when done
-  // (set to true if created for homology computation only)
-  bool _delimage;
-  
-  // get vertex 
-  MVertex* getVertex(int vertex) const { return _image->getVertex(vertex); }  
-   // get the number of vertices this cell has
-  int getNumVertices() const { return _image->getNumPrimaryVertices(); }
-  // get the number of facets of this cell
-  int getNumFacets() const;
-  // get the vertices on a facet of this cell
-  void getFacetVertices(const int num, std::vector<MVertex*> &v) const;
-  
+  std::vector<MVertex*> _v;
+
  public:
 
-  Cell() : _subdomain(false), _combined(false), _index(0), _immune(false), _image(NULL), 
-    _delimage(false) {}
-  Cell(MElement* image);
-  ~Cell();
-
-  // the mesh element this cell is represented by
-  MElement* getImageMElement() const { 
-    if(_image == NULL || _combined){
-      printf("ERROR: No image mesh element for cell. \n"); }
-    return _image; }
-  // get the cell dimension
-  virtual int getDim() const { return _image->getDim(); };
-  // get/set whether the image mesh element should be deleted when
-  // this cell gets deleted 
-  // (was the mesh element in the original mesh,
-  // is it needed after homology computation for visualization?) 
-  void setDeleteImage(bool delimage) { if(!_combined) _delimage = delimage; }
-  bool getDeleteImage() const { 
-    if(!_combined) return _delimage;
-    else return false; }
-
-  // find the cells on the boundary of this cell
-  bool findBoundaryElements(std::vector<MElement*>& bdElements);
-  // get boundary cell orientation
-  int findBoundaryCellOrientation(Cell* cell);
-  
+  Cell(int dim, std::vector<MVertex*>& v); 
+ Cell() : _num(0), _dim(0), _index(0), _type(0), _domain(0), _combined(false),  _immune(false) {}
+
+  Cell(MElement* element, int domain);
+  Cell(Cell* parent, int i);
+
+  int getDomain() const { return _domain; }
+  int getNum() const { return _num; }
+  void setNum(int num) { _num = num; };
+  int getType() const { return _type; }
+  int getTypeMSH() const; 
+  virtual int getDim() const { return _dim; }
+  int getParentNum() const { return _pnum; }
+  bool inSubdomain() const { if(_domain != 0) return true; else return false; }
+  void getMeshVertices(std::vector<MVertex*>& v) const { v = _v; }
+
   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; }
+
+  int getNumSortedVertices() const { return _v.size(); }
+  int getSortedVertex(int vertex) const { return _v.at(vertex)->getNum(); }
+  int getNumVertices() const { return _v.size(); }
+  int getVertex(int vertex) const { return _v.at(vertex)->getNum(); }
+  MVertex* getMeshVertex(int vertex) const { return _v.at(vertex); }
+  void clearSortedVertices() { _v.clear(); }
+
+  void findBdElement(int i, int& type, std::vector<MVertex*>& vertices) const;
+  int getNumBdElements() const;
+  int findBdCellOrientation(Cell* cell) const;
+
+  void revertGlobalNum() { _globalNum--; }
+  void increaseGlobalNum() { _globalNum++; }
 
   // restores the cell information to its original state before reduction
   void restoreCell();
@@ -121,36 +103,73 @@ class Cell
   virtual bool hasVertex(int vertex) const;
   
   // (co)boundary cell iterator
-  typedef std::map<Cell*, int, Less_Cell>::iterator biter;
+  typedef std::map<Cell*, BdInfo, Less_Cell>::iterator biter;
   // iterators to (first/last (co)boundary cells of this cell
   // (orig: to original (co)boundary cells of this cell)
   biter firstBoundary(bool orig=false){ 
-    return orig ? _obd.begin() : _bd.begin(); }
-  biter lastBoundary(bool orig=false){ 
-    return orig ? _obd.end() : _bd.end(); }
+    biter it = _bd.begin();
+    if(!orig) while(it->second.get() == 0 && it != _bd.end()) it++;
+    else while(it->second.geto() == 0 && it != _bd.end()) it++;
+    return it;
+  }
+  biter lastBoundary(){ return _bd.end(); }
   biter firstCoboundary(bool orig=false){ 
-    return orig ? _ocbd.begin() : _cbd.begin(); }
-  biter lastCoboundary(bool orig=false){ 
-    return orig ? _ocbd.end() : _cbd.end(); }
+    biter it = _cbd.begin();
+    if(!orig) while(it->second.get() == 0 && it != _cbd.end()) it++;
+    else while(it->second.geto() == 0 && it != _cbd.end()) it++;
+    return it; 
+  }
+  biter lastCoboundary(){ return _cbd.end(); }
 
-  int getBoundarySize() { return _bd.size(); }
-  int getCoboundarySize() { return _cbd.size(); }
+  int getBoundarySize(bool orig=false) { 
+    int size = 0;
+    for(biter bit = _bd.begin(); bit != _bd.end(); bit++){
+      if(!orig && bit->second.get() != 0) size++;
+      else if(orig && bit->second.geto() != 0) size++;
+    }
+    return size; }
+  int getCoboundarySize(bool orig=false) { 
+    int size = 0;
+    for(biter bit = _cbd.begin(); bit != _cbd.end(); bit++){
+      if(!orig && bit->second.get() != 0) size++;
+      else if(orig && bit->second.geto() != 0) size++;
+    }
+    return size; }
    
   // get the (orig: original) cell boundary
-  void getBoundary(std::map<Cell*, int, Less_Cell >& boundary, 
-		   bool orig=false){
-    orig ? boundary = _obd : boundary =  _bd; }
-  void getCoboundary(std::map<Cell*, int, Less_Cell >& coboundary,
-		     bool orig=false){
-    orig ? coboundary = _ocbd : coboundary = _cbd; }
+  void getBoundary(std::map<Cell*, short int, Less_Cell >& boundary, bool orig=false){
+    boundary.clear();
+    for(biter it = firstBoundary(); it != lastBoundary(); it++){
+      Cell* cell = it->first;
+      if(!orig && it->second.get() != 0) boundary[cell] = it->second.get();
+      if(orig && it->second.geto() != 0) boundary[cell] = it->second.geto();
+    }
+  }
+  void getCoboundary(std::map<Cell*, short int, Less_Cell >& coboundary, bool orig = false){
+    coboundary.clear();
+    for(biter it = firstCoboundary(); it != lastCoboundary(); it++){
+      Cell* cell = it->first;
+      if(!orig && it->second.get() != 0) coboundary[cell] = it->second.get();
+      if(orig && it->second.geto() != 0) coboundary[cell] = it->second.geto();
+    }
+  }
+
   
-  // add (co)boundary cell (orig: as original)
+  // add (co)boundary cell
   // (other: reciprocally also add this cell from the other cell's (co)boundary)
-  void addBoundaryCell(int orientation, Cell* cell, 
-		       bool orig, bool other); 
-  void addCoboundaryCell(int orientation, Cell* cell, 
-			 bool orig, bool other);
+  void addBoundaryCell(int orientation, Cell* cell, bool other); 
+  void addCoboundaryCell(int orientation, Cell* cell, bool other);
   
+  //void saveOriginalBd() { _obd = _bd; _ocbd = _cbd; }
+  void saveOriginalBd() { 
+    for(biter it = firstCoboundary(); it != lastCoboundary(); it++){
+      it->second.init();
+    }
+    for(biter it = firstBoundary(); it != lastBoundary(); it++){
+      it->second.init();
+    }
+  }
+
   // remove (co)boundary cell 
   // (other: reciprocally also revove this cell from the other cell's (co)boundary)
   void removeBoundaryCell(Cell* cell, bool other);
@@ -160,13 +179,13 @@ class Cell
   bool hasBoundary(Cell* cell, bool orig=false);
   bool hasCoboundary(Cell* cell, bool orig=false);
   
-  void clearBoundary() { _bd.clear(); }
-  void clearCoboundary() { _cbd.clear(); }
+  void clearBoundary() { }//_bd.clear(); }
+  void clearCoboundary() {}// _cbd.clear(); }
   
   // print cell debug info
   virtual void printCell();
-  virtual void printBoundary(bool orig=false);
-  virtual void printCoboundary(bool orig=false);   
+  //virtual void printBoundary(bool orig=false);
+  //virtual void printCoboundary(bool orig=false);   
   
   // tools for combined cells
   bool isCombined() const { return _combined; }
@@ -178,22 +197,14 @@ class Cell
   virtual int getNumCells() const {return 1;}
   typedef std::map<Cell*, int, Less_Cell>::iterator citer;
 
-  // equivalence
-  virtual bool operator==(const Cell& c2) const {  
-    if(this->isCombined() != c2.isCombined()) return false;
-    
-    if(this->getNumSortedVertices() != c2.getNumSortedVertices()){
-      return false;
-    }
-    for(int i=0; i < this->getNumSortedVertices();i++){
-      if(this->getSortedVertex(i) != c2.getSortedVertex(i)){
-	return false;
-      }
-    }
-    return true;
+  bool operator==(const Cell& c2) const {
+    //if(this->isCombined() != c2.isCombined()) return false;
+    return (this->getNum() == c2.getNum());
   }
+
 };
 
+
 // A cell that is a combination of cells of same dimension
 class CombinedCell : public Cell{
   
@@ -201,10 +212,6 @@ class CombinedCell : public Cell{
   // list of cells this cell is a combination of
   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);
@@ -212,20 +219,15 @@ class CombinedCell : public Cell{
   ~CombinedCell() {}
 
   int getDim() const { return _cells.begin()->first->getDim(); }
-  int getNumSortedVertices() const { return 0; }
-  int getSortedVertex(int vertex) const { return 0; }
   void getCells(std::map< Cell*, int, Less_Cell >& cells) { cells = _cells; }
   int getNumCells() const {return _cells.size();}  
-  int getNum() const { return _num; }
   bool hasVertex(int vertex) const;
 
   bool operator==(const Cell& c2) const {
-    if(this->isCombined() != c2.isCombined()) return false;
-    if(this->getNum() != c2.getNum()) return false;
-    else return true;
+    //if(this->isCombined() != c2.isCombined()) return false;
+    return (this->getNum() == c2.getNum());
   }
 };
-  
-#endif
+
 
 #endif
diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 9d34cf5539f095a7eea3836bd24d09254ccc6242..ffdf69bac958b0a4af29d3ebf0bc78987731fcad 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -1,88 +1,62 @@
-// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle
-//
-// See the LICENSE.txt file for license information. Please report all
-// bugs and problems to <gmsh@geuz.org>.
-//
-// Contributed by Matti Pellikka <matti.pellikka@tut.fi>.
-
 #include "CellComplex.h"
+#include "MElement.h"
 
-#if defined(HAVE_KBIPACK)
 
-CellComplex::CellComplex( std::vector<MElement*>& domainElements, 
-			  std::vector<MElement*>& subdomainElements)
+CellComplex::CellComplex(GModel* model,
+			 std::vector<MElement*>& domainElements, 
+			 std::vector<MElement*>& subdomainElements) :
+  _model(model), _dim(0), _simplicial(true), _saveorig(true)
 {
-  _dim = 0;
-  _simplicial = true;
 
-  // insert cells into cell complex
-  // subdomain needs to be inserted first!
-  if(!insert_cells(subdomainElements, true)){ panic_exit(); return; }
-  if(!insert_cells(domainElements, false)) { panic_exit(); return; }
+  _insertCells(subdomainElements, 1);
+  _insertCells(domainElements, 0);  
 
-  for(int i = 0; i < 4; i++){
-    _ocells[i] = _cells[i];
-    if(getSize(i) > _dim) _dim = i;
-  }  
-}
 
-void CellComplex::panic_exit(){
-  for(int i = 0; i < 4; i++){
-    for(citer cit = _ocells[i].begin(); cit != _ocells[i].end(); cit++){
+  int num = 0;
+  for(int dim = 0; dim < 4; dim++){
+    if(getSize(dim) != 0) _dim = dim;
+    _ocells[dim] = _cells[dim];
+    for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
       Cell* cell = *cit;
-      delete cell;
+      cell->setNum(++num);
+      cell->increaseGlobalNum();
+      cell->saveOriginalBd();
     }
   }
-  for(unsigned int i = 0; i < _newcells.size(); i++) delete _newcells.at(i);
-  printf("ERROR: No proper cell complex could be constructed!\n");
 }
 
-bool CellComplex::insert_cells(std::vector<MElement*>& elements,
-			       bool subdomain)
+bool CellComplex::_insertCells(std::vector<MElement*>& elements,
+			       int domain) 
 {
-  // add highest dimensional cells
   for(unsigned int i=0; i < elements.size(); i++){
-    MElement* element = elements.at(i);    
-    
-    int type = element->getType();   
-    if(_simplicial && !(type == TYPE_PNT || type == TYPE_LIN 
-		     || type == TYPE_TRI || type == TYPE_TET) ){
-      _simplicial = false;
-    }
-    //FIXME: no getFaceInfo methods for these MElements
-    if(type == TYPE_PYR){
-      printf("ERROR: mesh element %d not implemented! \n", type);
+    MElement* element = elements.at(i);
+    int type = element->getType();
+    if(type != TYPE_PNT && type != TYPE_LIN && 
+       type != TYPE_TRI && type != TYPE_TET) {
+      printf("Mesh element type %d not implemented in homology solver. \n", 
+	     type);
       return false;
     }
-    //Cell* cell = _cellPool.construct(element);
-    Cell* cell = new Cell(element);
-    cell->setInSubdomain(subdomain);
-    std::pair<citer, bool> insertInfo = _cells[cell->getDim()].insert(cell);
-    if(!insertInfo.second) delete cell; 
-    //if(!insertInfo.second) _cellPool.free(cell);  
+    Cell* cell = new Cell(element, domain);
+    bool insert = _cells[cell->getDim()].insert(cell).second;
+    if(!insert) delete cell;
   }
   
-  // add lower dimensional cells recursively
   for (int dim = 3; dim > 0; dim--){
     for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
       Cell* cell = *cit;
-      std::vector<MElement*> bdElements;
-      if(!cell->findBoundaryElements(bdElements)) return false;
-      for(unsigned int i = 0; i < bdElements.size(); i++){
-	//Cell* newCell = _cellPool.construct(bdElements.at(i));
-	Cell* newCell = new Cell(bdElements.at(i));
-	newCell->setInSubdomain(subdomain);
-	newCell->setDeleteImage(true);
-	std::pair<citer, bool> insertInfo = 
+      int numBdElements = cell->getNumBdElements();
+      for(int i = 0; i < numBdElements; i++){
+	Cell* newCell = new Cell(cell, i);
+	std::pair<citer, bool> insert = 
 	  _cells[newCell->getDim()].insert(newCell);
-	if(!insertInfo.second){ // the cell was already in the cell complex
-	  delete newCell; 
-	  //_cellPool.free(newCell); 
-	  newCell = *(insertInfo.first); 
+	if(!insert.second) { 
+	  delete newCell;
+	  newCell = *(insert.first);
 	}
-	if(!subdomain) {
-	  int ori = cell->findBoundaryCellOrientation(newCell);
-	  cell->addBoundaryCell( ori, newCell, true, true);
+	if(domain == 0) {
+	  int ori = cell->findBdCellOrientation(newCell);
+	  cell->addBoundaryCell( ori, newCell, true);
 	}
       }
     }
@@ -99,6 +73,7 @@ CellComplex::~CellComplex()
     }
   }
   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);
 }
 
 void CellComplex::insertCell(Cell* cell)
@@ -116,35 +91,40 @@ void CellComplex::insertCell(Cell* cell)
 void CellComplex::removeCell(Cell* cell, bool other)
 {  
   if(!hasCell(cell)) return;
-  std::map<Cell*, int, Less_Cell > coboundary;
+  std::map<Cell*, short int, Less_Cell > coboundary;
   cell->getCoboundary(coboundary);
-  std::map<Cell*, int, Less_Cell > boundary; 
+  std::map<Cell*, short int, Less_Cell > boundary; 
   cell->getBoundary(boundary);
-  
-  for(Cell::biter 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(Cell::biter 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);
+    }*/
+  
   _cells[cell->getDim()].erase(cell); 
+  if(!_saveorig && !cell->isCombined()) _removedcells.push_back(cell);
 }
 
-void CellComplex::removeCellQset(Cell* cell, 
-				 std::set<Cell*, Less_Cell>& Qset)
-{
-  Qset.erase(cell);
-}
-
-void CellComplex::enqueueCells(std::map<Cell*, int, Less_Cell>& cells, 
+void CellComplex::enqueueCells(std::map<Cell*, short int, Less_Cell>& cells, 
 			       std::queue<Cell*>& Q,
 			       std::set<Cell*, Less_Cell>& Qset)
 {
-  for(std::map<Cell*, int, Less_Cell>::iterator cit = cells.begin();
+  for(std::map<Cell*, short int, Less_Cell>::iterator cit = cells.begin();
       cit != cells.end(); cit++){
     Cell* cell = (*cit).first;
     citer it = Qset.find(cell);
@@ -166,14 +146,14 @@ int CellComplex::coreduction(Cell* startCell, bool omit,
   Q.push(startCell);
   Qset.insert(startCell);
   
-  std::map<Cell*, int, Less_Cell > bd_s;
-  std::map<Cell*, int, Less_Cell > cbd_c;
+  std::map<Cell*, short int, Less_Cell > bd_s;
+  std::map<Cell*, short int, Less_Cell > cbd_c;
 
   Cell* s;
   while( !Q.empty() ){
     s = Q.front();
     Q.pop();
-    removeCellQset(s, Qset);
+    Qset.erase(s);
     if(s->getBoundarySize() == 1 
        && inSameDomain(s, s->firstBoundary()->first) ){
       s->getBoundary(bd_s);
@@ -259,10 +239,10 @@ int CellComplex::coreduction(int dim, bool omit,
   }
   return count;
 }
-  
-int CellComplex::reduceComplex(bool omit)
+
+int CellComplex::reduceComplex(bool docombine, bool omit)
 {  
-  printf("Cell Complex: \n %d volumes, %d faces, %d edges and %d vertices. \n",
+  printf("Cell Complex reduction: \n %d volumes, %d faces, %d edges and %d vertices. \n",
 	 getSize(3), getSize(2), getSize(1), getSize(0));
 
   int count = 0;
@@ -293,13 +273,13 @@ int CellComplex::reduceComplex(bool omit)
   }
   
   printf(" %d volumes, %d faces, %d edges and %d vertices. \n",
-       getSize(3), getSize(2), getSize(1), getSize(0));
+	 getSize(3), getSize(2), getSize(1), getSize(0));
 
-  combine(3);
+  if(docombine) combine(3);
   reduction(2, false, empty);
-  combine(2);
+  if(docombine) combine(2);
   reduction(1, false, empty);
-  combine(1);
+  if(docombine) combine(1);
   
   printf(" %d volumes, %d faces, %d edges and %d vertices. \n",
 	 getSize(3), getSize(2), getSize(1), getSize(0));
@@ -319,9 +299,9 @@ void CellComplex::removeSubdomain()
   for(unsigned int i = 0; i < toRemove.size(); i++) removeCell(toRemove[i]);
 }
 
-int CellComplex::coreduceComplex(bool omit)
+int CellComplex::coreduceComplex(bool docombine, bool omit)
 {
-  printf("Cell Complex: \n %d volumes, %d faces, %d edges and %d vertices. \n",
+  printf("Cell Complex coreduction: \n %d volumes, %d faces, %d edges and %d vertices. \n",
 	 getSize(3), getSize(2), getSize(1), getSize(0));
   
   int count = 0;
@@ -359,11 +339,11 @@ int CellComplex::coreduceComplex(bool omit)
   printf(" %d volumes, %d faces, %d edges and %d vertices. \n",
 	 getSize(3), getSize(2), getSize(1), getSize(0));
   
-  cocombine(0);
+  if(docombine) cocombine(0);
   coreduction(1, false, empty);
-  cocombine(1);
+  if(docombine) cocombine(1);
   coreduction(2, false, empty);
-  cocombine(2);
+  if(docombine) cocombine(2);
   coreduction(3, false, empty);
   coherent();
 
@@ -373,116 +353,112 @@ int CellComplex::coreduceComplex(bool omit)
   return 0;
 }
 
-int CellComplex::cocombine(int dim)
-{ 
-  //printf("Cell complex before cocombining: %d volumes, %d faces, %d edges and %d vertices.\n", getSize(3), getSize(2), getSize(1), getSize(0));
-    
-  if(dim < 0 || dim > 2) return 0;
+int CellComplex::combine(int dim)
+{
+  //printf("Cell complex before combining: %d volumes, %d faces, %d edges and %d vertices.\n",  getSize(3), getSize(2), getSize(1), getSize(0));
+  if(dim < 1 || dim > 3) return 0;
   
   std::queue<Cell*> Q;
   std::set<Cell*, Less_Cell> Qset;
-  std::map<Cell*, int, Less_Cell> cbd_c;
+  std::map<Cell*, short int, Less_Cell> bd_c;
   int count = 0;
-  
+
   for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
     Cell* cell = *cit;
-    cell->getCoboundary(cbd_c);
-    enqueueCells(cbd_c, Q, Qset);
-
+    cell->getBoundary(bd_c);
+    enqueueCells(bd_c, Q, Qset);
     while(Q.size() != 0){
       Cell* s = Q.front();
-      Q.pop();
-      if(s->getBoundarySize() == 2){
-	Cell::biter it = s->firstBoundary();
-        int or1 = (*it).second;
-        Cell* c1 = (*it).first;
-        it++;
-        int or2 = (*it).second;
-        Cell* c2 = (*it).first;
+      Q.pop(); 
+
+      if(s->getCoboundarySize() == 2){
+	Cell::biter it = s->firstCoboundary();
+        int or1 = it->second.get();
+        Cell* c1 = it->first;
+	it++;
+        while (it->second.get() == 0) it++;
+        int or2 = it->second.get();
+        Cell* c2 = it->first;
 
         if(!(*c1 == *c2) && abs(or1) == abs(or2)
            && inSameDomain(s, c1) && inSameDomain(s, c2)){
-	  removeCell(s);
-          
-          c1->getCoboundary(cbd_c);
-          enqueueCells(cbd_c, Q, Qset);
-          c2->getCoboundary(cbd_c);
-          enqueueCells(cbd_c, Q, Qset);
-          
-          CombinedCell* newCell = new CombinedCell(c1, c2, 
-						   (or1 != or2), true );
+          removeCell(s);
+
+          c1->getBoundary(bd_c);
+          enqueueCells(bd_c, Q, Qset);
+          c2->getBoundary(bd_c);
+          enqueueCells(bd_c, Q, Qset);
+
+	  CombinedCell* newCell = new CombinedCell(c1, c2, (or1 != or2));
           removeCell(c1);
           removeCell(c2);
           insertCell(newCell);
-          
+	    
           cit = firstCell(dim);
           count++;
         }
       }
-      removeCellQset(s, Qset);
-      
+      Qset.erase(s);
     }
   }
-  //printf("Cell complex after cocombining: %d volumes, %d faces, %d edges and %d vertices.\n", getSize(3), getSize(2), getSize(1), getSize(0));
-  
+  //printf("Cell complex after combining: %d volumes, %d faces, %d edges and %d vertices.\n", getSize(3), getSize(2), getSize(1), getSize(0));
+    
   return count;
 }
 
-int CellComplex::combine(int dim)
-{
-  //printf("Cell complex before combining: %d volumes, %d faces, %d edges and %d vertices.\n",  getSize(3), getSize(2), getSize(1), getSize(0));
+
+int CellComplex::cocombine(int dim)
+{ 
+  //printf("Cell complex before cocombining: %d volumes, %d faces, %d edges and %d vertices.\n", getSize(3), getSize(2), getSize(1), getSize(0));
     
-  if(dim < 1 || dim > 3) return 0;
+  if(dim < 0 || dim > 2) return 0;
   
   std::queue<Cell*> Q;
   std::set<Cell*, Less_Cell> Qset;
-  std::map<Cell*, int, Less_Cell> bd_c;
+  std::map<Cell*, short int, Less_Cell> cbd_c;
   int count = 0;
-
   
   for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
     Cell* cell = *cit;
-    cell->getBoundary(bd_c);
-    enqueueCells(bd_c, Q, Qset);
+    cell->getCoboundary(cbd_c);
+    enqueueCells(cbd_c, Q, Qset);
 
     while(Q.size() != 0){
-      
       Cell* s = Q.front();
-      Q.pop(); 
+      Q.pop();
+      if(s->getBoundarySize() == 2){
 
-      if(s->getCoboundarySize() == 2){
-	Cell::biter it = s->firstCoboundary();
-        int or1 = (*it).second;
-        Cell* c1 = (*it).first;
-        it++;
-        int or2 = (*it).second;
-        Cell* c2 = (*it).first;
+	Cell::biter it = s->firstBoundary();
+        int or1 = it->second.get();
+        Cell* c1 = it->first;
+	it++;
+        while (it->second.get() == 0) it++;
+        int or2 = it->second.get();
+        Cell* c2 = it->first;
 
         if(!(*c1 == *c2) && abs(or1) == abs(or2)
            && inSameDomain(s, c1) && inSameDomain(s, c2)){
-          removeCell(s);
+	  removeCell(s);
           
-          c1->getBoundary(bd_c);
-          enqueueCells(bd_c, Q, Qset);
-          c2->getBoundary(bd_c);
-          enqueueCells(bd_c, Q, Qset);
-
-	  CombinedCell* newCell = new CombinedCell(c1, c2, (or1 != or2));
-          //CombinedCell* newCell = _ccellPool.construct(c1, c2, (or1 != or2));
+          c1->getCoboundary(cbd_c);
+          enqueueCells(cbd_c, Q, Qset);
+          c2->getCoboundary(cbd_c);
+          enqueueCells(cbd_c, Q, Qset);
+          
+          CombinedCell* newCell = new CombinedCell(c1, c2, 
+						   (or1 != or2), true );
           removeCell(c1);
           removeCell(c2);
           insertCell(newCell);
-	  
+          
           cit = firstCell(dim);
           count++;
         }
       }
-      removeCellQset(s, Qset);
-      
+      Qset.erase(s);
     }
   }
-  //printf("Cell complex after combining: %d volumes, %d faces, %d edges and %d vertices.\n", getSize(3), getSize(2), getSize(1), getSize(0));
-  
+  //printf("Cell complex after cocombining: %d volumes, %d faces, %d edges and %d vertices.\n", getSize(3), getSize(2), getSize(1), getSize(0));
   
   return count;
 }
@@ -493,9 +469,9 @@ bool CellComplex::coherent()
   for(int i = 0; i < 4; i++){
     for(citer cit = firstCell(i); cit != lastCell(i); cit++){
       Cell* cell = *cit;
-      std::map<Cell*, int, Less_Cell> boundary;
+      std::map<Cell*, short int, Less_Cell> boundary;
       cell->getBoundary(boundary);
-      for(Cell::biter it = boundary.begin();
+      for(std::map<Cell*, short int, Less_Cell>::iterator it = boundary.begin();
 	  it != boundary.end(); it++){
         Cell* bdCell = (*it).first;
         int ori = (*it).second;
@@ -507,14 +483,14 @@ bool CellComplex::coherent()
         }
         if(!bdCell->hasCoboundary(cell)){
           printf("Warning! Incoherent boundary/coboundary pair! Fixed. \n");
-	  bdCell->addCoboundaryCell(ori, cell, false, false);
+	  bdCell->addCoboundaryCell(ori, cell, false);
           coherent = false;
         }
         
       }
-      std::map<Cell*, int, Less_Cell> coboundary;
+      std::map<Cell*, short int, Less_Cell> coboundary;
       cell->getCoboundary(coboundary);
-      for(Cell::biter it = coboundary.begin();
+      for(std::map<Cell*, short int, Less_Cell>::iterator it = coboundary.begin();
 	  it != coboundary.end(); it++){
         Cell* cbdCell = (*it).first;
         int ori = (*it).second;
@@ -526,7 +502,7 @@ bool CellComplex::coherent()
         }
         if(!cbdCell->hasBoundary(cell)){
           printf("Warning! Incoherent coboundary/boundary pair! Fixed. \n");
-	  cbdCell->addBoundaryCell(ori, cell, false, false);
+	  cbdCell->addBoundaryCell(ori, cell, false);
           coherent = false;
         }
         
@@ -558,31 +534,147 @@ void  CellComplex::getCells(std::set<Cell*, Less_Cell>& cells,
   }
 }
 
-void CellComplex::restoreComplex()
+bool CellComplex::restoreComplex()
 {
-  for(int i = 0; i < 4; i++){
-    _cells[i] = _ocells[i];
-    for(citer cit = firstCell(i); cit != lastCell(i); cit++){
-      Cell* cell = *cit;
-      cell->restoreCell();
+  if(_saveorig){
+    for(int i = 0; i < 4; i++){
+      _cells[i] = _ocells[i];
+      for(citer cit = firstCell(i); cit != lastCell(i); cit++){
+	Cell* cell = *cit;
+	cell->restoreCell();
+      }
     }
+    for(unsigned int i = 0; i < _newcells.size(); i++){
+      Cell* cell = _newcells.at(i);
+      delete cell;
+    }
+    _newcells.clear();
+    return true;
   }
-  for(unsigned int i = 0; i < _newcells.size(); i++){
-    Cell* cell = _newcells.at(i);
-    delete cell;
-  }
-  _newcells.clear();
+  else return false;
 }
 
 void CellComplex::printComplex(int dim)
 {
+  if(getSize(dim) == 0) printf("Cell complex dimension %d is empty. \n", dim);
   for (citer cit = firstCell(dim); cit != lastCell(dim); cit++){
     Cell* cell = *cit;
     cell->printCell();
-    cell->printBoundary();
-    cell->printCoboundary();
+    //cell->printBoundary();
+    //cell->printCoboundary();
     //printf("--- \n");
   }
 }
 
-#endif
+int CellComplex::saveComplex(std::string filename)
+{
+  /*FILE *fp = fopen (filename.c_str(), "w");
+  if(!fp){
+    printf("\nUnable to open file '%s' \n", filename.c_str());
+    return 0;
+  }
+  printf("\nWriting file '%s' ... \n", filename.c_str());
+  
+  fprintf(fp, "$Cells\n");
+  fprintf(fp, "%d\n", getSize(0)+getSize(1)+getSize(2)+getSize(3));
+  for(int dim = 0; dim < 4; dim++){
+    for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
+      Cell* cell = *cit;
+      fprintf(fp, "%d %d %d %d %d", cell->getNum(), cell->getType(), 
+	      1, cell->getDomain(), cell->getNumVertices());
+      for(int i = 0; i < cell->getNumVertices(); i++){
+	fprintf(fp, " %d", cell->getVertex(i));
+      }
+      fprintf(fp, " %d", cell->getBoundarySize());
+      for(Cell::biter bit = cell->firstBoundary();
+	  bit != cell->lastBoundary(); bit++){
+	fprintf(fp, " %d %d", bit->first->getNum(), bit->second);
+      }
+      fprintf(fp, " %d", cell->getCoboundarySize());
+      for(Cell::biter bit = cell->firstCoboundary();
+	  bit != cell->lastCoboundary(); bit++){
+	fprintf(fp, " %d %d", bit->first->getNum(), bit->second);
+      }
+      fprintf(fp, "\n");
+    }
+  }
+
+  fclose(fp);
+
+  printf("Wrote %d cells to '%s' \n", 
+	 getSize(0)+getSize(1)+getSize(2)+getSize(3), filename.c_str());
+
+	 return 1;*/
+}
+
+int CellComplex::loadComplex(std::string filename)
+{
+  /*  FILE *fp = fopen (filename.c_str(), "r");
+  if(!fp){
+    printf("\nUnable to open file '%s' \n", filename.c_str());
+    return 0;
+  }
+
+  std::map<int, Cell*> numToCell;
+
+  char str[256] = "XXX";
+  while(1) {
+
+    while(str[0] != '$'){
+      if(!fgets(str, sizeof(str), fp) || feof(fp))
+        break;
+    }
+    
+    if(feof(fp))
+      break;
+    
+    if(!strncmp(&str[1], "Cells", 5)) {
+      if(!fgets(str, sizeof(str), fp)) return 0;
+      int numCells;
+      sscanf(str, "%d", &numCells);
+      for(int i = 0; i < numCells; i++) {
+	int num, type, numTags;
+	std::vector<int> domain;
+	int tag;
+	if(fscanf(fp, "%d %d %d", &num, &type, &numTags) != 3) return 0;
+	for(int j = 0; j < numTags; j++){
+	  if(fscanf(fp, "%d", &tag) != 1) return 0;
+	  domain.push_back(tag);
+	}
+	
+	std::vector<int> vertices;
+	if(fscanf(fp, "%d", &numTags) != 1) return 0;
+	for(int j = 0; j < numTags; j++){
+	  if(fscanf(fp, "%d", &tag) != 1) return 0;
+	  vertices.push_back(tag);
+	}
+
+	int dim = 0;
+	if(type == 1){
+	  if(vertices.size() == 2) dim = 1;
+	  else if(vertices.size() == 3) dim = 2;
+	  else if(vertices.size() == 4) dim = 3;
+	}
+
+	Cell* cell = new Cell(num, dim, type, domain, vertices);
+	numToCell[num] = cell;
+	
+
+	int numCell;
+	if(fscanf(fp, "%d", &numTags) != 1) return 0;
+	for(int j = 0; j < numTags; j++){
+	  if(fscanf(fp, "%d %d", &numCell, &tag) != 1) return 0;
+	}
+	if(fscanf(fp, "%d", &numTags) != 1) return 0;
+	for(int j = 0; j < numTags; j++){
+	  if(fscanf(fp, "%d %d", &numCell, &tag) != 1) return 0;
+	}
+
+      }
+    }
+    
+  }
+
+  fclose(fp);
+  return 1;*/
+}
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index 07c9180612bc9c1557119429fc4f3f32a9662b72..6a601372c497f747c39f87469eacabb6853c670c 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -1,59 +1,57 @@
-// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle
+// Homology Solver 
+// 
+//  Copyright (C) 2010 Matti Pellikka and Saku Suuriniemi
+//  Tampere University of Technology, Electromagnetics
 //
-// See the LICENSE.txt file for license information. Please report all
-// bugs and problems to <gmsh@geuz.org>.
-//
-// Contributed by Matti Pellikka <matti.pellikka@tut.fi>.
+//  See LICENCE.txt for license information.
+//  See README.txt for more information.
 
 #ifndef _CELLCOMPLEX_H_
 #define _CELLCOMPLEX_H_
 
-#include "GmshConfig.h"
-
-#if defined(HAVE_KBIPACK)
-
-#include <stdio.h>
-#include <string>
-#include <algorithm>
+#include <map>
+#include <string.h>
 #include <set>
+#include <algorithm>
 #include <queue>
+#include <string>
 #include "Cell.h"
 #include "MElement.h"
-#include "ChainComplex.h"
-//#include "GModel.h"
-//#include <boost/pool/object_pool.hpp>
+#include "GModel.h"
 
 class Cell;
+class BdInfo;
+
 
-// A class representing a cell complex made out of a finite element mesh.
 class CellComplex
 {
+
  private:
+
+  GModel* _model;
+  
   // sorted containers of unique cells in this cell complex 
   // one for each dimension
   std::set<Cell*, Less_Cell> _cells[4];
-    
+
   // original cells of this cell complex
-  std::set<Cell*, Less_Cell>  _ocells[4];
-  
+  std::set<Cell*, Less_Cell> _ocells[4];
+
   // new cells created during reductions
   std::vector<Cell*> _newcells;
-  
-  int _dim;
-  
-  // is the cell complex simplicial
+  std::vector<Cell*> _removedcells;
+
   bool _simplicial;
+  int _dim;
+  bool _saveorig;
+
+  // for constructor 
+  bool _insertCells(std::vector<MElement*>& elements,  int domain);
   
   // enqueue cells in queue if they are not there already
-  void enqueueCells(std::map<Cell*, int, Less_Cell>& cells, 
+  void enqueueCells(std::map<Cell*, short 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);
   
-  // for constructor 
-  bool insert_cells(std::vector<MElement*>& elements, bool subdomain);
-  void panic_exit();
-
   // insert/remove a cell from this cell complex
   void removeCell(Cell* cell, bool other=true);
   void insertCell(Cell* cell);
@@ -61,27 +59,29 @@ class CellComplex
   // queued coreduction
   int coreduction(Cell* startCell, bool omit, 
 		  std::vector<Cell*>& omittedCells);
-  
- public: 
-  
-  CellComplex( std::vector<MElement*>& domainElements, 
-	       std::vector<MElement*>& subdomainElements);
+
+ public:
+  CellComplex(GModel* model,
+	      std::vector<MElement*>& domainElements, 
+	      std::vector<MElement*>& subdomainElements);
   ~CellComplex();
+  
+
+  GModel* getModel() const { return _model; }
+  int getDim() { return _dim; } 
+  bool simplicial() { return _simplicial; }
 
   // get the number of certain dimensional cells
   int getSize(int dim, bool orig=false){ 
     if(!orig) return _cells[dim].size();
     else return _ocells[dim].size(); }
   
-  int getDim() {return _dim; } 
-  bool simplicial() { return _simplicial; }
-  
   // get dim-dimensional cells
   // domain = 0: cells in domain relative to subdomain
   // domain = 1: cells in domain
   // domain = 2: cells in subdomain
   void getCells(std::set<Cell*, Less_Cell>& cells, int dim, int domain=0);
-  std::set<Cell*, Less_Cell> getOrigCells(int dim){ return _ocells[dim]; }
+  //std::set<Cell*, Less_Cell> getOrigCells(int dim){ return _ocells[dim]; }
   
   // iterator for the cells of same dimension
   typedef std::set<Cell*, Less_Cell>::iterator citer;
@@ -97,8 +97,7 @@ class CellComplex
   
   // check whether two cells both belong to subdomain or if neither one does
   bool inSameDomain(Cell* c1, Cell* c2) const { 
-    return ( (!c1->inSubdomain() && !c2->inSubdomain()) 
-	     || (c1->inSubdomain() && c2->inSubdomain() )); }
+    return (c1->getDomain() ==  c2->getDomain()); }
   
   // remove cells in subdomain from this cell complex
   void removeSubdomain();
@@ -120,8 +119,8 @@ class CellComplex
 
   // full (co)reduction of this cell complex (all dimensions, with combining)
   // (with highest dimensional cell omitting?)
-  int reduceComplex(bool omit=true);
-  int coreduceComplex(bool imit=true);
+  int reduceComplex(bool docombine=true, bool omit=true);
+  int coreduceComplex(bool docombine=true, bool omit=true);
    
   int eulerCharacteristic(){ 
     return getSize(0) - getSize(1) + getSize(2) - getSize(3);}
@@ -129,12 +128,16 @@ class CellComplex
     printf("Euler characteristic: %d. \n", eulerCharacteristic()); }
   
   // restore the cell complex to its original state before (co)reduction
-  void restoreComplex();
+  bool restoreComplex();
 
   // print the vertices of cells of certain dimension
   void printComplex(int dim);
-};
 
-#endif
+
+  // experimental
+  int saveComplex(std::string filename);
+  int loadComplex(std::string filename);
+
+};
 
 #endif
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index 10f87d86d6bcae6027e2ff63649f8d03ded50b04..4c73f7ee054ba7c5704d83902e23bdd257ac4094 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -1,15 +1,13 @@
-
-// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle
-//
-// See the LICENSE.txt file for license information. Please report all
-// bugs and problems to <gmsh@geuz.org>.
+// Homology Solver 
+// 
+//  Copyright (C) 2010 Matti Pellikka and Saku Suuriniemi
+//  Tampere University of Technology, Electromagnetics
 //
-// Contributed by Matti Pellikka <matti.pellikka@tut.fi>.
+//  See LICENCE.txt for license information.
+//  See README.txt for more information.
 
 #include "ChainComplex.h"
 
-#if defined(HAVE_KBIPACK)
-
 ChainComplex::ChainComplex(CellComplex* cellComplex, int domain)
 { 
   _dim = cellComplex->getDim();
@@ -68,7 +66,8 @@ ChainComplex::ChainComplex(CellComplex* cellComplex, int domain)
 	    || (domain == 2 && cell->inSubdomain()) ){
           for(Cell::biter it = cell->firstBoundary();
 	      it != cell->lastBoundary(); it++){
-            Cell* bdCell = (*it).first;
+            Cell* bdCell = it->first;
+	    if(it->second.get() == 0) continue;
             if((domain == 0 && !bdCell->inSubdomain()) || domain == 1 
 	       || (domain == 2 && cell->inSubdomain()) ){
               int old_elem = 0;
@@ -83,8 +82,8 @@ ChainComplex::ChainComplex(CellComplex* cellComplex, int domain)
                 gmp_matrix_get_elem(elem, bdCell->getIndex(), 
 				    cell->getIndex(), _HMatrix[dim]);
                 old_elem = mpz_get_si(elem);
-                mpz_set_si(elem, old_elem + (*it).second);
-                if( abs((old_elem + (*it).second)) > 1){
+                mpz_set_si(elem, old_elem + it->second.get());
+                if( abs((old_elem + it->second.get())) > 1){
 		  //printf("Incidence index: %d, in HMatrix: %d. \n", (old_elem + (*it).second), dim);
                 }
                 gmp_matrix_set_elem(elem, bdCell->getIndex(), 
@@ -307,7 +306,6 @@ void ChainComplex::computeHomology(bool dual)
       //KerCod(highDim);
     }
     
-    //printf("Homology computation process: step %d of 4 \n", i+1);
     
     KerCod(highDim);
     
@@ -369,7 +367,7 @@ void ChainComplex::computeHomology(bool dual)
         
         gmp_matrix_right_mult(getHbasis(setDim), getQMatrix(lowDim));
       } 
-    } 
+    }
 
     //destroy_gmp_matrix(getKerHMatrix(lowDim));
     //destroy_gmp_matrix(getCodHMatrix(lowDim));
@@ -380,7 +378,7 @@ void ChainComplex::computeHomology(bool dual)
     //setCodHMatrix(lowDim, NULL);
     setJMatrix(lowDim, NULL);
     setQMatrix(lowDim, NULL);  
-  } 
+  }
   return;
 }
 
@@ -566,17 +564,18 @@ bool ChainComplex::deformChain(std::map<Cell*, int, Less_Cell>& cells,
 {  
   Cell* c1 = cell.first;
   int dim = c1->getDim();
-  for(citer cit = c1->firstCoboundary(true); cit != c1->lastCoboundary(true);
+  for(Cell::biter cit = c1->firstCoboundary(true); cit != c1->lastCoboundary();
       cit++){
     
     std::map<Cell*, int, Less_Cell> cellsInChain;
     std::map<Cell*, int, Less_Cell> cellsNotInChain;
-    Cell* c1CbdCell = (*cit).first;
+    Cell* c1CbdCell = cit->first;
 
-    for(citer cit2 = c1CbdCell->firstBoundary(true);
-	cit2 != c1CbdCell->lastBoundary(true); cit2++){
-      Cell* c1CbdBdCell = (*cit2).first;
-      int coeff = (*cit2).second;
+    for(Cell::biter cit2 = c1CbdCell->firstBoundary(true);
+	cit2 != c1CbdCell->lastBoundary(); cit2++){
+      Cell* c1CbdBdCell = cit2->first;
+      int coeff = cit2->second.geto();
+      if(coeff == 0) continue;
       if( (cells.find(c1CbdBdCell) != cells.end() && cells[c1CbdBdCell] != 0) 
 	  || c1CbdBdCell->inSubdomain()){
 	cellsInChain.insert(std::make_pair(c1CbdBdCell, coeff));
@@ -633,7 +632,7 @@ void ChainComplex::smoothenChain(std::map<Cell*, int, Less_Cell>& cells)
 {
   if(!_cellComplex->simplicial() || cells.empty()) return;
   int dim = cells.begin()->first->getDim();
-  int start = cells.size();
+  int start = cells.size(); 
 
   int useless = 0;
   for(int i = 0; i < 20; i++){
@@ -1031,5 +1030,3 @@ void HomologySequence::blockHBasis(gmp_matrix* block1T,
   destroy_gmp_matrix(temp2);
 }
 
-
-#endif
diff --git a/Geo/ChainComplex.h b/Geo/ChainComplex.h
index 40489a540e92256f28a59dd6bf907ac8a249a3e6..eb2c4c367cde873525a7f2616dacf58a2e8e70ff 100644
--- a/Geo/ChainComplex.h
+++ b/Geo/ChainComplex.h
@@ -1,15 +1,7 @@
-// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle
-//
-// See the LICENSE.txt file for license information. Please report all
-// bugs and problems to <gmsh@geuz.org>.
-//
-// Contributed by Matti Pellikka <matti.pellikka@tut.fi>.
-
 #ifndef _CHAINCOMPLEX_H_
 #define _CHAINCOMPLEX_H_
 
 #include "GmshConfig.h"
-
 #if defined(HAVE_KBIPACK)
 
 #include <cstdio>
@@ -172,6 +164,7 @@ class ChainComplex
       gmp_matrix_right_mult(_Hbasis[dim], T);
     }
   }
+  //void printBasisChain(std::map<Cell*, int, Less_Cell>& chain);
 
   // debugging aid
   int printMatrix(gmp_matrix* matrix){ 
@@ -182,7 +175,7 @@ class ChainComplex
 };
 
 
-// An experimental class to modify computed bases for homnology spaces
+// An experimental class to modify computed bases for homology spaces
 // so that the basis chains are decomposed according to the long
 // exact homology sequence.
 class HomologySequence
@@ -237,6 +230,5 @@ class HomologySequence
 };
 
 #endif
-   
-#endif
+#endif   
    
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index 6914e37f0ad426f00c85280d167ad252fd40b0d2..45df3123c05b5f65bdf06033e5acd6cc69cfcb93 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -14,12 +14,12 @@
 #if defined(HAVE_KBIPACK)
 
 Homology::Homology(GModel* model, std::vector<int> physicalDomain, 
-		   std::vector<int> physicalSubdomain)
+		   std::vector<int> physicalSubdomain, 
+		   bool combine, bool omit, bool smoothen) :
+  _model(model), _domain(physicalDomain), _subdomain(physicalSubdomain), 
+  _combine(combine), _omit(omit), _smoothen(smoothen)
 { 
-  _model = model; 
-  _domain = physicalDomain;
-  _subdomain = physicalSubdomain;
-  _fileName = "";
+  _fileName = "";  
 
   // default to the whole model
   if(_domain.empty()){
@@ -58,7 +58,8 @@ Homology::Homology(GModel* model, std::vector<int> physicalDomain,
       }
     }
   }
-  
+  _maxdomain = 0;
+  _maxnum = 0;
 }
 
 CellComplex* Homology::createCellComplex(std::vector<GEntity*>& domainEntities,
@@ -86,7 +87,8 @@ CellComplex* Homology::createCellComplex(std::vector<GEntity*>& domainEntities,
     }
   }
 
-  CellComplex* cellComplex =  new CellComplex(domainElements, 
+  CellComplex* cellComplex =  new CellComplex(_model, 
+					      domainElements, 
 					      subdomainElements);
 
   if(cellComplex->getSize(0) == 0){ 
@@ -119,7 +121,7 @@ void Homology::findGenerators(CellComplex* cellComplex)
   Msg::StatusBar(2, true, "Reducing cell complex...");
 
   double t1 = Cpu();
-  int omitted = cellComplex->reduceComplex();
+  int omitted = cellComplex->reduceComplex(_combine, _omit);
   
   double t2 = Cpu();
   Msg::StatusBar(2, true, "Done reducing cell complex (%g s)", t2 - t1);
@@ -129,8 +131,8 @@ void Homology::findGenerators(CellComplex* cellComplex)
   
   Msg::StatusBar(2, true, "Computing homology spaces...");
   t1 = Cpu();
-  ChainComplex* chains = new ChainComplex(cellComplex);
-  chains->computeHomology();
+  ChainComplex chains = ChainComplex(cellComplex);
+  chains.computeHomology();
   t2 = Cpu();
   Msg::StatusBar(2, true, "Done computing homology spaces (%g s)", t2 - t1);
   
@@ -139,18 +141,19 @@ void Homology::findGenerators(CellComplex* cellComplex)
     HRank[j] = 0;
     std::string dimension = "";
     convert(j, dimension);
-    for(int i = 1; i <= chains->getBasisSize(j, 3); i++){
+    for(int i = 1; i <= chains.getBasisSize(j, 3); i++){
       
       std::string generator = "";
       convert(i, generator);
       
       std::string name = "H" + dimension + domainString + generator;
       std::map<Cell*, int, Less_Cell> protoChain;
-      chains->getBasisChain(protoChain, i, j, 3, true);
-      Chain* chain = new Chain(protoChain, cellComplex, _model, name, 
-			       chains->getTorsion(j,i));      
+      chains.getBasisChain(protoChain, i, j, 3, _smoothen);
+      Chain* chain = new Chain(protoChain, cellComplex, ++_maxdomain, name, 
+			       chains.getTorsion(j,i));      
       if(chain->getSize() == 0) {
 	delete chain;
+	_maxdomain--;
 	continue;
       }
       HRank[j] = HRank[j] + 1;
@@ -158,15 +161,13 @@ void Homology::findGenerators(CellComplex* cellComplex)
 	Msg::Warning("H%d %d has torsion coefficient %d!", 
 		     j, i, chain->getTorsion());
       }
-      _generators.push_back(chain->createPGroup());
-      delete chain;
+      _basisChains[chain->createPGroup()] = chain;
     }
   }
   
   if(_fileName != "") writeGeneratorsMSH();
-  
+
   if(ownComplex) delete cellComplex;
-  delete chains;
   
   Msg::Info("Ranks of homology spaces for primal cell complex:");
   Msg::Info("H0 = %d", HRank[0]);
@@ -191,7 +192,7 @@ void Homology::findDualGenerators(CellComplex* cellComplex)
   Msg::StatusBar(2, true, "Reducing cell complex...");
   
   double t1 = Cpu();
-  int omitted = cellComplex->coreduceComplex();
+  int omitted = cellComplex->coreduceComplex(_combine, _omit);
   double t2 = Cpu();
   
   Msg::StatusBar(2, true, "Done reducing cell complex (%g s)", t2 - t1);
@@ -201,9 +202,9 @@ void Homology::findDualGenerators(CellComplex* cellComplex)
    
   Msg::StatusBar(2, true, "Computing homology spaces...");
   t1 = Cpu();
-  ChainComplex* chains = new ChainComplex(cellComplex);
-  chains->transposeHMatrices();
-  chains->computeHomology(true);
+  ChainComplex chains = ChainComplex(cellComplex);
+  chains.transposeHMatrices();
+  chains.computeHomology(true);
   t2 = Cpu();
   Msg::StatusBar(2, true, "Done computing homology spaces (%g s)", t2- t1);
   
@@ -215,7 +216,7 @@ void Homology::findDualGenerators(CellComplex* cellComplex)
     std::string dimension = "";
     convert(dim-j, dimension);
 
-    for(int i = 1; i <= chains->getBasisSize(j, 3); i++){
+    for(int i = 1; i <= chains.getBasisSize(j, 3); i++){
       
       std::string generator = "";
       convert(i, generator);
@@ -223,26 +224,27 @@ void Homology::findDualGenerators(CellComplex* cellComplex)
       std::string name = "H" + dimension + "*" + 
 	getDomainString(_domain, _subdomain) + generator;
       std::map<Cell*, int, Less_Cell> protoChain;
-      chains->getBasisChain(protoChain, i, j, 3);
-      Chain* chain = new Chain(protoChain, cellComplex, _model, name, 
-			       chains->getTorsion(j,i));
+      chains.getBasisChain(protoChain, i, j, 3, _smoothen);
+      Chain* chain = new Chain(protoChain, cellComplex, ++_maxdomain, name, 
+			       chains.getTorsion(j,i));
       if(chain->getSize() == 0) {
 	delete chain;
+	--_maxdomain;
 	continue;
       }
-      _generators.push_back(chain->createPGroup());
+
       HRank[dim-j] = HRank[dim-j] + 1;
       if(chain->getTorsion() != 1){ 
 	Msg::Warning("H%d* %d has torsion coefficient %d!", 
 		     dim-j, i, chain->getTorsion());
-      }     
+      }
+      _basisChains[chain->createPGroup()] = chain;
     }
   }
 
   if(_fileName != "") writeGeneratorsMSH();
 
   if(ownComplex) delete cellComplex;
-  delete chains;
   
   Msg::Info("Ranks of homology spaces for the dual cell complex:");
   Msg::Info("H0* = %d", HRank[0]);
@@ -255,7 +257,7 @@ void Homology::findDualGenerators(CellComplex* cellComplex)
 		 HRank[0], HRank[1], HRank[2], HRank[3]);
 }
 
-void Homology::findHomSequence(){
+/*void Homology::findHomSequence(){
   CellComplex* cellComplex = createCellComplex(_domainEntities, 
 					       _subdomainEntities);
   Msg::StatusBar(2, true, "Reducing cell complex...");
@@ -362,7 +364,7 @@ void Homology::findHomSequence(){
   delete cellComplex;
   delete seq;
 }
-
+*/
 std::string Homology::getDomainString(const std::vector<int>& domain,
 				      const std::vector<int>& subdomain) 
 {
@@ -403,39 +405,11 @@ bool Homology::writeGeneratorsMSH(bool binary)
   Msg::Info("Wrote homology computation results to %s", _fileName.c_str());
   return true;
 }
-Chain::Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs, 
-	     CellComplex* cellComplex, GModel* model,
-	     std::string name, int torsion)
-{  
-  int i = 0;
-  for(std::set<Cell*, Less_Cell>::iterator cit = cells.begin();
-      cit != cells.end(); cit++){
-    Cell* cell = *cit;
-    _dim = cell->getDim();
-    if((int)coeffs.size() > i){
-      if(coeffs.at(i) != 0){
-        std::map<Cell*, int, Less_Cell > subCells;
-	cell->getCells(subCells);
-        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));
-        }
-      }
-      i++;
-    }
-    
-  }
-  _name = name;
-  _cellComplex = cellComplex;
-  _torsion = torsion;
-  _model = model;
-  
-}
 
 void Homology::storeCells(CellComplex* cellComplex, int dim)
 {
   std::vector<MElement*> elements;
+  MElementFactory factory;
 
   for(CellComplex::citer cit = cellComplex->firstCell(dim); 
       cit != cellComplex->lastCell(dim); cit++){
@@ -446,8 +420,10 @@ void Homology::storeCells(CellComplex* cellComplex, int dim)
     for(Cell::citer it = cells.begin(); it != cells.end(); it++){
       Cell* subCell = it->first;
       
-      MElement* e = subCell->getImageMElement();
-      subCell->setDeleteImage(false);
+      std::vector<MVertex*> v;
+      cell->getMeshVertices(v);
+
+      MElement* e = factory.create(cell->getTypeMSH(), v);
       elements.push_back(e);
     }
   }
@@ -469,20 +445,6 @@ void Homology::storeCells(CellComplex* cellComplex, int dim)
   _model->setPhysicalName("Cell Complex", dim, physicalNum);
 }
 
-Chain::Chain(std::map<Cell*, int, Less_Cell>& chain, 
-	     CellComplex* cellComplex, GModel* model,
-	     std::string name, int torsion)
-{  
-  _cells = chain;
-  if(!_cells.empty()) _dim = firstCell()->first->getDim();
-  else _dim = 0;
-  _name = name;
-  _cellComplex = cellComplex;
-  _torsion = torsion;
-  _model = model; 
-  eraseNullCells();
-}
-
 int Chain::writeChainMSH(const std::string &name)
 {  
   if(getSize() == 0) return 1;
@@ -508,7 +470,7 @@ int Chain::writeChainMSH(const std::string &name)
   for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
     Cell* cell = (*cit).first;
     int coeff = (*cit).second;
-    fprintf(fp, "%d %d \n", cell->getImageMElement()->getNum(), coeff );
+    fprintf(fp, "%d %d \n", cell->getIndex(), coeff );
   }
   
   fprintf(fp, "$EndElementData\n");
@@ -526,33 +488,32 @@ int Chain::createPGroup()
   for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
     Cell* cell = (*cit).first;
     int coeff = (*cit).second;
-
-    MElement* e = cell->getImageMElement();
     std::vector<MVertex*> v;
-    for(int i = 0; i < e->getNumVertices(); i++){
-      v.push_back(e->getVertex(i));
-    }
-    MElement* ne = factory.create(e->getTypeForMSH(), v, 0, e->getPartition());
+    cell->getMeshVertices(v);
+    MElement* e = factory.create(cell->getTypeMSH(), v);
+    
     
-    if(cell->getDim() > 0 && coeff < 0) ne->revert(); // flip orientation
-    for(int i = 0; i < abs(coeff); i++) elements.push_back(ne);    
+    if(cell->getDim() > 0 && coeff < 0) e->revert(); // flip orientation
+    for(int i = 0; i < abs(coeff); i++) elements.push_back(e);    
 
     std::vector<double> coeffs (1,abs(coeff));
-    data[ne->getNum()] = coeffs;
+    data[e->getNum()] = coeffs;
   }
   
   int max[4];
-  for(int i = 0; i < 4; i++) max[i] = _model->getMaxElementaryNumber(i);
+  for(int i = 0; i < 4; i++) 
+    max[i] = this->getCellComplex()->getModel()->getMaxElementaryNumber(i);
   int entityNum = *std::max_element(max,max+4) + 1;
-  for(int i = 0; i < 4; i++) max[i] = _model->getMaxPhysicalNumber(i);
+  for(int i = 0; i < 4; i++) 
+    max[i] = this->getCellComplex()->getModel()->getMaxPhysicalNumber(i);
   int physicalNum = *std::max_element(max,max+4) + 1;
   setNum(physicalNum);
-  
+
   std::map<int, std::vector<MElement*> > entityMap;
   entityMap[entityNum] = elements;
   std::map<int, std::map<int, std::string> > physicalMap;
   std::map<int, std::string> physicalInfo;
-  physicalInfo[physicalNum]=getName();
+  physicalInfo[physicalNum] = getName();
   physicalMap[entityNum] = physicalInfo;
   
   // hide mesh
@@ -571,11 +532,14 @@ int Chain::createPGroup()
   //opt_view_show_element(0, GMSH_SET, 1);
   
   if(!data.empty()){
-    _model->storeChain(getDim(), entityMap, physicalMap);
-    _model->setPhysicalName(getName(), getDim(), physicalNum);
+    this->getCellComplex()->getModel()->storeChain(getDim(), 
+						   entityMap, physicalMap);
+    this->getCellComplex()->getModel()->setPhysicalName(getName(), 
+							getDim(), physicalNum);
 #if defined(HAVE_POST)    
-    // create PView for visualization
-    new PView(getName(), "ElementData", getGModel(), data, 0, 1);
+    // create PView for instant visualization
+    new PView(getName(), "ElementData", this->getCellComplex()->getModel(), 
+	      data, 0, 1);
 #endif
   }
    
@@ -583,46 +547,18 @@ int Chain::createPGroup()
 }
 
 
-void Chain::removeCell(Cell* cell) 
-{
-  citer it = _cells.find(cell);
-  if(it != _cells.end()){
-    (*it).second = 0;
-  }
-  return;
-}
-
-void Chain::addCell(Cell* cell, int coeff) 
-{
-  std::pair<citer,bool> insert = _cells.insert( std::make_pair( cell, coeff));
-  if(!insert.second && (*insert.first).second == 0){
-    (*insert.first).second = coeff; 
-  }
-  else if (!insert.second && (*insert.first).second != 0){
-    Msg::Debug("Error: invalid chain smoothening add!");
-  }
-  return;
-}
-
-bool Chain::hasCell(Cell* c)
-{
-  citer it = _cells.find(c);
-  if(it != _cells.end() && (*it).second != 0) return true;
-  return false;
-}
-   
-Cell* Chain::findCell(Cell* c)
-{
-  citer it = _cells.find(c);
-  if(it != _cells.end() && (*it).second != 0) return (*it).first;
-  return NULL;
-}
-
-int Chain::getCoeff(Cell* c)
-{
-  citer it = _cells.find(c);
-  if(it != _cells.end()) return (*it).second;
-  return 0;
+Chain::Chain(std::map<Cell*, int, Less_Cell>& chain, 
+	     CellComplex* cellComplex, int num,
+	     std::string name, int torsion)
+{  
+  _cells = chain;
+  if(!_cells.empty()) _dim = firstCell()->first->getDim();
+  else _dim = 0;
+  _name = name;
+  _num = num;
+  _cellComplex = cellComplex;
+  _torsion = torsion;
+  eraseNullCells();
 }
 
 void Chain::eraseNullCells()
@@ -635,13 +571,5 @@ void Chain::eraseNullCells()
   return;
 }
 
-void Chain::deImmuneCells()
-{
-  for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
-    Cell* cell = (*cit).first;
-    cell->setImmune(false);
-  }
-}
-
 #endif
 
diff --git a/Geo/Homology.h b/Geo/Homology.h
index 5f497f8bb70c1a33ce3108a18d204a2031be4f46..65ebee142cceffccba09c6b7fb6b9cb573278bd4 100644
--- a/Geo/Homology.h
+++ b/Geo/Homology.h
@@ -43,16 +43,26 @@ class Homology
   std::vector<GEntity*> _domainEntities;
   std::vector<GEntity*> _subdomainEntities;  
 
-  // physical group numbers of resulted generator chains in model
-  std::vector<int> _generators;
+  // use cell combining
+  bool _combine;
+  // use cell omit
+  bool _omit;
+  // use chain smoothning
+  bool _smoothen;
+
+  int _maxdomain;
+  int _maxnum;
 
   // file name to store the results
   std::string _fileName;
  
+  std::map<int, Chain*> _basisChains;
+
  public:
   
   Homology(GModel* model, std::vector<int> physicalDomain,
-	   std::vector<int> physicalSubdomain);
+	   std::vector<int> physicalSubdomain, 
+	   bool combine=true, bool omit=true, bool smoothen=true);
   ~Homology();
   
   // create a cell complex from a mesh in geometrical entities of Gmsh
@@ -67,9 +77,8 @@ class Homology
   void findGenerators(CellComplex* cellComplex=NULL);
   void findDualGenerators(CellComplex* cellComplex=NULL);
 
-
   // experimental
-  void findHomSequence();
+  void findHomSequence() {}
   void computeRanks() {}  
    
   // create a string describing the generator
@@ -82,11 +91,13 @@ class Homology
   // in _model, for debugging
   void storeCells(CellComplex* cellComplex, int dim);
 
+  GModel* getModel() const { return _model; }
+
 };
 
 // A class representing a chain.
-// Used to visualize and post-process generators of the homology spaces in Gmsh.
-class Chain{
+// Used to store generators of the homology spaces and visualize them in Gmsh.
+class Chain {
   
  private:
   // cells and their coefficients in this chain
@@ -97,7 +108,6 @@ class Chain{
   int _num;
   // cell complex this chain belongs to
   CellComplex* _cellComplex;
-  GModel* _model;
    
   // torsion coefficient
   int _torsion;
@@ -105,50 +115,22 @@ class Chain{
   int _dim;
   
  public:
-  Chain() {}
-  Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs, 
-	CellComplex* cellComplex, GModel* model,
-	std::string name="Chain", int torsion=0);
   Chain(std::map<Cell*, int, Less_Cell>& chain,
-        CellComplex* cellComplex, GModel* model,
+        CellComplex* cellComplex, int num,
         std::string name="Chain", int torsion=0);
-  Chain(Chain* chain){ 
-    chain->getCells(_cells);
-    _torsion = chain->getTorsion();
-    _name = chain->getName();
-    _cellComplex = chain->getCellComplex();
-    _dim = chain->getDim();
-    _model = chain->getGModel();
-    _num = chain->getNum();
-  }  
+
   typedef std::map<Cell*, int, Less_Cell>::iterator citer;
   citer firstCell() {return _cells.begin(); }
   citer lastCell() {return _cells.end(); }
 
-  ~Chain() {} 
-
-  // remove a cell from this chain
-  void removeCell(Cell* cell);
-   
-  // add a cell to this chain
-  void addCell(Cell* cell, int coeff);
-  
-  bool hasCell(Cell* c);
-  Cell* findCell(Cell* c);
-  int getCoeff(Cell* c);
-  
-  
   int getTorsion() const { return _torsion; }
   int getDim() const { return _dim; }
   CellComplex* getCellComplex() const { return _cellComplex; }
-  GModel* getGModel() const { return _model; }
-  void getCells(std::map<Cell*, int, Less_Cell> cells) const{ cells = _cells; }
+  void getCells(std::map<Cell*, int, Less_Cell> cells) const { 
+    cells = _cells; }
   
   // erase cells from the chain with zero coefficient
   void eraseNullCells();
-
-  // set all cell in chain not immune
-  void deImmuneCells();
   
   // number of cells in this chain 
   int getSize() const { return _cells.size();}
@@ -164,7 +146,7 @@ class Chain{
   // for debugging only
   int writeChainMSH(const std::string &name);
 
-  // create a physical group of this chain.
+  // create a Gmsh physical group from this chain.
   int createPGroup();
   
 };