diff --git a/Geo/Cell.cpp b/Geo/Cell.cpp index b7a7a662db3a392e61b97e9b77d28c8bef15e3fe..cd238e73fcc85dc84582008c82a7f1641a987670 100755 --- a/Geo/Cell.cpp +++ b/Geo/Cell.cpp @@ -17,7 +17,7 @@ bool Less_Cell::operator()(const Cell* c1, const Cell* c2) const { if(c1->getNumVertices() != c2->getNumVertices()){ return (c1->getNumVertices() < c2->getNumVertices()); } - + // "natural number" -like ordering // (the number of a vertice corresponds a digit) @@ -115,7 +115,8 @@ void Cell::restoreCell(){ _immune = false; } -bool Cell::addBoundaryCell(int orientation, Cell* cell) { +bool Cell::addBoundaryCell(int orientation, Cell* cell, bool org) { + if(org) _obd.insert( std::make_pair(cell, orientation ) ); biter it = _boundary.find(cell); if(it != _boundary.end()){ (*it).second = (*it).second + orientation; @@ -130,7 +131,8 @@ bool Cell::addBoundaryCell(int orientation, Cell* cell) { return true; } -bool Cell::addCoboundaryCell(int orientation, Cell* cell) { +bool Cell::addCoboundaryCell(int orientation, Cell* cell, bool org) { + if(org) _ocbd.insert( std::make_pair(cell, orientation ) ); biter it = _coboundary.find(cell); if(it != _coboundary.end()){ (*it).second = (*it).second + orientation; @@ -198,43 +200,28 @@ void Cell::makeDualCell(){ _dim = 3-_dim; } -void Cell::printBoundary() { - for(biter it = _boundary.begin(); it != _boundary.end(); it++){ +void Cell::printBoundary(bool org) { + for(biter it = firstBoundary(org); it != lastBoundary(org); it++){ printf("Boundary cell orientation: %d ", (*it).second); Cell* cell2 = (*it).first; cell2->printCell(); } - if(_boundary.empty()) printf("Cell boundary is empty. \n"); -} -void Cell::printCoboundary() { - for(biter it = _coboundary.begin(); it != _coboundary.end(); it++){ - printf("Coboundary cell orientation: %d, ", (*it).second); - Cell* cell2 = (*it).first; - cell2->printCell(); - if(_coboundary.empty()) printf("Cell coboundary is empty. \n"); - } -} - -void Cell::printOrgBd() { - for(biter it = _obd.begin(); it != _obd.end(); it++){ - printf("Boundary cell orientation: %d ", (*it).second); - Cell* cell2 = (*it).first; - cell2->printCell(); + if(firstBoundary(org) == lastBoundary(org)){ + printf("Cell boundary is empty. \n"); } - if(_obd.empty()) printf("Cell boundary is empty. \n"); } -void Cell::printOrgCbd() { - for(biter it = _ocbd.begin(); it != _ocbd.end(); it++){ +void Cell::printCoboundary(bool org) { + for(biter it = firstCoboundary(org); it != lastCoboundary(org); it++){ printf("Coboundary cell orientation: %d, ", (*it).second); Cell* cell2 = (*it).first; cell2->printCell(); - if(_ocbd.empty()) printf("Cell coboundary is empty. \n"); + if(firstCoboundary(org) == lastCoboundary(org)){ + printf("Cell coboundary is empty. \n"); + } } -} - +} -CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell() { - +CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell() { // use "smaller" cell as c2 if(c1->getNumVertices() < c2->getNumVertices()){ Cell* temp = c1; @@ -244,27 +231,23 @@ CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell() { _index = c1->getIndex(); _dim = c1->getDim(); - _num = c1->getNum(); _inSubdomain = c1->inSubdomain(); _onDomainBoundary = c1->onDomainBoundary(); _combined = true; _image = NULL; - _v.reserve(c1->getNumVertices() + c2->getNumVertices()); + // vertices _vs.reserve(c1->getNumVertices() + c2->getNumVertices()); - - - for(int i = 0; i < c1->getNumVertices(); i++) _v.push_back(c1->getVertex(i)); + for(int i = 0; i < c1->getNumVertices(); i++){ + _vs.push_back(c1->getSortedVertex(i)); + } + std::vector<int> v; for(int i = 0; i < c2->getNumVertices(); i++){ - if(!this->hasVertex(c2->getVertex(i)->getNum())){ - _v.push_back(c2->getVertex(i)); + if(!this->hasVertex(c2->getSortedVertex(i))){ + v.push_back(c2->getSortedVertex(i)); } } - - // sorted vertices - for(unsigned int i = 0; i < _v.size(); i++){ - _vs.push_back(_v.at(i)->getNum()); - } + for(unsigned int i = 0; i < v.size(); i++) _vs.push_back(v.at(i)); std::sort(_vs.begin(), _vs.end()); // cells diff --git a/Geo/Cell.h b/Geo/Cell.h index b5df67f2fce115e49cd5e5d24b4d4a0bc9534089..85a916da201d22fa77a21dbd274f557f02258667 100644 --- a/Geo/Cell.h +++ b/Geo/Cell.h @@ -78,7 +78,10 @@ class Cell // sorted vertices of this cell (used for ordering of the cells) std::vector<int> _vs; - + + virtual MVertex* getVertex(int vertex) const { + return _image->getVertex(vertex); } + public: Cell() : _combined(false), _index(0), _immune(false), _image(NULL), @@ -92,10 +95,10 @@ class Cell virtual int getDim() const { return _dim; }; virtual int getIndex() const { return _index; }; virtual void setIndex(int index) { _index = index; }; - virtual int getNum() { return _image->getNum(); } - virtual int getType() { return _image->getType(); } - virtual int getTypeForMSH() { return _image->getTypeForMSH(); } - virtual int getPartition() { return _image->getPartition(); } + virtual int getNum() const { return _image->getNum(); } + virtual int getType() const { return _image->getType(); } + virtual int getTypeForMSH() const { return _image->getTypeForMSH(); } + virtual int getPartition() const { return _image->getPartition(); } virtual void setImmune(bool immune) { _immune = immune; }; virtual bool getImmune() const { return _immune; }; virtual void setDeleteImage(bool deleteImage) { @@ -104,8 +107,6 @@ class Cell // get the number of vertices this cell has virtual int getNumVertices() const { return _image->getNumVertices(); } - virtual MVertex* getVertex(int vertex) const { - return _image->getVertex(vertex); } virtual int getSortedVertex(int vertex) const { return _vs.at(vertex); } // restores the cell information to its original state before reduction @@ -116,23 +117,29 @@ class Cell // (co)boundary cell iterator typedef std::map<Cell*, int, Less_Cell>::iterator biter; - biter firstBoundary(){ return _boundary.begin(); } - biter lastBoundary(){ return _boundary.end(); } - biter firstCoboundary(){ return _coboundary.begin(); } - biter lastCoboundary(){ return _coboundary.end(); } + biter firstBoundary(bool org=false){ + return org ? _obd.begin() : _boundary.begin(); } + biter lastBoundary(bool org=false){ + return org ? _obd.end() : _boundary.end(); } + biter firstCoboundary(bool org=false){ + return org ? _ocbd.begin() : _coboundary.begin(); } + biter lastCoboundary(bool org=false){ + return org ? _ocbd.end() : _coboundary.end(); } virtual int getBoundarySize() { return _boundary.size(); } virtual int getCoboundarySize() { return _coboundary.size(); } // get the cell boundary - virtual void getBoundary(std::map<Cell*, int, Less_Cell >& boundary){ - boundary = _boundary; } - virtual void getCoboundary(std::map<Cell*, int, Less_Cell >& coboundary){ - coboundary = _coboundary; } + virtual void getBoundary(std::map<Cell*, int, Less_Cell >& boundary, + bool org=false){ + org ? boundary = _obd : boundary = _boundary; } + virtual void getCoboundary(std::map<Cell*, int, Less_Cell >& coboundary, + bool org=false){ + org ? coboundary = _ocbd : coboundary = _coboundary; } // add (co)boundary cell - virtual bool addBoundaryCell(int orientation, Cell* cell); - virtual bool addCoboundaryCell(int orientation, Cell* cell); + virtual bool addBoundaryCell(int orientation, Cell* cell, bool org=false); + virtual bool addCoboundaryCell(int orientation, Cell* cell, bool org=false); // remove (co)boundary cell virtual int removeBoundaryCell(Cell* cell, bool other=true); @@ -150,20 +157,8 @@ class Cell // print cell info virtual void printCell(); - virtual void printBoundary(); - virtual void printCoboundary(); - - // original (co)boundary - virtual void addOrgBdCell(int orientation, Cell* cell) { - _obd.insert( std::make_pair(cell, orientation ) ); }; - virtual void addOrgCbdCell(int orientation, Cell* cell) { - _ocbd.insert( std::make_pair(cell, orientation ) ); }; - virtual void getOrgBd( std::map<Cell*, int, Less_Cell >& obd) { - obd = _obd; } - virtual void getOrgCbd( std::map<Cell*, int, Less_Cell >& obdc) { - obdc = _ocbd; } - virtual void printOrgBd(); - virtual void printOrgCbd(); + virtual void printBoundary(bool org=false); + virtual void printCoboundary(bool org=false); virtual bool inSubdomain() const { return _inSubdomain; } virtual void setInSubdomain(bool subdomain) { _inSubdomain = subdomain; } @@ -214,28 +209,26 @@ class Cell class CombinedCell : public Cell{ private: - // vertices - std::vector<MVertex*> _v; - // sorted list of all vertices + // sorted list of vertices std::vector<int> _vs; // list of cells this cell is a combination of std::list< std::pair<int, Cell*> > _cells; - int _num; + + MVertex* getVertex(int vertex) const { return NULL; } public: CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co=false); ~CombinedCell() {} - int getNum() { return _num; } - int getType() { return 0; } - int getTypeForMSH() { return 0; } - int getPartition() { return 0; } + int getNum() const { return 0; } + int getType() const { return 0; } + int getTypeForMSH() const { return 0; } + int getPartition() const { return 0; } - int getNumVertices() const { return _v.size(); } - MVertex* getVertex(int vertex) const { return _v.at(vertex); } + int getNumVertices() const { return _vs.size(); } int getSortedVertex(int vertex) const { return _vs.at(vertex); } - + void getCells(std::list< std::pair<int, Cell*> >& cells) { cells = _cells; } int getNumCells() {return _cells.size();} diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp index e7f059b0b63b463f692032d198ac7c3ea1154cf1..9daf3e83ecabfb113384220de28de14900bf8240 100644 --- a/Geo/CellComplex.cpp +++ b/Geo/CellComplex.cpp @@ -217,19 +217,15 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){ Cell* oldCell = *(insertInfo.first); if(!subdomain && !boundary){ int ori = cell->getFacetOri(oldCell); - oldCell->addCoboundaryCell( ori, cell ); - oldCell->addOrgCbdCell( ori, cell ); - cell->addBoundaryCell( ori, oldCell); - cell->addOrgBdCell( ori, oldCell); + oldCell->addCoboundaryCell( ori, cell, true); + cell->addBoundaryCell( ori, oldCell, true); } } else if(!subdomain && !boundary) { int ori = cell->getFacetOri(vertices); - cell->addBoundaryCell( ori, newCell ); - cell->addOrgBdCell( ori, newCell ); - newCell->addCoboundaryCell( ori, cell); - newCell->addOrgCbdCell( ori, cell); - } + cell->addBoundaryCell( ori, newCell, true); + newCell->addCoboundaryCell( ori, cell, true); + } } } } @@ -327,9 +323,6 @@ int CellComplex::coreduction(Cell* generator, int omitted){ Cell* s; int round = 0; while( !Q.empty() ){ - round++; - //printf("%d ", round); - s = Q.front(); Q.pop(); removeCellQset(s, Qset); @@ -349,10 +342,7 @@ int CellComplex::coreduction(Cell* generator, int omitted){ s->getCoboundary(cbd_c); enqueueCells(cbd_c, Q, Qset); } - - } - //printf("Coreduction: %d loops with %d coreductions\n", round, coreductions); return coreductions; } @@ -368,8 +358,6 @@ int CellComplex::reduction(int dim, int omitted){ reduced = false; citer cit = firstCell(dim-1); while(cit != lastCell(dim-1)){ - - Cell* cell = *cit; if( cell->getCoboundarySize() == 1 && inSameDomain(cell, cell->firstCoboundary()->first)){ @@ -381,15 +369,12 @@ int CellComplex::reduction(int dim, int omitted){ removeCell(cell); count++; reduced = true; - } if(getSize(dim) == 0 || getSize(dim-1) == 0) break; cit++; } - } - return count; } /* @@ -482,7 +467,7 @@ int CellComplex::reduceComplex(){ double t1 = Cpu(); Msg::Debug("Cell complex before reduction: %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)); int count = 0; for(int i = 3; i > 0; i--) count = count + reduction(i); @@ -503,23 +488,11 @@ int CellComplex::reduceComplex(){ _store.push_back(omittedCells); _store.at(omitted-1).insert(cell); for(int j = 3; j > 0; j--) reduction(j, omitted); - - } - - + double t2 = Cpu(); Msg::Debug("Cell complex after reduction: %d volumes, %d faces, %d edges and %d vertices (%g s).\n", getSize(3), getSize(2), getSize(1), getSize(0), t2 - t1); - - /* - for(int i = 0; i < _store.size(); i++){ - printf("omitted generator %d: \n", i+1); - for(citer cit = _store.at(i).begin(); cit != _store.at(i).end(); cit++){ - (*cit)->printCell(); - } - }*/ - return 0; } @@ -532,7 +505,6 @@ void CellComplex::removeSubdomain(){ if(cell->inSubdomain()) { removeCell(cell); ++cit; - //cit = firstCell(i); } } for(citer cit = firstCell(i); cit != lastCell(i); cit++){ @@ -580,19 +552,9 @@ int CellComplex::coreduceComplex(){ _store.at(omitted-1).insert(cell); coreduction(cell, omitted); } - - /* - for(citer cit = generatorCells.begin(); cit != generatorCells.end(); cit++){ - Cell* cell = *cit; - cell->clearBoundary(); - cell->clearCoboundary(); - _cells[0].insert(cell); - } - */ - Msg::Debug("Cell complex after coreduction: %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)); return omitted; } @@ -616,7 +578,8 @@ void CellComplex::computeBettiNumbers(){ } } Msg::Debug("Cell complex Betti numbers: \nH0 = %d \nH1 = %d \nH2 = %d \nH3 = %d \n", - getBettiNumber(0), getBettiNumber(1), getBettiNumber(2), getBettiNumber(3)); + getBettiNumber(0), getBettiNumber(1), + getBettiNumber(2), getBettiNumber(3)); return; } @@ -654,7 +617,8 @@ int CellComplex::cocombine(int dim){ if(!(*c1 == *c2) && abs(or1) == abs(or2) && inSameDomain(s, c1) && inSameDomain(s, c2) - && c1->getNumVertices() < getSize(dim) // heuristics for mammoth cell birth control + && c1->getNumVertices() < getSize(dim) + // heuristics for mammoth cell birth control && c2->getNumVertices() < getSize(dim)){ removeCell(s); @@ -694,7 +658,6 @@ int CellComplex::combine(int dim){ std::queue<Cell*> Q; std::set<Cell*, Less_Cell> Qset; - //std::map<Cell*, int, Less_Cell> cbd_c; std::map<Cell*, int, Less_Cell> bd_c; int count = 0; @@ -719,7 +682,8 @@ int CellComplex::combine(int dim){ if(!(*c1 == *c2) && abs(or1) == abs(or2) && inSameDomain(s, c1) && inSameDomain(s, c2) - && c1->getNumVertices() < getSize(dim) // heuristics for mammoth cell birth control + && c1->getNumVertices() < getSize(dim) + // heuristics for mammoth cell birth control && c2->getNumVertices() < getSize(dim)){ removeCell(s); @@ -847,24 +811,24 @@ bool CellComplex::checkCoherence(){ Cell* cell = *cit; std::map<Cell*, int, Less_Cell> boundary; cell->getBoundary(boundary); - for(std::map<Cell*, int, Less_Cell>::iterator it = boundary.begin(); + for(Cell::biter 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("Warning! Boundary cell not in cell complex! Boundary removed. \n"); - //printf(" "); cell->printCell(); - //printf(" "); bdCell->printCell(); + //cell->printCell(); + //bdCell->printCell(); cell->removeBoundaryCell(bdCell); coherent = false; } if(!bdCell->hasCoboundary(cell)){ Msg::Debug("Warning! Incoherent boundary/coboundary pair! Fixed. \n"); - //printf(" "); cell->printCell(); - //printf(" "); cell->printBoundary(); - //printf(" "); bdCell->printCell(); - //printf(" "); bdCell->printCoboundary(); + /*cell->printCell(); + cell->printBoundary(); + bdCell->printCell(); + bdCell->printCoboundary();*/ bdCell->addCoboundaryCell(ori, cell); coherent = false; } @@ -872,24 +836,24 @@ bool CellComplex::checkCoherence(){ } std::map<Cell*, int, Less_Cell> coboundary; cell->getCoboundary(coboundary); - for(std::map<Cell*, int, Less_Cell>::iterator it = coboundary.begin(); + for(Cell::biter 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("Warning! Coboundary cell not in cell complex! Coboundary removed. \n"); - //printf(" "); cell->printCell(); - //printf(" "); cbdCell->printCell(); + cell->printCell(); + cbdCell->printCell(); cell->removeCoboundaryCell(cbdCell); coherent = false; } if(!cbdCell->hasBoundary(cell)){ Msg::Debug("Warning! Incoherent coboundary/boundary pair! Fixed. \n"); - //printf(" "); cell->printCell(); - //printf(" "); cell->printCoboundary(); - //printf(" "); cbdCell->printCell(); - //printf(" "); cbdCell->printBoundary(); + cell->printCell(); + cell->printCoboundary(); + cbdCell->printCell(); + cbdCell->printBoundary(); cbdCell->addBoundaryCell(ori, cell); coherent = false; } @@ -931,8 +895,10 @@ bool CellComplex::swapSubdomain(){ for(int i = 0; i < 4; i++){ for(citer cit = firstCell(i); cit != lastCell(i); cit++){ Cell* cell = *cit; - if(cell->onDomainBoundary() && cell->inSubdomain()) cell->setInSubdomain(false); - else if(cell->onDomainBoundary() && !cell->inSubdomain()) cell->setInSubdomain(true); + if(cell->onDomainBoundary() + && cell->inSubdomain()) cell->setInSubdomain(false); + else if(cell->onDomainBoundary() + && !cell->inSubdomain()) cell->setInSubdomain(true); } } return true; diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp index 7edb78d8c41eb27a6ae888cd2cd856f3bf578117..cbd580629870472829bd3cece777194001d2187e 100644 --- a/Geo/ChainComplex.cpp +++ b/Geo/ChainComplex.cpp @@ -33,7 +33,8 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){ int index = 1; // ignore subdomain cells - for(std::set<Cell*, Less_Cell>::iterator cit = cellComplex->firstCell(dim); cit != cellComplex->lastCell(dim); cit++){ + for(CellComplex::citer cit = cellComplex->firstCell(dim); + cit != cellComplex->lastCell(dim); cit++){ Cell* cell = *cit; cell->setIndex(index); index++; @@ -44,7 +45,8 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){ } index = 1; if(dim > 0){ - for(std::set<Cell*, Less_Cell>::iterator cit = cellComplex->firstCell(dim-1); cit != cellComplex->lastCell(dim-1); cit++){ + for(CellComplex::citer cit = cellComplex->firstCell(dim-1); + cit != cellComplex->lastCell(dim-1); cit++){ Cell* cell = *cit; cell->setIndex(index); index++; @@ -70,7 +72,6 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){ mpz_init(elem); _HMatrix[dim] = create_gmp_matrix_zero(rows, cols); - //printMatrix(_HMatrix[dim]); for( std::set<Cell*, Less_Cell>::iterator cit = cellComplex->firstCell(dim); cit != cellComplex->lastCell(dim); cit++){ @@ -81,12 +82,13 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){ Cell* bdCell = (*it).first; if(!bdCell->inSubdomain()){ int old_elem = 0; - //printf("cell1: %d, cell2: %d \n", bdCell->getIndex(), cell->getIndex()); + if(bdCell->getIndex() > (int)gmp_matrix_rows( _HMatrix[dim]) || bdCell->getIndex() < 1 || cell->getIndex() > (int)gmp_matrix_cols( _HMatrix[dim]) || cell->getIndex() < 1){ - Msg::Debug("Warning: Index out of bound! HMatrix: %d. \n", dim); + Msg::Debug("Warning: Index out of bound! HMatrix: %d. \n", + dim); } else{ gmp_matrix_get_elem(elem, bdCell->getIndex(), @@ -121,15 +123,19 @@ void ChainComplex::KerCod(int dim){ if(dim < 0 || dim > 3 || _HMatrix[dim] == NULL) return; - gmp_matrix* HMatrix = copy_gmp_matrix(_HMatrix[dim], 1, 1, - gmp_matrix_rows(_HMatrix[dim]), gmp_matrix_cols(_HMatrix[dim])); + gmp_matrix* HMatrix + = copy_gmp_matrix(_HMatrix[dim], 1, 1, + gmp_matrix_rows(_HMatrix[dim]), + gmp_matrix_cols(_HMatrix[dim])); - gmp_normal_form* normalForm = create_gmp_Hermite_normal_form(HMatrix, NOT_INVERTED, INVERTED); + gmp_normal_form* normalForm + = create_gmp_Hermite_normal_form(HMatrix, NOT_INVERTED, INVERTED); //printMatrix(normalForm->left); //printMatrix(normalForm->canonical); //printMatrix(normalForm->right); - int minRowCol = std::min(gmp_matrix_rows(normalForm->canonical), gmp_matrix_cols(normalForm->canonical)); + int minRowCol = std::min(gmp_matrix_rows(normalForm->canonical), + gmp_matrix_cols(normalForm->canonical)); int rank = 0; mpz_t elem; mpz_init(elem); @@ -142,14 +148,17 @@ void ChainComplex::KerCod(int dim){ } if(rank != (int)gmp_matrix_cols(normalForm->canonical)){ - _kerH[dim] = copy_gmp_matrix(normalForm->right, 1, rank+1, - gmp_matrix_rows(normalForm->right), gmp_matrix_cols(normalForm->right)); + _kerH[dim] + = copy_gmp_matrix(normalForm->right, 1, rank+1, + gmp_matrix_rows(normalForm->right), + gmp_matrix_cols(normalForm->right)); } if(rank > 0){ - _codH[dim] = copy_gmp_matrix(normalForm->canonical, 1, 1, - gmp_matrix_rows(normalForm->canonical), rank); - gmp_matrix_left_mult(normalForm->left, _codH[dim]); + _codH[dim] = + copy_gmp_matrix(normalForm->canonical, 1, 1, + gmp_matrix_rows(normalForm->canonical), rank); + gmp_matrix_left_mult(normalForm->left, _codH[dim]); } mpz_clear(elem); @@ -161,12 +170,18 @@ void ChainComplex::KerCod(int dim){ //j:B_(k+1)->Z_k void ChainComplex::Inclusion(int lowDim, int highDim){ - if(getKerHMatrix(lowDim) == NULL || getCodHMatrix(highDim) == NULL || abs(lowDim-highDim) != 1) return; + if(getKerHMatrix(lowDim) == NULL + || getCodHMatrix(highDim) == NULL + || abs(lowDim-highDim) != 1) return; - gmp_matrix* Zbasis = copy_gmp_matrix(_kerH[lowDim], 1, 1, - gmp_matrix_rows(_kerH[lowDim]), gmp_matrix_cols(_kerH[lowDim])); - gmp_matrix* Bbasis = copy_gmp_matrix(_codH[highDim], 1, 1, - gmp_matrix_rows(_codH[highDim]), gmp_matrix_cols(_codH[highDim])); + gmp_matrix* Zbasis = + copy_gmp_matrix(_kerH[lowDim], 1, 1, + gmp_matrix_rows(_kerH[lowDim]), + gmp_matrix_cols(_kerH[lowDim])); + gmp_matrix* Bbasis + = copy_gmp_matrix(_codH[highDim], 1, 1, + gmp_matrix_rows(_codH[highDim]), + gmp_matrix_cols(_codH[highDim])); int rows = gmp_matrix_rows(Bbasis); @@ -179,7 +194,8 @@ void ChainComplex::Inclusion(int lowDim, int highDim){ // A*inv(V) = U*S - gmp_normal_form* normalForm = create_gmp_Smith_normal_form(Zbasis, INVERTED, INVERTED); + gmp_normal_form* normalForm + = create_gmp_Smith_normal_form(Zbasis, INVERTED, INVERTED); mpz_t elem; mpz_init(elem); @@ -198,7 +214,9 @@ void ChainComplex::Inclusion(int lowDim, int highDim){ - gmp_matrix* LB = copy_gmp_matrix(Bbasis, 1, 1, gmp_matrix_cols(Zbasis), gmp_matrix_cols(Bbasis)); + gmp_matrix* LB = copy_gmp_matrix(Bbasis, 1, 1, + gmp_matrix_cols(Zbasis), + gmp_matrix_cols(Bbasis)); destroy_gmp_matrix(Bbasis); rows = gmp_matrix_rows(LB); @@ -240,12 +258,15 @@ void ChainComplex::Quotient(int dim){ if(dim < 0 || dim > 4 || _JMatrix[dim] == NULL) return; - gmp_matrix* JMatrix = copy_gmp_matrix(_JMatrix[dim], 1, 1, - gmp_matrix_rows(_JMatrix[dim]), gmp_matrix_cols(_JMatrix[dim])); + gmp_matrix* JMatrix = + copy_gmp_matrix(_JMatrix[dim], 1, 1, + gmp_matrix_rows(_JMatrix[dim]), + gmp_matrix_cols(_JMatrix[dim])); int rows = gmp_matrix_rows(JMatrix); int cols = gmp_matrix_cols(JMatrix); - gmp_normal_form* normalForm = create_gmp_Smith_normal_form(JMatrix, NOT_INVERTED, NOT_INVERTED); + gmp_normal_form* normalForm = + create_gmp_Smith_normal_form(JMatrix, NOT_INVERTED, NOT_INVERTED); //printMatrix(normalForm->left); //printMatrix(normalForm->canonical); @@ -266,7 +287,8 @@ void ChainComplex::Quotient(int dim){ int rank = cols - _torsion[dim].size(); if(rows - rank > 0){ - gmp_matrix* Hbasis = copy_gmp_matrix(normalForm->left, 1, rank+1, rows, rows); + gmp_matrix* Hbasis = + copy_gmp_matrix(normalForm->left, 1, rank+1, rows, rows); _QMatrix[dim] = Hbasis; } @@ -302,11 +324,17 @@ void ChainComplex::computeHomology(bool dual){ KerCod(highDim); // 1) no edges, but zero cells - if(lowDim == 0 && !dual && gmp_matrix_cols(getHMatrix(lowDim)) > 0 && getHMatrix(highDim) == NULL) { - setHbasis( setDim, create_gmp_matrix_identity(gmp_matrix_cols(getHMatrix(lowDim))) ); + if(lowDim == 0 && !dual + && gmp_matrix_cols(getHMatrix(lowDim)) > 0 + && getHMatrix(highDim) == NULL) { + setHbasis( setDim, + create_gmp_matrix_identity(gmp_matrix_cols(getHMatrix(lowDim))) ); } - else if(highDim == 0 && dual && gmp_matrix_rows(getHMatrix(highDim)) > 0 && getHMatrix(lowDim) == NULL) { - setHbasis( setDim, create_gmp_matrix_identity(gmp_matrix_rows(getHMatrix(highDim))) ); + else if(highDim == 0 && dual + && gmp_matrix_rows(getHMatrix(highDim)) > 0 + && getHMatrix(lowDim) == NULL) { + setHbasis( setDim, + create_gmp_matrix_identity(gmp_matrix_rows(getHMatrix(highDim))) ); } // 2) this dimension is empty @@ -315,9 +343,10 @@ void ChainComplex::computeHomology(bool dual){ } // 3) No higher dimension cells -> none of the cycles are boundaries else if(getHMatrix(highDim) == NULL){ - setHbasis( setDim, copy_gmp_matrix(getKerHMatrix(lowDim), 1, 1, - gmp_matrix_rows(getKerHMatrix(lowDim)), - gmp_matrix_cols(getKerHMatrix(lowDim))) ); + setHbasis( setDim, + copy_gmp_matrix(getKerHMatrix(lowDim), 1, 1, + gmp_matrix_rows(getKerHMatrix(lowDim)), + gmp_matrix_cols(getKerHMatrix(lowDim))) ); } @@ -329,23 +358,26 @@ void ChainComplex::computeHomology(bool dual){ // 4) No lower dimension cells -> all chains are cycles if(getHMatrix(lowDim) == NULL){ - setKerHMatrix(lowDim, create_gmp_matrix_identity(gmp_matrix_rows(getHMatrix(highDim))) ); + setKerHMatrix(lowDim, + create_gmp_matrix_identity(gmp_matrix_rows(getHMatrix(highDim))) ); } Inclusion(lowDim, highDim); Quotient(lowDim); if(getCodHMatrix(highDim) == NULL){ - setHbasis(setDim, copy_gmp_matrix(getKerHMatrix(lowDim), 1, 1, - gmp_matrix_rows(getKerHMatrix(lowDim)), - gmp_matrix_cols(getKerHMatrix(lowDim))) ); + setHbasis(setDim, + copy_gmp_matrix(getKerHMatrix(lowDim), 1, 1, + gmp_matrix_rows(getKerHMatrix(lowDim)), + gmp_matrix_cols(getKerHMatrix(lowDim))) ); } else if(getJMatrix(lowDim) == NULL || getQMatrix(lowDim) == NULL){ setHbasis(setDim, NULL); } else{ - setHbasis(setDim, copy_gmp_matrix(getKerHMatrix(lowDim), 1, 1, - gmp_matrix_rows(getKerHMatrix(lowDim)), - gmp_matrix_cols(getKerHMatrix(lowDim))) ); + setHbasis(setDim, + copy_gmp_matrix(getKerHMatrix(lowDim), 1, 1, + gmp_matrix_rows(getKerHMatrix(lowDim)), + gmp_matrix_cols(getKerHMatrix(lowDim))) ); gmp_matrix_right_mult(getHbasis(setDim), getQMatrix(lowDim)); } @@ -392,7 +424,8 @@ std::vector<int> ChainComplex::getCoeffVector(int dim, int chainNumber){ std::vector<int> coeffVector; if(dim < 0 || dim > 4) return coeffVector; - if(_Hbasis[dim] == NULL || (int)gmp_matrix_cols(_Hbasis[dim]) < chainNumber) return coeffVector; + if(_Hbasis[dim] == NULL + || (int)gmp_matrix_cols(_Hbasis[dim]) < chainNumber) return coeffVector; int rows = gmp_matrix_rows(_Hbasis[dim]); @@ -416,10 +449,11 @@ std::vector<int> ChainComplex::getCoeffVector(int dim, int chainNumber){ int ChainComplex::getTorsion(int dim, int chainNumber){ if(dim < 0 || dim > 4) return 0; - if(_Hbasis[dim] == NULL || (int)gmp_matrix_cols(_Hbasis[dim]) < chainNumber) return 0; - if(_torsion[dim].empty() || (int)_torsion[dim].size() < chainNumber) return 1; + if(_Hbasis[dim] == NULL + || (int)gmp_matrix_cols(_Hbasis[dim]) < chainNumber) return 0; + if(_torsion[dim].empty() + || (int)_torsion[dim].size() < chainNumber) return 1; else return _torsion[dim].at(chainNumber-1); - } Chain::Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs, @@ -435,7 +469,8 @@ Chain::Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs, if(coeffs.at(i) != 0){ std::list< std::pair<int, Cell*> > subCells; cell->getCells(subCells); - for(std::list< std::pair<int, Cell*> >::iterator it = subCells.begin(); it != subCells.end(); it++){ + for(std::list< std::pair<int, Cell*> >::iterator it = + subCells.begin(); it != subCells.end(); it++){ Cell* subCell = (*it).second; int coeff = (*it).first; _cells.insert( std::make_pair(subCell, coeffs.at(i)*coeff)); @@ -478,7 +513,8 @@ bool Chain::deform(std::map<Cell*, int, Less_Cell> &cellsInChain, } int n = 1; - for(citer cit = cellsNotInChain.begin(); cit != cellsNotInChain.end(); cit++){ + for(citer cit = cellsNotInChain.begin(); cit != cellsNotInChain.end(); + cit++){ Cell* c = (*cit).first; if(n == 2) c->setImmune(true); else c->setImmune(false); @@ -493,17 +529,15 @@ bool Chain::deform(std::map<Cell*, int, Less_Cell> &cellsInChain, bool Chain::deformChain(std::pair<Cell*, int> cell, bool bend){ Cell* c1 = cell.first; - std::map<Cell*, int, Less_Cell> c1Cbd; - c1->getOrgCbd(c1Cbd); - for(citer cit = c1Cbd.begin(); cit != c1Cbd.end(); cit++){ + for(citer cit = c1->firstCoboundary(true); cit != c1->lastCoboundary(true); + cit++){ std::map<Cell*, int, Less_Cell> cellsInChain; std::map<Cell*, int, Less_Cell> cellsNotInChain; Cell* c1CbdCell = (*cit).first; - std::map<Cell*, int, Less_Cell> c1CbdBd; - c1CbdCell->getOrgBd(c1CbdBd); - for(citer cit2 = c1CbdBd.begin(); cit2 != c1CbdBd.end(); cit2++){ + for(citer cit2 = c1CbdCell->firstBoundary(true); + cit2 != c1CbdCell->lastBoundary(true); cit2++){ Cell* c1CbdBdCell = (*cit2).first; int coeff = (*cit2).second; if( (hasCell(c1CbdBdCell) && getCoeff(c1CbdBdCell) != 0) @@ -527,12 +561,9 @@ bool Chain::deformChain(std::pair<Cell*, int> cell, bool bend){ cit2++){ Cell* c = (*cit2).first; if(c->inSubdomain()) next = true; - } - + } if(next) continue; - //printf("dim: %d, in chain: %d, not in chain: %d \n", getDim(), cellsInChain.size(), cellsNotInChain.size()); - if( (getDim() == 1 && cellsInChain.size() == 2 && cellsNotInChain.size() == 1) || (getDim() == 2 && cellsInChain.size() == 3 @@ -581,11 +612,14 @@ void Chain::smoothenChain(){ } deImmuneCells(); - for(citer cit = _cells.begin(); cit != _cells.end(); cit++) deformChain(*cit, false); + for(citer cit = _cells.begin(); cit != _cells.end(); cit++){ + deformChain(*cit, false); + } eraseNullCells(); double t2 = Cpu(); - Msg::Debug("Smoothened a %d-chain from %d cells to %d cells (%g s).\n", getDim(), start, getSize(), t2-t1); + Msg::Debug("Smoothened a %d-chain from %d cells to %d cells (%g s).\n", + getDim(), start, getSize(), t2-t1); return; } @@ -678,7 +712,8 @@ void Chain::createPView(){ _model->setPhysicalName(getName(), getDim(), physicalNum); // only for visualization - PView* chain = new PView(getName(), "ElementData", getGModel(), data, 0, 1); + PView* chain = new PView(getName(), "ElementData", getGModel(), + data, 0, 1); } return; @@ -695,8 +730,12 @@ void Chain::removeCell(Cell* cell) { 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! \n"); + 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! \n"); + } return; } diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp index 7736b9ee8ce0751a21db3a37e36df46aa4372f82..1ca2f803e7dbc39e20ad73ab3a7dfec99270ffe7 100644 --- a/Geo/Homology.cpp +++ b/Geo/Homology.cpp @@ -11,7 +11,8 @@ #include "OS.h" #if defined(HAVE_KBIPACK) -Homology::Homology(GModel* model, std::vector<int> physicalDomain, std::vector<int> physicalSubdomain){ +Homology::Homology(GModel* model, std::vector<int> physicalDomain, + std::vector<int> physicalSubdomain){ _model = model; @@ -69,9 +70,11 @@ Homology::Homology(GModel* model, std::vector<int> physicalDomain, std::vector<i double t2 = Cpu(); Msg::Info("Cell Complex complete (%g s).", t2 - t1); Msg::Info("%d volumes, %d faces, %d edges and %d vertices.", - _cellComplex->getSize(3), _cellComplex->getSize(2), _cellComplex->getSize(1), _cellComplex->getSize(0)); + _cellComplex->getSize(3), _cellComplex->getSize(2), + _cellComplex->getSize(1), _cellComplex->getSize(0)); Msg::StatusBar(2, false, "%d V, %d F, %d E, %d V.", - _cellComplex->getSize(3), _cellComplex->getSize(2), _cellComplex->getSize(1), _cellComplex->getSize(0)); + _cellComplex->getSize(3), _cellComplex->getSize(2), + _cellComplex->getSize(1), _cellComplex->getSize(0)); } Homology::~Homology(){ @@ -90,10 +93,8 @@ void Homology::findGenerators(std::string fileName){ Msg::Info("Reducing the Cell Complex..."); Msg::StatusBar(1, false, "Reducing..."); double t1 = Cpu(); - //_cellComplex->printEuler(); + int omitted = _cellComplex->reduceComplex(); - //_cellComplex->printEuler(); - _cellComplex->combine(3); _cellComplex->reduction(2); @@ -101,16 +102,16 @@ void Homology::findGenerators(std::string fileName){ _cellComplex->reduction(1); _cellComplex->combine(1); - _cellComplex->checkCoherence(); - //_cellComplex->printEuler(); - + _cellComplex->checkCoherence(); double t2 = Cpu(); Msg::Info("Cell Complex reduction complete (%g s).", t2 - t1); Msg::Info("%d volumes, %d faces, %d edges and %d vertices.", - _cellComplex->getSize(3), _cellComplex->getSize(2), _cellComplex->getSize(1), _cellComplex->getSize(0)); + _cellComplex->getSize(3), _cellComplex->getSize(2), + _cellComplex->getSize(1), _cellComplex->getSize(0)); Msg::StatusBar(2, false, "%d V, %d F, %d E, %d N.", - _cellComplex->getSize(3), _cellComplex->getSize(2), _cellComplex->getSize(1), _cellComplex->getSize(0)); + _cellComplex->getSize(3), _cellComplex->getSize(2), + _cellComplex->getSize(1), _cellComplex->getSize(0)); Msg::Info("Computing homology spaces..."); Msg::StatusBar(1, false, "Computing..."); @@ -131,15 +132,22 @@ void Homology::findGenerators(std::string fileName){ convert(i, generator); std::string name = "H" + dimension + getDomainString() + generator; - Chain* chain = new Chain(_cellComplex->getCells(j), chains->getCoeffVector(j,i), _cellComplex, _model, name, chains->getTorsion(j,i)); + Chain* chain = new Chain(_cellComplex->getCells(j), + chains->getCoeffVector(j,i), + _cellComplex, _model, name, + chains->getTorsion(j,i)); t1 = Cpu(); int start = chain->getSize(); chain->smoothenChain(); t2 = Cpu(); - Msg::Info("Smoothened H%d %d from %d cells to %d cells (%g s).", j, i, start, chain->getSize(), t2 - t1); + Msg::Info("Smoothened H%d %d from %d cells to %d cells (%g s).", + j, i, start, chain->getSize(), t2 - t1); if(chain->getSize() != 0) { HRank[j] = HRank[j] + 1; - if(chain->getTorsion() != 1) Msg::Warning("H%d %d has torsion coefficient %d!", j, i, chain->getTorsion()); + if(chain->getTorsion() != 1){ + Msg::Warning("H%d %d has torsion coefficient %d!", + j, i, chain->getTorsion()); + } } _generators[j].push_back(chain); } @@ -149,7 +157,8 @@ void Homology::findGenerators(std::string fileName){ convert(i+1, generator); std::string name = "H" + dimension + getDomainString() + generator; std::vector<int> coeffs (_cellComplex->getOmitted(i).size(),1); - Chain* chain = new Chain(_cellComplex->getOmitted(i), coeffs, _cellComplex, _model, name, 1); + Chain* chain = new Chain(_cellComplex->getOmitted(i), coeffs, + _cellComplex, _model, name, 1); if(chain->getSize() != 0) HRank[j] = HRank[j] + 1; _generators[j].push_back(chain); } @@ -176,31 +185,19 @@ void Homology::findGenerators(std::string fileName){ Msg::Debug("H3 = %d \n", HRank[3]); Msg::StatusBar(1, false, "Homology"); - Msg::StatusBar(2, false, "H0: %d, H1: %d, H2: %d, H3: %d.", HRank[0], HRank[1], HRank[2], HRank[3]); + Msg::StatusBar(2, false, "H0: %d, H1: %d, H2: %d, H3: %d.", + HRank[0], HRank[1], HRank[2], HRank[3]); return; } void Homology::findDualGenerators(std::string fileName){ - //for(int i = 0; i < 4; i++) { printf("Dim %d: \n", i); _cellComplex->printComplex(i); } - Msg::Info("Reducing Cell Complex..."); Msg::StatusBar(1, false, "Reducing..."); double t1 = Cpu(); int omitted = _cellComplex->coreduceComplex(); - /* - _cellComplex->makeDualComplex(); - int omitted = _cellComplex->reduceComplex(_omit); - if(getCombine()){ - _cellComplex->combine(3); - _cellComplex->combine(2); - _cellComplex->combine(1); - } - _cellComplex->makeDualComplex(); - */ - _cellComplex->cocombine(0); _cellComplex->cocombine(1); _cellComplex->cocombine(2); @@ -209,9 +206,11 @@ void Homology::findDualGenerators(std::string fileName){ double t2 = Cpu(); Msg::Info("Cell Complex reduction complete (%g s).", t2 - t1); Msg::Info("%d volumes, %d faces, %d edges and %d vertices.", - _cellComplex->getSize(3), _cellComplex->getSize(2), _cellComplex->getSize(1), _cellComplex->getSize(0)); + _cellComplex->getSize(3), _cellComplex->getSize(2), + _cellComplex->getSize(1), _cellComplex->getSize(0)); Msg::StatusBar(2, false, "%d V, %d F, %d E, %d N.", - _cellComplex->getSize(3), _cellComplex->getSize(2), _cellComplex->getSize(1), _cellComplex->getSize(0)); + _cellComplex->getSize(3), _cellComplex->getSize(2), + _cellComplex->getSize(1), _cellComplex->getSize(0)); Msg::Info("Computing homology spaces..."); @@ -239,11 +238,16 @@ void Homology::findDualGenerators(std::string fileName){ convert(i, generator); std::string name = "H" + dimension + "*" + getDomainString() + generator; - Chain* chain = new Chain(_cellComplex->getCells(j), chains->getCoeffVector(j,i), _cellComplex, _model, name, chains->getTorsion(j,i)); + Chain* chain = new Chain(_cellComplex->getCells(j), + chains->getCoeffVector(j,i), _cellComplex, + _model, name, chains->getTorsion(j,i)); _generators[dim-j].push_back(chain); if(chain->getSize() != 0){ 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()); + if(chain->getTorsion() != 1){ + Msg::Warning("H%d* %d has torsion coefficient %d!", + dim-j, i, chain->getTorsion()); + } } } @@ -252,9 +256,11 @@ void Homology::findDualGenerators(std::string fileName){ for(int i = 0; i < _cellComplex->getNumOmitted(); i++){ std::string generator; convert(i+1, generator); - std::string name = "H" + dimension + "*" + getDomainString() + generator; + std::string name + = "H" + dimension + "*" + getDomainString() + generator; std::vector<int> coeffs (_cellComplex->getOmitted(i).size(),1); - Chain* chain = new Chain(_cellComplex->getOmitted(i), coeffs, _cellComplex, _model, name, 1); + Chain* chain = new Chain(_cellComplex->getOmitted(i), coeffs, + _cellComplex, _model, name, 1); _generators[dim-j].push_back(chain); if(chain->getSize() != 0) HRank[dim-j] = HRank[dim-j] + 1; } @@ -281,7 +287,8 @@ void Homology::findDualGenerators(std::string fileName){ Msg::Debug("H3* = %d \n", HRank[3]); Msg::StatusBar(1, false, "Homology"); - Msg::StatusBar(2, false, "H0*: %d, H1*: %d, H2*: %d, H3*: %d.", HRank[0], HRank[1], HRank[2], HRank[3]); + Msg::StatusBar(2, false, "H0*: %d, H1*: %d, H2*: %d, H3*: %d.", + HRank[0], HRank[1], HRank[2], HRank[3]); return; } @@ -301,8 +308,11 @@ void Homology::computeBettiNumbers(){ Msg::Info("H3 = %d", _cellComplex->getBettiNumber(3)); Msg::StatusBar(1, false, "Homology"); - Msg::StatusBar(2, false, "H0: %d, H1: %d, H2: %d, H3: %d.", _cellComplex->getBettiNumber(0), _cellComplex->getBettiNumber(1), _cellComplex->getBettiNumber(2), _cellComplex->getBettiNumber(3)); - + Msg::StatusBar(2, false, "H0: %d, H1: %d, H2: %d, H3: %d.", + _cellComplex->getBettiNumber(0), + _cellComplex->getBettiNumber(1), + _cellComplex->getBettiNumber(2), + _cellComplex->getBettiNumber(3)); return; } @@ -353,7 +363,8 @@ void Homology::createPViews(){ bool Homology::writeGeneratorsMSH(std::string fileName, bool binary){ if(!_model->writeMSH(fileName, 2.0, binary)) return false; Msg::Info("Wrote homology computation results to %s.", fileName.c_str()); - Msg::Debug("Wrote homology computation results to %s. \n", fileName.c_str()); + Msg::Debug("Wrote homology computation results to %s. \n", + fileName.c_str()); return true; }