Skip to content
Snippets Groups Projects
Select Git revision
  • be1117bc68966c7df6014f0fc4bc138fcb51bfca
  • master default protected
  • alphashapes
  • quadMeshingTools
  • cygwin_conv_path
  • macos_arm64
  • add-transfiniteautomatic-to-geo
  • patch_releases_4_10
  • HierarchicalHDiv
  • isuruf-master-patch-63355
  • hyperbolic
  • hexdom
  • hxt_update
  • jf
  • 1618-pythonocc-and-gmsh-api-integration
  • octreeSizeField
  • hexbl
  • alignIrregularVertices
  • getEdges
  • patch_releases_4_8
  • isuruf-master-patch-51992
  • gmsh_4_11_0
  • gmsh_4_10_5
  • gmsh_4_10_4
  • gmsh_4_10_3
  • gmsh_4_10_2
  • gmsh_4_10_1
  • gmsh_4_10_0
  • gmsh_4_9_5
  • gmsh_4_9_4
  • gmsh_4_9_3
  • gmsh_4_9_2
  • gmsh_4_9_1
  • gmsh_4_9_0
  • gmsh_4_8_4
  • gmsh_4_8_3
  • gmsh_4_8_2
  • gmsh_4_8_1
  • gmsh_4_8_0
  • gmsh_4_7_1
  • gmsh_4_7_0
41 results

CellComplex.cpp

Blame
  • Forked from gmsh / gmsh
    Source project has a limited visibility.
    CellComplex.cpp 29.84 KiB
    // Gmsh - Copyright (C) 1997-2016 C. Geuzaine, J.-F. Remacle
    //
    // See the LICENSE.txt file for license information. Please report all
    // bugs and problems to the public mailing list <gmsh@onelab.info>.
    //
    // Contributed by Matti Pellikka <matti.pellikka@microsoft.com>.
    
    #include "CellComplex.h"
    #include "MElement.h"
    #include "OS.h"
    
    double CellComplex::_patience = 10;
    
    CellComplex::CellComplex(GModel* model,
    			 std::vector<MElement*>& domainElements,
    			 std::vector<MElement*>& subdomainElements,
                             std::vector<MElement*>& nondomainElements,
                             std::vector<MElement*>& nonsubdomainElements,
                             std::vector<MElement*>& immuneElements,
                             bool saveOriginalComplex) :
      _model(model), _dim(0), _simplicial(true), _saveorig(saveOriginalComplex),
      _relative(false)
    {
      _smallestCell.second = -1.;
      _biggestCell.second = -1.;
      _deleteCount = 0;
      _createCount = 0;
      _insertCells(subdomainElements, 1);
      if(getSize(0) > 0) _relative = true;
      for(int i = 0; i < 4; i++) _numSubdomainCells[i] = getSize(i);
    
      _insertCells(domainElements, 0);
      for(int i = 0; i < 4; i++)
        _numRelativeCells[i] = getSize(i) - _numSubdomainCells[i];
    
      _removeCells(nonsubdomainElements, 1);
      _removeCells(nondomainElements, 0);
      _immunizeCells(immuneElements);
      int num = 0;
      for(int dim = 0; dim < 4; dim++){
        if(getSize(dim) != 0) _dim = dim;
        if(_saveorig) _ocells[dim] = _cells[dim];
        for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
          Cell* cell = *cit;
          cell->setNum(++num);
          cell->increaseGlobalNum();
          cell->saveCellBoundary();
        }
      }
    
      _reduced = false;
    
      Msg::Debug("Cells in domain:");
      Msg::Debug(" %d volumes, %d faces %d edges, and %d vertices",
                 getNumCells(3, 1), getNumCells(2, 1),
                 getNumCells(1, 1), getNumCells(0, 1));
      Msg::Debug("Cells in subdomain:");
      Msg::Debug(" %d volumes, %d faces %d edges, and %d vertices",
                 getNumCells(3, 2), getNumCells(2, 2),
                 getNumCells(1, 2), getNumCells(0, 2));
      Msg::Debug("Cells in relative domain:");
      Msg::Debug(" %d volumes, %d faces %d edges, and %d vertices",
                 getNumCells(3, 0), getNumCells(2, 0),
                 getNumCells(1, 0), getNumCells(0, 0));
    }
    
    bool CellComplex::_insertCells(std::vector<MElement*>& elements,
    			       int domain)
    {
      std::pair<Cell*, double> smallestElement[4];
      std::pair<Cell*, double> biggestElement[4];
      for(int i = 0; i < 4; i++) {
        smallestElement[i].second = -1.;
        biggestElement[i].second = -1.;
      }
      _dim = 0;
    
      double t1 = Cpu();
    
      for(unsigned int i=0; i < elements.size(); i++){
    
        MElement* element = elements.at(i);
        int dim = element->getDim();
        int type = element->getType();
        if(type == TYPE_POLYG || type == TYPE_POLYH) {
          Msg::Error("Mesh element type %d not implemented in homology solver", type);
        }
        if(type == TYPE_QUA || type == TYPE_HEX ||
           type == TYPE_PYR || type == TYPE_PRI)
          _simplicial = false;
        std::pair<Cell*, bool> maybeCell = Cell::createCell(element, domain);
        if(!maybeCell.second) {
          delete maybeCell.first;
          continue;
        }
    
        if(_dim < dim) _dim = dim;
        Cell* cell = maybeCell.first;
        std::pair<citer, bool> insert =
          _cells[cell->getDim()].insert(cell);
        if(!insert.second) {
          delete cell;
          cell = *(insert.first);
          if(domain) cell->setDomain(domain);
        }
        else _createCount++;
    
        if(domain == 0) {
          double size = fabs(element->getVolume());
          if(smallestElement[dim].second < 0. || smallestElement[dim].second > size)
            smallestElement[dim] = std::make_pair(cell, size);
          if(biggestElement[dim].second < 0. || biggestElement[dim].second < size)
            biggestElement[dim] = std::make_pair(cell, size);
        }
      }
      _smallestCell = smallestElement[_dim];
      _biggestCell = biggestElement[_dim];
    
      for (int dim = 3; dim > 0; dim--){
    
        double t2 = Cpu();
        if(t2-t1 > CellComplex::_patience && dim > 1) {
          if(domain == 0)
            Msg::Info(" ... creating domain %d-cells", dim);
          else if(domain == 1)
            Msg::Info(" ... creating subdomain %d-cells", dim);
        }
    
        for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
          Cell* cell = *cit;
          for(int i = 0; i < cell->getNumBdElements(); i++){
            std::pair<Cell*, bool> maybeCell = Cell::createCell(cell, i);
            if(!maybeCell.second) {
              delete maybeCell.first;
              continue;
            }
            Cell* newCell = maybeCell.first;
    	std::pair<citer, bool> insert =
    	  _cells[newCell->getDim()].insert(newCell);
    	if(!insert.second) {
    	  delete newCell;
    	  newCell = *(insert.first);
              if(domain) newCell->setDomain(domain);
    	}
            else _createCount++;
    	if(domain == 0) {
    	  int ori = cell->findBdCellOrientation(newCell, i);
    	  cell->addBoundaryCell( ori, newCell, true);
              if(_smallestCell.first == cell)
                _smallestCell = std::make_pair(newCell, _smallestCell.second);
              if(_biggestCell.first == cell)
                _biggestCell = std::make_pair(newCell, _biggestCell.second);
    	}
          }
        }
      }
      return true;
    }
    
    bool CellComplex::_removeCells(std::vector<MElement*>& elements,
    			       int domain)
    {
      if(elements.empty()) return true;
      Msg::Debug("Removing %d elements and their subcells from the cell complex.",
                 (int)elements.size());
      std::set<Cell*, Less_Cell> removed[4];
    
      for(unsigned int i=0; i < elements.size(); i++){
        MElement* element = elements.at(i);
        int type = element->getType();
        if(type == TYPE_PYR || type == TYPE_PRI ||
           type == TYPE_POLYG || type == TYPE_POLYH) {
          Msg::Error("Mesh element type %d not implemented in homology solver", type);
          return false;
        }
        Cell* cell = new Cell(element, domain);
        int dim = cell->getDim();
        citer cit = _cells[dim].find(cell);
        if(cit != lastCell(dim)) {
          removeCell(*cit);
          removed[dim].insert(cell);
        }
        else delete cell;
      }
    
      for (int dim = 3; dim > 0; dim--){
        for(citer cit = removed[dim].begin(); cit != removed[dim].end(); cit++){
          Cell* cell = *cit;
          for(int i = 0; i < cell->getNumBdElements(); i++){
    	Cell* newCell = new Cell(cell, i);
    
            citer cit2 = _cells[dim-1].find(newCell);
            if(cit2 != lastCell(dim-1)) {
              removeCell(*cit2);
              removed[dim-1].insert(newCell);
            }
            else delete newCell;
          }
        }
      }
      for (int dim = 3; dim >= 0; dim--){
        for(citer cit = removed[dim].begin(); cit != removed[dim].end(); cit++){
          delete *cit;
        }
      }
      Msg::Debug("Removed %d volumes, %d faces, %d edges, and %d vertices from the cell complex.",
                 (int)removed[3].size(), (int)removed[2].size(),
                 (int)removed[1].size(), (int)removed[0].size());
      return true;
    }
    
    bool CellComplex::_immunizeCells(std::vector<MElement*>& elements)
    {
      for(unsigned int i=0; i < elements.size(); i++){
        MElement* element = elements.at(i);
        Cell* cell = new Cell(element, 0);
        int dim = cell->getDim();
        citer cit = _cells[dim].find(cell);
        if(cit != lastCell(dim)) (*cit)->setImmune(true);
        delete cell;
      }
      return true;
    }
    
    CellComplex::~CellComplex()
    {
      for(int i = 0; i < 4; i++){
        for(citer cit = _cells[i].begin(); cit != _cells[i].end(); cit++){
          Cell* cell = *cit;
          delete cell;
          _deleteCount++;
        }
      }
    
      for(unsigned int i = 0; i < _removedcells.size(); i++) {
        delete _removedcells.at(i);
        _deleteCount++;
      }
    
      Msg::Debug("Total number of cells created: %d", _createCount);
      Msg::Debug("Total number of cells deleted: %d", _deleteCount);
    }
    
    void CellComplex::insertCell(Cell* cell)
    {
      std::pair<citer, bool> insertInfo = _cells[cell->getDim()].insert(cell);
      if(!insertInfo.second){
        Msg::Debug("Cell not inserted");
        Cell* oldCell = (*insertInfo.first);
        cell->printCell();
        oldCell->printCell();
      }
    }
    
    void CellComplex::removeCell(Cell* cell, bool other, bool del)
    {
      std::map<Cell*, short int, Less_Cell > coboundary;
      cell->getCoboundary(coboundary);
      std::map<Cell*, short int, Less_Cell > boundary;
      cell->getBoundary(boundary);
    
      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(std::map<Cell*, short int, Less_Cell>::iterator it =
            boundary.begin(); it != boundary.end(); it++){
        Cell* bdCell = (*it).first;
        bdCell->removeCoboundaryCell(cell, other);
      }
    
      int dim = cell->getDim();
      int erased = _cells[dim].erase(cell);
      if(relative()) {
        if(cell->inSubdomain()) _numSubdomainCells[dim] -= 1;
        else _numRelativeCells[dim] -= 1;
      }
      if(!erased) Msg::Debug("Tried to remove a cell from the cell complex \n");
      else if(!del) _removedcells.push_back(cell);
    }
    
    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*, short int, Less_Cell>::iterator cit = cells.begin();
          cit != cells.end(); cit++){
        Cell* cell = (*cit).first;
        citer it = Qset.find(cell);
        if(it == Qset.end()){
          Qset.insert(cell);
          Q.push(cell);
        }
      }
    }
    
    int CellComplex::coreduction(Cell* startCell, int omit,
    			     std::vector<Cell*>& omittedCells)
    {
      int coreductions = 0;
    
      std::queue<Cell*> Q;
      std::set<Cell*, Less_Cell> Qset;
    
      Q.push(startCell);
      Qset.insert(startCell);
    
      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();
        Qset.erase(s);
        if(s->getBoundarySize() == 1 &&
           inSameDomain(s, s->firstBoundary()->first) &&
           !s->getImmune() && !s->firstBoundary()->first->getImmune() &&
           abs(s->firstBoundary()->second.get()) < 2){
          s->getBoundary(bd_s);
          removeCell(s);
          bd_s.begin()->first->getCoboundary(cbd_c);
          enqueueCells(cbd_c, Q, Qset);
          removeCell(bd_s.begin()->first);
          if(bd_s.begin()->first->getDim() == omit){
    	omittedCells.push_back(bd_s.begin()->first);
          }
          coreductions++;
        }
        else if(s->getBoundarySize() == 0){
          s->getCoboundary(cbd_c);
          enqueueCells(cbd_c, Q, Qset);
        }
      }
      _reduced = true;
      return coreductions;
    }
    
    int CellComplex::reduction(int dim, int omit,
    			   std::vector<Cell*>& omittedCells)
    {
      if(dim < 1 || dim > 3) return 0;
    
      int numCells[4];
      for(int i = 0; i < 4; i++) numCells[i] = getSize(i);
    
      int count = 0;
    
      bool reduced = true;
      while (reduced){
    
        reduced = false;
        citer cit = firstCell(dim-1);
        while(cit != lastCell(dim-1)){
          Cell* cell = *cit;
          if(cell->getCoboundarySize() == 1 &&
             inSameDomain(cell, cell->firstCoboundary()->first) &&
             !cell->getImmune() && !cell->firstCoboundary()->first->getImmune() &&
             !cell->firstCoboundary()->first->getImmune() &&
             abs(cell->firstCoboundary()->second.get()) < 2){
    	cit++;
    	if(dim == omit){
    	  omittedCells.push_back(cell->firstCoboundary()->first);
    	}
    	removeCell(cell->firstCoboundary()->first);
    	removeCell(cell);
    	count++;
    	reduced = true;
          }
    
          if(getSize(dim) == 0 || getSize(dim-1) == 0) break;
          if(cit != lastCell(dim-1)) cit++;
        }
      }
      _reduced = true;
      Msg::Debug("Cell complex %d-reduction removed %dv, %df, %de, %dn", dim,
                 numCells[3]-getSize(3), numCells[2]-getSize(2),
                 numCells[1]-getSize(1), numCells[0]-getSize(0));
      return count;
    }
    
    int CellComplex::coreduction(int dim, int omit,
    			     std::vector<Cell*>& omittedCells)
    {
      if(dim < 1 || dim > 3) return 0;
    
      int numCells[4];
      for(int i = 0; i < 4; i++) numCells[i] = getSize(i);
    
      int count = 0;
    
      bool reduced = true;
      while (reduced){
    
        reduced = false;
        citer cit = firstCell(dim);
        while(cit != lastCell(dim)){
          Cell* cell = *cit;
          if(cell->getBoundarySize() == 1 &&
             inSameDomain(cell, cell->firstBoundary()->first) &&
             !cell->getImmune() && !cell->firstBoundary()->first->getImmune() &&
             abs(cell->firstBoundary()->second.get()) < 2) {
            ++cit;
    	if(dim-1 == omit){
    	  omittedCells.push_back(cell->firstBoundary()->first);
    	}
            removeCell(cell->firstBoundary()->first);
            removeCell(cell);
            count++;
            reduced = true;
          }
    
          if(getSize(dim) == 0 || getSize(dim-1) == 0) break;
          if(cit != lastCell(dim)) cit++;
        }
      }
      _reduced = true;
      Msg::Debug("Cell complex %d-coreduction removed %dv, %df, %de, %dn", dim,
                 numCells[3]-getSize(3), numCells[2]-getSize(2),
                 numCells[1]-getSize(1), numCells[0]-getSize(0));
      return count;
    }
    
    int CellComplex::getSize(int dim, bool orig) {
      if(dim == -1) {
        unsigned int size = 0;
        if(!orig) for(int i = 0; i < 4; i++) size += _cells[i].size();
        else for(int i = 0; i < 4; i++) size += _ocells[i].size();
        return size;
      }
      if(!orig) return _cells[dim].size();
      else return _ocells[dim].size();
    }
    
    int CellComplex::getDomain(Cell* cell, std::string& str)
    {
      int domain = 0;
      if(cell->inSubdomain()) {
        str = "subdomain";
        domain = 2;
      }
      if(!cell->inSubdomain()) {
        if(relative()) {
          str = "relative domain";
          domain = 0;
        }
        else {
          str = "domain";
          domain = 1;
        }
      }
      return domain;
    }
    
    Cell* CellComplex::_omitCell(Cell* cell, bool dual)
    {
      Msg::Debug("Omitting %d-cell from the cell complex", cell->getDim());
      removeCell(cell, false);
      std::vector<Cell*> omittedCells;
      omittedCells.push_back(cell);
      int count = 0;
    
      int numCells[4];
      for(int i = 0; i < 4; i++) numCells[i] = getSize(i);
    
      if(!dual) {
        for(int j = 3; j > 0; j--)
          count += reduction(j, cell->getDim(), omittedCells);
      }
      else {
        count += coreduction(cell, cell->getDim(), omittedCells);
        for(int j = 1; j <= getDim(); j++)
          count += coreduction(j, cell->getDim(), omittedCells);
      }
    
      CombinedCell* newcell = new CombinedCell(omittedCells);
      _createCount++;
    
      std::string domainstr = "";
      int domain = getDomain(cell, domainstr);
    
      Msg::Debug("Cell complex %d-omit removed %dv, %df, %de, %dn",
                 cell->getDim(),
                 numCells[3]-getSize(3), numCells[2]-getSize(2),
                 numCells[1]-getSize(1), numCells[0]-getSize(0));
      Msg::Debug(" - number of %d-cells left in %s: %d",
                 cell->getDim(), domainstr.c_str(),
                 getNumCells(cell->getDim(), domain));
    
      return newcell;
    }
    
    int CellComplex::reduceComplex(int combine, bool omit, bool homseq)
    {
      if(!getSize(0)) return 0;
    
      double t1 = Cpu();
      int count = 0;
      if(relative() && !homseq) removeSubdomain();
      std::vector<Cell*> empty;
      for(int i = 3; i > 0; i--) count = count + reduction(i, -1, empty);
    
      if(omit && !homseq){
    
        std::vector<Cell*> newCells;
    
        while (getSize(getDim()) != 0){
    
          citer cit = firstCell(getDim());
          Cell* cell = *cit;
    
          newCells.push_back(_omitCell(cell, false));
        }
    
        for(unsigned int i = 0; i < newCells.size(); i++){
          insertCell(newCells.at(i));
        }
      }
    
      double t2 = Cpu();
      if(t2-t1 > CellComplex::_patience) {
        Msg::Info(" .. %d volumes, %d faces, %d edges, and %d vertices",
                  getSize(3), getSize(2), getSize(1), getSize(0));
      }
    
      if(combine > 0) this->combine(3);
    
      if(combine > 2) for(int i = 3; i > 0; i--) reduction(i, -1, empty);
      else if(combine > 1) reduction(2, -1, empty);
    
      if(combine > 0) this->combine(2);
    
      if(combine > 2) for(int i = 3; i > 0; i--) reduction(i, -1, empty);
      else if(combine > 1) reduction(1, -1, empty);
    
      if(combine > 0) this->combine(1);
    
      if(combine > 2) for(int i = 3; i > 0; i--) reduction(i, -1, empty);
      else if(combine > 1) reduction(0, -1, empty);
    
      _reduced = true;
      return count;
    }
    
    void CellComplex::removeSubdomain()
    {
      std::vector<Cell*> toRemove;
      for(int i = 0; i < 4; i++){
        for(citer cit = firstCell(i); cit != lastCell(i); ++cit){
          Cell *cell = *cit;
          if(cell->inSubdomain()) toRemove.push_back(cell);
        }
      }
      for(unsigned int i = 0; i < toRemove.size(); i++) removeCell(toRemove[i]);
      _reduced = true;
    }
    
    void CellComplex::removeCells(int dim)
    {
      if(dim < 0 || dim > 3) return;
      std::vector<Cell*> toRemove;
      for(citer cit = firstCell(dim); cit != lastCell(dim); ++cit){
        toRemove.push_back(*cit);
      }
      for(unsigned int i = 0; i < toRemove.size(); i++) removeCell(toRemove[i]);
      _reduced = true;
    }
    
    int CellComplex::coreduceComplex(int combine, bool omit, int heuristic)
    {
      if(!getSize(0)) return 0;
    
      double t1 = Cpu();
    
      int count = 0;
      if(relative()) removeSubdomain();
      std::vector<Cell*> empty;
      for(int dim = 0; dim < 4; dim++){
        citer cit = firstCell(dim);
        while(cit != lastCell(dim)){
          Cell* cell = *cit;
          int count =+ coreduction(cell, -1, empty);
          if(count != 0) break;
          cit++;
        }
      }
    
      for(int j = 1; j <= getDim(); j++)
        count += coreduction(j, -1, empty);
    
      if(omit){
    
        std::vector<Cell*> newCells;
        while (getSize(0) != 0){
    
          citer cit = firstCell(0);
          Cell* cell = *cit;
    
          if(heuristic == -1 && _smallestCell.second != 0. &&
             hasCell(_smallestCell.first)) {
            Msg::Debug("Omitted a cell in the smallest mesh element with volume %g",
                       _smallestCell.second);
            cell = _smallestCell.first;
          }
          else if(heuristic == 1 && _biggestCell.second != 0. &&
                  hasCell(_biggestCell.first)) {
            Msg::Debug("Omitted a cell in the biggest mesh element with volume %g",
                       _biggestCell.second);
            cell = _biggestCell.first;
          }
    
          newCells.push_back(_omitCell(cell, true));
        }
        for(unsigned int i = 0; i < newCells.size(); i++){
          insertCell(newCells.at(i));
        }
    
      }
    
      double t2 = Cpu();
      if(t2-t1 > CellComplex::_patience) {
        Msg::Info(" .. %d volumes, %d faces, %d edges, and %d vertices",
                  getSize(3), getSize(2), getSize(1), getSize(0));
      }
    
      if(combine > 0) this->cocombine(0);
    
      if(combine > 2) for(int i = 1; i < 4; i++) coreduction(i, -1, empty);
      else if(combine > 1) coreduction(1, -1, empty);
    
      if(combine > 0) this->cocombine(1);
    
      if(combine > 2)  for(int i = 1; i < 4; i++) coreduction(i, -1, empty);
      else if(combine > 1) coreduction(2, -1, empty);
    
      if(combine > 0) this->cocombine(2);
    
      if(combine > 2) for(int i = 1; i < 4; i++) coreduction(i, -1, empty);
      else if(combine > 1) coreduction(3, -1, empty);
    
      coherent();
    
      _reduced = true;
      return count;
    }
    
    void CellComplex::bettiReduceComplex()
    {
      reduceComplex(3, true);
      for(int i = 1; i <= 3; i++) cocombine(i-1);
      return;
    }
    
    int CellComplex::combine(int dim)
    {
      if(dim < 1 || dim > 3) return 0;
    
      int numCells[4];
      for(int i = 0; i < 4; i++) numCells[i] = getSize(i);
    
      double t1 = Cpu();
    
      std::queue<Cell*> Q;
      std::set<Cell*, Less_Cell> Qset;
      std::map<Cell*, short int, Less_Cell> bd_c;
      int count = 0;
    
      for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
    
        double t2 = Cpu();
        if(t2-t1 > CellComplex::_patience) {
          t1 = Cpu();
          Msg::Info(" ... %d volumes, %d faces, %d edges, and %d vertices",
                    getSize(3), getSize(2), getSize(1), getSize(0));
        }
    
        Cell* cell = *cit;
        cell->getBoundary(bd_c);
        enqueueCells(bd_c, Q, Qset);
    
        while(Q.size() != 0){
          Cell* s = Q.front();
          Q.pop();
    
          if(s->getCoboundarySize() == 2 && !s->getImmune()){
    	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) &&
               c1->getImmune() == c2->getImmune() ) {
    
              removeCell(s, true, false);
    
              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));
              _createCount++;
              removeCell(c1, true, c1->isCombined());
              removeCell(c2, true, c2->isCombined());
              insertCell(newCell);
    
              cit = firstCell(dim);
              count++;
    
              if(c1->isCombined()) {
                delete c1;
                _deleteCount++;
              }
              if(c2->isCombined()) {
                delete c2;
                _deleteCount++;
              }
    
            }
          }
          Qset.erase(s);
        }
      }
    
      Msg::Debug("Cell complex %d-combine removed %dv, %df, %de, %dn", dim,
                 numCells[3]-getSize(3), numCells[2]-getSize(2),
                 numCells[1]-getSize(1), numCells[0]-getSize(0));
      _reduced = true;
      return count;
    }
    
    int CellComplex::cocombine(int dim)
    {
      if(dim < 0 || dim > 2) return 0;
    
      int numCells[4];
      for(int i = 0; i < 4; i++) numCells[i] = getSize(i);
    
      double t1 = Cpu();
    
      std::queue<Cell*> Q;
      std::set<Cell*, Less_Cell> Qset;
      std::map<Cell*, short int, Less_Cell> cbd_c;
      int count = 0;
    
      for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
    
        double t2 = Cpu();
        if(t2-t1 > CellComplex::_patience) {
          t1 = Cpu();
          Msg::Info(" ... %d volumes, %d faces, %d edges, and %d vertices",
                    getSize(3), getSize(2), getSize(1), getSize(0));
        }
    
        Cell* cell = *cit;
    
        cell->getCoboundary(cbd_c);
        enqueueCells(cbd_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.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)
               && c1->getImmune() == c2->getImmune()){
    
              removeCell(s, true, false);
    
              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 );
              _createCount++;
              removeCell(c1, true, c1->isCombined());
              removeCell(c2, true, c2->isCombined());
              insertCell(newCell);
    
              cit = firstCell(dim);
              count++;
    
              if(c1->isCombined()) {
                delete c1;
                _deleteCount++;
              }
              if(c2->isCombined()) {
                delete c2;
                _deleteCount++;
              }
    
            }
          }
          Qset.erase(s);
        }
      }
    
      Msg::Debug("Cell complex %d-cocombine removed %dv, %df, %de, %dn", dim,
                 numCells[3]-getSize(3), numCells[2]-getSize(2),
                 numCells[1]-getSize(1), numCells[0]-getSize(0));
      _reduced = true;
      return count;
    }
    
    bool CellComplex::coherent()
    {
      bool coherent = true;
      for(int i = 0; i < 4; i++){
        for(citer cit = firstCell(i); cit != lastCell(i); cit++){
          Cell* cell = *cit;
          std::map<Cell*, short int, Less_Cell> boundary;
          cell->getBoundary(boundary);
          for(std::map<Cell*, short int, Less_Cell>::iterator it = boundary.begin();
    	  it != boundary.end(); it++){
            Cell* bdCell = (*it).first;
            int ori = (*it).second;
            citer cit = _cells[bdCell->getDim()].find(bdCell);
            if(cit == lastCell(bdCell->getDim())){
              Msg::Debug("Boundary cell not in cell complex! Boundary removed");
              cell->removeBoundaryCell(bdCell, false);
              coherent = false;
            }
            if(!bdCell->hasCoboundary(cell)){
              Msg::Debug("Incoherent boundary/coboundary pair! Fixed");
    	  bdCell->addCoboundaryCell(ori, cell, false);
              coherent = false;
            }
    
          }
          std::map<Cell*, short int, Less_Cell> coboundary;
          cell->getCoboundary(coboundary);
          for(std::map<Cell*, short int, Less_Cell>::iterator it = coboundary.begin();
    	  it != coboundary.end(); it++){
            Cell* cbdCell = (*it).first;
            int ori = (*it).second;
            citer cit = _cells[cbdCell->getDim()].find(cbdCell);
            if(cit == lastCell(cbdCell->getDim())){
              Msg::Debug("Coboundary cell not in cell complex! Coboundary removed");
              cell->removeCoboundaryCell(cbdCell, false);
              coherent = false;
            }
            if(!cbdCell->hasBoundary(cell)){
              Msg::Debug("Incoherent coboundary/boundary pair! Fixed");
    	  cbdCell->addBoundaryCell(ori, cell, false);
              coherent = false;
            }
    
          }
    
        }
      }
      return coherent;
    }
    
    bool CellComplex::hasCell(Cell* cell, bool orig)
    {
      citer cit;
      if(!orig) cit = _cells[cell->getDim()].find(cell);
      else cit = _ocells[cell->getDim()].find(cell);
      if( cit == lastCell(cell->getDim(), orig) ) return false;
      else return true;
    }
    
    void CellComplex::getCells(std::set<Cell*, Less_Cell>& cells,
                               int dim, int domain){
      cells.clear();
      for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
        Cell* cell = *cit;
        if( (domain == 0 && !cell->inSubdomain()) || domain == 1
    	|| (domain == 2 && cell->inSubdomain()) ){
          cells.insert(cell);
        }
      }
    }
    
    int CellComplex::getNumCells(int dim, int domain)
    {
      if(domain == 0) return _numRelativeCells[dim];
      else if(domain == 1) return getSize(dim);
      else if(domain == 2) return _numSubdomainCells[dim];
      return 0;
    }
    
    Cell* CellComplex::getACell(int dim, int domain)
    {
      int num = getNumCells(dim, domain);
      if(num < 0) Msg::Debug("Domain cell counts not in sync.");
    
      if(num <= 0) {
        if(domain == 0) Msg::Warning("%d cells in relative domain", num);
        else if(domain == 1) Msg::Warning("%d cells in domain", num);
        else if(domain == 2) Msg::Warning("%d cells in subdomain", num);
        return NULL;
      }
    
      for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
        Cell* cell = *cit;
        if( (domain == 1) ||
            (domain == 0 && !cell->inSubdomain()) ||
            (domain == 2 && cell->inSubdomain()) ) return cell;
      }
      Msg::Debug("Domain cell counts not in sync.");
      return NULL;
    }
    
    bool CellComplex::restoreComplex()
    {
      if(_saveorig){
    
        for(unsigned int i = 0; i < _removedcells.size(); i++) {
          Cell* cell = _removedcells.at(i);
          if(cell->isCombined()) {
            delete cell;
            _deleteCount++;
          }
        }
        _removedcells.clear();
    
        for(int i = 0; i < 4; i++){
    
          for(citer cit = _cells[i].begin(); cit != _cells[i].end(); cit++){
            Cell* cell = *cit;
            if(cell->isCombined()) {
              delete cell;
              _deleteCount++;
            }
          }
    
          _cells[i] = _ocells[i];
          for(citer cit = firstCell(i); cit != lastCell(i); cit++){
    	Cell* cell = *cit;
    	cell->restoreCellBoundary();
            if(relative()) {
              if(cell->inSubdomain()) _numSubdomainCells[i] += 1;
              else _numRelativeCells[i] += 1;
            }
          }
        }
    
        Msg::Info("Restored Cell Complex:");
        Msg::Info("%d volumes, %d faces, %d edges, and %d vertices",
                   getSize(3), getSize(2), getSize(1), getSize(0));
        _reduced = false;
        return true;
      }
      else {
        Msg::Error("Cannot restore cell complex");
        return false;
      }
    }
    
    void CellComplex::printComplex(int dim)
    {
      if(getSize(dim) == 0)
        Msg::Info("Cell complex dimension %d is empty", dim);
      for (citer cit = firstCell(dim); cit != lastCell(dim); cit++){
        Cell* cell = *cit;
        cell->printCell();
        cell->printBoundary();
        cell->printCoboundary();
        printf("--- \n");
      }
    }
    
    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;
    }