diff --git a/Geo/Cell.cpp b/Geo/Cell.cpp
index 9a013fb44894aba7c84a11fe7d2f0436c8d937be..906262739f5b77b7ea8e471c681d58a7f875e963 100755
--- a/Geo/Cell.cpp
+++ b/Geo/Cell.cpp
@@ -30,12 +30,10 @@ bool Less_Cell::operator()(const Cell* c1, const Cell* c2) const
}
-Cell::Cell(MElement* image, bool subdomain, bool boundary) :
+Cell::Cell(MElement* image) :
_combined(false), _index(0), _immune(false), _image(NULL),
- _deleteImage(false), _deleteWithCellComplex(true)
+ _deleteImage(false), _inSubdomain(false)
{
- _onDomainBoundary = boundary;
- _inSubdomain = subdomain;
_dim = image->getDim();
_image = image;
for(int i = 0; i < image->getNumVertices(); i++)
@@ -251,7 +249,6 @@ CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell()
_index = c1->getIndex();
_dim = c1->getDim();
_inSubdomain = c1->inSubdomain();
- _onDomainBoundary = c1->onDomainBoundary();
_combined = true;
_image = NULL;
diff --git a/Geo/Cell.h b/Geo/Cell.h
index fe02f850c076b178404f4e8779e2915f7595df5b..60d0b1742c201d6e3f3edf9bb0b473e35d577d4f 100644
--- a/Geo/Cell.h
+++ b/Geo/Cell.h
@@ -43,9 +43,6 @@ class Cell
// used in relative homology computation
bool _inSubdomain;
- // whether this cell belongs to the boundary of a cell complex
- bool _onDomainBoundary;
-
// whether this cell a combinded cell of elemetary cells
bool _combined;
@@ -70,7 +67,6 @@ class Cell
// Whether to delete the mesh element when done
// (created for homology computation only)
bool _deleteImage;
- bool _deleteWithCellComplex;
// sorted vertices of this cell (used for ordering of the cells)
std::vector<int> _vs;
@@ -81,8 +77,8 @@ class Cell
public:
Cell() : _combined(false), _index(0), _immune(false), _image(NULL),
- _deleteImage(false), _deleteWithCellComplex(true) {}
- Cell(MElement* image, bool subdomain, bool boundary);
+ _deleteImage(false), _inSubdomain(false) {}
+ Cell(MElement* image);
virtual ~Cell();
@@ -97,13 +93,11 @@ class Cell
virtual int getPartition() const { return _image->getPartition(); }
virtual void setImmune(bool immune) { _immune = immune; };
virtual bool getImmune() const { return _immune; };
+ virtual bool inSubdomain() const { return _inSubdomain; }
+ virtual void setInSubdomain(bool subdomain) { _inSubdomain = subdomain; }
virtual void setDeleteImage(bool deleteImage) {
_deleteImage = deleteImage; };
virtual bool getDeleteImage() const { return _deleteImage; };
- virtual void setDeleteWithCellComplex(bool deleteWithCellComplex) {
- _deleteWithCellComplex = deleteWithCellComplex; };
- virtual bool getDeleteWithCellComplex() const {
- return _deleteWithCellComplex; };
// get the number of vertices this cell has
virtual int getNumVertices() const { return _image->getNumVertices(); }
@@ -160,13 +154,6 @@ class Cell
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; }
-
- virtual bool onDomainBoundary() const { return _onDomainBoundary; }
- virtual void setOnDomainBoundary(bool domainboundary) {
- _onDomainBoundary = domainboundary; }
-
// get the number of facets of this cell
virtual int getNumFacets() const;
// get the vertices on a facet of this cell
diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 7d99b7f51b9c5173f60c1bb4e8c0fbd890910fbb..c38226ccad226d9f54462d8711469fae03c676da 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -15,9 +15,6 @@ CellComplex::CellComplex( std::vector<GEntity*> domain,
_domain = domain;
_subdomain = subdomain;
- if(!_domain.empty()) _model = _domain.at(0)->model();
- else _model = NULL;
-
_dim = 0;
_simplicial = true;
@@ -28,96 +25,39 @@ CellComplex::CellComplex( std::vector<GEntity*> domain,
if(i == 0) dim = entity->dim();
if(dim != entity->dim()){
_multidim = true;
- Msg::Warning("Domain is not a manifold.");
+ //printf("Warning: Domain is not a manifold.");
break;
}
}
- // find boundary entities
- //if(!_multidim) find_boundary(domain, subdomain);
-
// insert cells into cell complex
// subdomain need to be inserted first!
- insert_cells(_subdomain, true, true);
- //insert_cells(_boundary, false, true);
- insert_cells(_domain, false, false);
+ if(!insert_cells(_subdomain, true)){ panic_exit(); return; }
+ if(!insert_cells(_domain, false)) { panic_exit(); return; }
for(int i = 0; i < 4; i++){
_ocells[i] = _cells[i];
- _betti[i] = 0;
if(getSize(i) > _dim) _dim = i;
}
}
-void CellComplex::find_boundary(std::vector<GEntity*>& domain,
- std::vector<GEntity*>& subdomain)
-{
- // determine mesh entities on boundary of the domain
- bool duplicate = false;
- for(unsigned int j=0; j < _domain.size(); j++){
- if(_domain.at(j)->dim() == 3){
- std::list<GFace*> faces = _domain.at(j)->faces();
- for(std::list<GFace*>::iterator it = faces.begin(); it != faces.end();
- it++){
- GFace* face = *it;
- duplicate = false;
- for(std::vector<GEntity*>::iterator itb = _boundary.begin();
- itb != _boundary.end(); itb++){
- GEntity* entity = *itb;
- if(face->tag() == entity->tag()){
- _boundary.erase(itb);
- duplicate = true;
- break;
- }
- }
- if(!duplicate) _boundary.push_back(face);
- }
- }
- else if(_domain.at(j)->dim() == 2){
- std::list<GEdge*> edges = _domain.at(j)->edges();
-
- for(std::list<GEdge*>::iterator it = edges.begin(); it != edges.end();
- it++){
- GEdge* edge = *it;
- duplicate = false;
- std::vector<GEntity*>::iterator erase;
- for(std::vector<GEntity*>::iterator itb = _boundary.begin();
- itb != _boundary.end(); itb++){
- GEntity* entity = *itb;
- if(edge->tag() == entity->tag()){
- _boundary.erase(itb);
- duplicate = true;
- break;
- }
- }
- if(!duplicate) _boundary.push_back(edge);
- }
- }
- else if(_domain.at(j)->dim() == 1){
- std::list<GVertex*> vertices = _domain.at(j)->vertices();
- for(std::list<GVertex*>::iterator it = vertices.begin();
- it != vertices.end(); it++){
- GVertex* vertex = *it;
- duplicate = false;
- for(std::vector<GEntity*>::iterator itb = _boundary.begin();
- itb != _boundary.end(); itb++){
- GEntity* entity = *itb;
- if(vertex->tag() == entity->tag()){
- _boundary.erase(itb);
- duplicate = true;
- break;
- }
- }
-
- if(!duplicate) _boundary.push_back(vertex);
- }
+void CellComplex::panic_exit(){
+ for(int i = 0; i < 4; i++){
+ _cells[i].clear();
+ for(citer cit = _ocells[i].begin(); cit != _ocells[i].end(); cit++){
+ Cell* cell = *cit;
+ delete cell;
}
+ _ocells[i].clear();
}
+ for(unsigned int i = 0; i < _newcells.size(); i++) delete _newcells.at(i);
+ _newcells.clear();
+ _domain.clear();
+ _subdomain.clear();
}
-void CellComplex::insert_cells(std::vector<GEntity*>& domain,
- bool subdomain, bool boundary)
+bool CellComplex::insert_cells(std::vector<GEntity*>& domain, bool subdomain)
{
std::vector<MVertex*> vertices;
std::pair<citer, bool> insertInfo;
@@ -146,18 +86,18 @@ void CellComplex::insert_cells(std::vector<GEntity*>& domain,
|| type == MSH_LIN_6 || type == MSH_TET_20 || type == MSH_TET_35
|| type == MSH_TET_56
|| type == MSH_TET_34 || type == MSH_TET_52 ){
- cell = new Cell(element, subdomain, boundary);
+ cell = new Cell(element);
}
else if(type == MSH_QUA_4 || type == MSH_QUA_8 || type == MSH_QUA_9){
- cell = new Cell(element, subdomain, boundary);
+ cell = new Cell(element);
_simplicial = false;
}
else if(type == MSH_HEX_8 || type == MSH_HEX_27 || type == MSH_HEX_20){
- cell = new Cell(element, subdomain, boundary);
+ cell = new Cell(element);
_simplicial = false;
}
else if(type == MSH_PRI_6 || type == MSH_PRI_18 || type == MSH_PRI_15){
- cell = new Cell(element, subdomain, boundary);
+ cell = new Cell(element);
_simplicial = false;
}/* FIXME: no getFaceInfo methods for these MElements
else if(type == MSH_PYR_5 || type == MSH_PYR_14 || type == MSH_PYR_13){
@@ -165,10 +105,11 @@ void CellComplex::insert_cells(std::vector<GEntity*>& domain,
_simplicial = false;
}*/
else {
- Msg::Error("Error: mesh element %d not implemented yet! \n", type);
- continue;
+ //printf("Error: mesh element %d not implemented yet! \n", type);
+ return false;
}
cell->setImmune(false);
+ cell->setInSubdomain(subdomain);
insertInfo = _cells[dim].insert(cell);
if(!insertInfo.second) delete cell;
}
@@ -208,24 +149,27 @@ void CellComplex::insert_cells(std::vector<GEntity*>& domain,
type == MSH_LIN_5 || type == MSH_LIN_6) newtype = MSH_PNT;*/
}
if(newtype == 0){
- Msg::Error("Error: mesh element %d not implemented yet! \n", type); }
+ //printf("Error: mesh element %d not implemented yet! \n", type);
+ return false;
+ }
MElement* element = factory.create(newtype, vertices, 0,
cell->getPartition());
- Cell* newCell = new Cell(element, subdomain, boundary);
+ Cell* newCell = new Cell(element);
newCell->setImmune(cell->getImmune());
+ newCell->setInSubdomain(subdomain);
newCell->setDeleteImage(true);
insertInfo = _cells[dim-1].insert(newCell);
if(!insertInfo.second){
delete newCell;
Cell* oldCell = *(insertInfo.first);
- if(!subdomain && !boundary){
+ if(!subdomain){
int ori = cell->getFacetOri(oldCell);
oldCell->addCoboundaryCell( ori, cell, true);
cell->addBoundaryCell( ori, oldCell, true);
}
}
- else if(!subdomain && !boundary) {
+ else if(!subdomain) {
int ori = cell->getFacetOri(vertices);
cell->addBoundaryCell( ori, newCell, true);
newCell->addCoboundaryCell( ori, cell, true);
@@ -233,6 +177,7 @@ void CellComplex::insert_cells(std::vector<GEntity*>& domain,
}
}
}
+ return true;
}
CellComplex::~CellComplex()
@@ -242,7 +187,7 @@ CellComplex::~CellComplex()
for(citer cit = _ocells[i].begin(); cit != _ocells[i].end(); cit++){
Cell* cell = *cit;
- if(cell->getDeleteWithCellComplex()) delete cell;
+ delete cell;
}
_ocells[i].clear();
@@ -255,7 +200,9 @@ void CellComplex::insertCell(Cell* cell)
{
_newcells.push_back(cell);
std::pair<citer, bool> insertInfo = _cells[cell->getDim()].insert(cell);
- if(!insertInfo.second) Msg::Debug("Warning: Cell not inserted! \n");
+ if(!insertInfo.second){
+ //printf("Warning: Cell not inserted! \n");
+ }
}
void CellComplex::removeCell(Cell* cell, bool other)
@@ -406,10 +353,7 @@ int CellComplex::coreduction(int dim, int omitted)
int CellComplex::reduceComplex(bool omit)
{
- 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));
+ //printf("Cell complex before reduction: %d volumes, %d faces, %d edges and %d vertices.\n", getSize(3), getSize(2), getSize(1), getSize(0));
int count = 0;
for(int i = 3; i > 0; i--) count = count + reduction(i);
@@ -418,7 +362,7 @@ int CellComplex::reduceComplex(bool omit)
int omitted = 0;
_store.clear();
- //removeSubdomain();
+ removeSubdomain();
while (getSize(getDim()) != 0){
@@ -434,9 +378,7 @@ int CellComplex::reduceComplex(bool omit)
}
}
- 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);
+ //printf("Cell complex after reduction: %d volumes, %d faces, %d edges and %d vertices.\n", getSize(3), getSize(2), getSize(1), getSize(0));
return 0;
}
@@ -454,9 +396,7 @@ void CellComplex::removeSubdomain()
int CellComplex::coreduceComplex()
{
-
- Msg::Debug("Cell complex before coreduction: %d volumes, %d faces, %d edges and %d vertices.\n",
- getSize(3), getSize(2), getSize(1), getSize(0));
+ //printf("Cell complex before coreduction: %d volumes, %d faces, %d edges and %d vertices.\n", getSize(3), getSize(2), getSize(1), getSize(0));
int count = 0;
@@ -487,52 +427,14 @@ int CellComplex::coreduceComplex()
coreduction(cell, omitted);
}
- Msg::Debug("Cell complex after coreduction: %d volumes, %d faces, %d edges and %d vertices.\n",
- getSize(3), getSize(2), getSize(1), getSize(0));
+ //printf("Cell complex after coreduction: %d volumes, %d faces, %d edges and %d vertices.\n", getSize(3), getSize(2), getSize(1), getSize(0));
return omitted;
}
-void CellComplex::computeBettiNumbers()
-{
- for(int i = 0; i < 4; i++){
- _betti[i] = 0;
- Msg::Debug("Betti number computation process: step %d of 4 \n", i+1);
- while (getSize(i) != 0){
- /*removeSubdomain();
- citer cit = firstCell(i);
- while(cit != lastCell(i)){
- if(coreduction(*cit)) cit = firstCell(i);
- else cit++;
- }
- if(getSize(i) != 0) {
- _betti[i] = _betti[i] + 1;
- cit = firstCell(i);
- removeCell(*cit, false);
- }*/
- citer cit = firstCell(i);
- Cell* cell = *cit;
- while(!cell->inSubdomain() && cit != lastCell(i)){
- cell = *cit;
- cit++;
- }
- if(!cell->inSubdomain()) _betti[i] = _betti[i] + 1;
- removeCell(cell, false);
- coreduction(cell);
- }
- }
- Msg::Debug("Cell complex Betti numbers: \nH0 = %d \nH1 = %d \nH2 = %d \nH3 = %d \n",
- getBettiNumber(0), getBettiNumber(1),
- getBettiNumber(2), getBettiNumber(3));
- return;
-}
-
-
int CellComplex::cocombine(int dim)
{
- double t1 = Cpu();
- Msg::Debug("Cell complex before cocombining: %d volumes, %d faces, %d edges and %d vertices.\n",
- getSize(3), getSize(2), getSize(1), getSize(0));
+ //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;
@@ -584,18 +486,14 @@ int CellComplex::cocombine(int dim)
}
}
- double t2 = Cpu();
- Msg::Debug("Cell complex after cocombining: %d volumes, %d faces, %d edges and %d vertices (%g s).\n",
- getSize(3), getSize(2), getSize(1), getSize(0), t2-t1);
+ //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;
}
int CellComplex::combine(int dim)
{
- double t1 = Cpu();
- Msg::Debug("Cell complex before combining: %d volumes, %d faces, %d edges and %d vertices.\n",
- getSize(3), getSize(2), getSize(1), getSize(0));
+ //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;
@@ -649,9 +547,7 @@ int CellComplex::combine(int dim)
}
}
- double t2 = Cpu();
- Msg::Debug("Cell complex after combining: %d volumes, %d faces, %d edges and %d vertices (%g s).\n",
- getSize(3), getSize(2), getSize(1), getSize(0), t2 - t1);
+ //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;
@@ -664,7 +560,7 @@ void CellComplex::printComplex(int dim)
cell->printCell();
cell->printBoundary();
cell->printCoboundary();
- printf("--- \n");
+ //printf("--- \n");
}
}
@@ -682,14 +578,14 @@ bool CellComplex::checkCoherence()
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("Warning! Boundary cell not in cell complex! Boundary removed. \n");
//cell->printCell();
//bdCell->printCell();
cell->removeBoundaryCell(bdCell);
coherent = false;
}
if(!bdCell->hasCoboundary(cell)){
- Msg::Debug("Warning! Incoherent boundary/coboundary pair! Fixed. \n");
+ //printf("Warning! Incoherent boundary/coboundary pair! Fixed. \n");
/*cell->printCell();
cell->printBoundary();
bdCell->printCell();
@@ -707,14 +603,14 @@ bool CellComplex::checkCoherence()
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("Warning! Coboundary cell not in cell complex! Coboundary removed. \n");
cell->printCell();
cbdCell->printCell();
cell->removeCoboundaryCell(cbdCell);
coherent = false;
}
if(!cbdCell->hasBoundary(cell)){
- Msg::Debug("Warning! Incoherent coboundary/boundary pair! Fixed. \n");
+ //printf("Warning! Incoherent coboundary/boundary pair! Fixed. \n");
cell->printCell();
cell->printCoboundary();
cbdCell->printCell();
@@ -730,19 +626,6 @@ bool CellComplex::checkCoherence()
return coherent;
}
-void CellComplex::restoreComplex()
-{
- for(int i = 0; i < 4; i++){
- _betti[i] = 0;
- _cells[i] = _ocells[i];
- for(citer cit = firstCell(i); cit != lastCell(i); cit++){
- Cell* cell = *cit;
- cell->restoreCell();
- }
- }
- _store.clear();
-}
-
bool CellComplex::hasCell(Cell* cell, bool org)
{
if(!org){
@@ -768,95 +651,4 @@ void CellComplex::getCells(std::set<Cell*, Less_Cell>& cells,
}
}
}
-
-bool CellComplex::swapSubdomain()
-{
- if(_multidim) return false;
- 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);
- }
- }
- return true;
-}
-
-void CellComplex::makeDualComplex()
-{
- std::set<Cell*, Less_Cell> temp = _cells[0];
- _cells[0] = _cells[3];
- _cells[3] = temp;
- temp = _cells[1];
- _cells[1] = _cells[2];
- _cells[2] = temp;
-
- for(int i = 0; i < 4; i++){
- for(citer cit = firstCell(i); cit != lastCell(i); cit++){
- Cell* cell = *cit;
- cell->makeDualCell();
- }
- }
-}
-
-bool CellComplex::writeBettiNumbers(std::string fileName)
-{
- if(fileName.empty()) return false;
- FILE *fp = fopen(fileName.c_str(), "w");
- if(!fp){
- Msg::Error("Unable to open file '%s'", fileName.c_str());
- return false;
- }
-
- fprintf(fp, "$BettiNumbers\n");
- fprintf(fp, "H0 %d \n", getBettiNumber(0));
- fprintf(fp, "H1 %d \n", getBettiNumber(1));
- fprintf(fp, "H2 %d \n", getBettiNumber(2));
- fprintf(fp, "H3 %d \n", getBettiNumber(3));
- fprintf(fp, "$EndBettiNumbers\n");
-
- fclose(fp);
- Msg::Info("Wrote Betti numbers to %s.", fileName.c_str());
- return true;
-}
-
-void CellComplex::storeCells(int dim)
-{
- std::vector<MElement*> elements;
-
- for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
- Cell* cell = *cit;
-
- std::list< std::pair<int, Cell*> > cells;
- cell->getCells(cells);
- for(std::list< std::pair<int, Cell*> >::iterator it = cells.begin();
- it != cells.end(); it++){
- Cell* subCell = it->second;
-
- MElement* e = subCell->getImageMElement();
- subCell->setDeleteImage(false);
- elements.push_back(e);
- }
- }
-
- int max[4];
- for(int i = 0; i < 4; i++) max[i] = _model->getMaxElementaryNumber(i);
- int entityNum = *std::max_element(max,max+4) + 1;
- for(int i = 0; i < 4; i++) max[i] = _model->getMaxPhysicalNumber(i);
- int physicalNum = *std::max_element(max,max+4) + 1;
-
- 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]="Cell Complex";
- physicalMap[entityNum] = physicalInfo;
-
- _model->storeChain(dim, entityMap, physicalMap);
- _model->setPhysicalName("Cell Complex", dim, physicalNum);
-
-}
-
#endif
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index 8a77c8cc38622217ee1b35139cbd1ce0e0bf8a18..05586ac5ea7d909acf28dc865b08fa8a632ba31f 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -20,7 +20,6 @@
#include "Cell.h"
#include "MElement.h"
#include "ChainComplex.h"
-#include "OS.h"
#include "GModel.h"
class Cell;
@@ -29,8 +28,6 @@ class Cell;
class CellComplex
{
private:
-
- GModel* _model;
// the domain in the model which this cell complex covers
std::vector<GEntity*> _domain;
@@ -54,9 +51,6 @@ class CellComplex
// new cells created during reductions
std::vector<Cell*> _newcells;
- // Betti numbers of this cell complex (ranks of homology groups)
- int _betti[4];
-
int _dim;
// is the cell complex simplicial
@@ -72,11 +66,9 @@ class CellComplex
void removeCellQset(Cell* cell, std::set<Cell*, Less_Cell>& Qset);
// for constructor
- void insert_cells(std::vector<GEntity*>& domain,
- bool subdomain, bool boundary);
- void find_boundary(std::vector<GEntity*>& domain,
- std::vector<GEntity*>& subdomain);
-
+ bool insert_cells(std::vector<GEntity*>& domain, bool subdomain);
+ void panic_exit();
+
// insert/remove a cell from this cell complex
void removeCell(Cell* cell, bool other=true);
void insertCell(Cell* cell);
@@ -93,18 +85,15 @@ class CellComplex
std::vector<GEntity*> getDomain() const { return _domain; }
std::vector<GEntity*> getSubdomain() const { return _subdomain; }
- // restore this cell complex to its original state
- void restoreComplex();
-
// get the number of certain dimensional cells
- int getSize(int dim){ return _cells[dim].size(); }
- int getOrgSize(int dim){ return _ocells[dim].size(); }
+ int getSize(int dim, bool org=false){
+ if(!org) return _cells[dim].size();
+ else return _ocells[dim].size(); }
int getDim() {return _dim; }
bool simplicial() { return _simplicial; }
- std::set<Cell*, Less_Cell> getCells(int dim){ return _cells[dim]; }
void getCells(std::set<Cell*, Less_Cell>& cells, int dim, int domain=0);
std::set<Cell*, Less_Cell> getOrgCells(int dim){ return _ocells[dim]; }
@@ -123,11 +112,8 @@ 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()))
- && ((!c1->getImmune() && !c2->getImmune())
- || (c1->getImmune() && c2->getImmune())) ); }
-
+ return ( (!c1->inSubdomain() && !c2->inSubdomain())
+ || (c1->inSubdomain() && c2->inSubdomain() )); }
// (co)reduction of this cell complex
// removes (co)reduction pairs of cell of dimension dim and dim-1
int reduction(int dim, int omitted=0);
@@ -139,7 +125,6 @@ class CellComplex
// remove cells in subdomain from this cell complex
void removeSubdomain();
- void deImmuneCells(bool subdomain=false);
// print the vertices of cells of certain dimension
void printComplex(int dim);
@@ -148,39 +133,20 @@ class CellComplex
int combine(int dim);
int cocombine(int dim);
- // Compute betti numbers of this cell complex
- void computeBettiNumbers();
- int getBettiNumber(int i) {
- if(i > -1 && i < 4) return _betti[i]; else return 0; }
- bool writeBettiNumbers(std::string fileName);
-
// check whether all boundary cells of a cell has this cell
// as coboundary cell and vice versa
// also check whether all (co)boundary cells of a cell
// belong to this cell complex
bool checkCoherence();
- // make cells on the domain boundary that are not in the subdomain
- // subdomain cells and vice versa
- bool swapSubdomain();
-
int eulerCharacteristic(){
return getSize(0) - getSize(1) + getSize(2) - getSize(3);}
void printEuler(){
printf("Euler characteristic: %d. \n", eulerCharacteristic()); }
- int getNumOmitted() { return _store.size(); };
+ int getNumOmitted() { return _store.size(); }
std::set<Cell*, Less_Cell> getOmitted(int i) { return _store.at(i); }
- // change roles of boundaries and coboundaries of the cells in this
- // cell complex
- // equivalent to transposing boundary operator matrices
- void makeDualComplex();
-
- // store cells of dimension dim to the GModel as a
- // physical group of MElements
- // for debugging purposes
- void storeCells(int dim);
};
#endif
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index 1dd04bb090ac6cf978e76af155d740b37c86fc5e..21f722e307ad6691ce7453e60b80c9646052968f 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -28,32 +28,34 @@ ChainComplex::ChainComplex(CellComplex* cellComplex, int domain)
if(dim > 0) rows = cellComplex->getSize(dim-1);
int index = 1;
- // ignore subdomain cells
+ // ignore cells depending on domain
for(CellComplex::citer cit = cellComplex->firstCell(dim);
cit != cellComplex->lastCell(dim); cit++){
Cell* cell = *cit;
- cell->setIndex(index);
- index++;
- if((domain == 0 && cell->inSubdomain())
- || (domain == 2 && !cell->inSubdomain()) ){
- index--;
- cols--;
+ cell->setIndex(0);
+ cols--;
+ if((domain == 0 && !cell->inSubdomain())
+ || (domain == 2 && cell->inSubdomain()) ){
+ cols++;
+ cell->setIndex(index);
+ index++;
+ _cellIndices[dim][cell->getIndex()] = cell;
}
- else _cellIndices[dim].insert( std::make_pair(index, cell));
}
index = 1;
if(dim > 0){
for(CellComplex::citer cit = cellComplex->firstCell(dim-1);
cit != cellComplex->lastCell(dim-1); cit++){
Cell* cell = *cit;
- cell->setIndex(index);
- index++;
- if( (domain == 0 && cell->inSubdomain())
- || (domain == 2 && !cell->inSubdomain()) ){
- index--;
- rows--;
- }
- else _cellIndices[dim-1].insert( std::make_pair(index, cell));
+ cell->setIndex(0);
+ rows--;
+ if((domain == 0 && !cell->inSubdomain())
+ || (domain == 2 && cell->inSubdomain()) ){
+ rows++;
+ cell->setIndex(index);
+ index++;
+ _cellIndices[dim-1][cell->getIndex()] = cell;
+ }
}
}
@@ -87,8 +89,7 @@ ChainComplex::ChainComplex(CellComplex* cellComplex, int domain)
|| 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);
+ //printf("Warning: Index out of bound! HMatrix: %d. \n", dim);
}
else{
gmp_matrix_get_elem(elem, bdCell->getIndex(),
@@ -96,8 +97,7 @@ ChainComplex::ChainComplex(CellComplex* cellComplex, int domain)
old_elem = mpz_get_si(elem);
mpz_set_si(elem, old_elem + (*it).second);
if( abs((old_elem + (*it).second)) > 1){
- Msg::Debug("Incidence index: %d! HMatrix: %d. \n",
- (old_elem + (*it).second), dim);
+ //printf("Incidence index: %d, in HMatrix: %d. \n", (old_elem + (*it).second), dim);
}
gmp_matrix_set_elem(elem, bdCell->getIndex(),
cell->getIndex(), _HMatrix[dim]);
@@ -312,7 +312,7 @@ void ChainComplex::computeHomology(bool dual)
//KerCod(highDim);
}
- Msg::Debug("Homology computation process: step %d of 4 \n", i+1);
+ //printf("Homology computation process: step %d of 4 \n", i+1);
KerCod(highDim);
@@ -433,7 +433,8 @@ std::vector<int> ChainComplex::getCoeffVector(int dim, int chainNumber)
return coeffVector;
}
-void ChainComplex::getChain(std::map<Cell*, int, Less_Cell>& chain, int dim, int chainNumber)
+void ChainComplex::getChain(std::map<Cell*, int, Less_Cell>& chain,
+ int dim, int chainNumber)
{
chain.clear();
if(dim < 0 || dim > 4) return;
@@ -452,9 +453,20 @@ void ChainComplex::getChain(std::map<Cell*, int, Less_Cell>& chain, int dim, int
elemli = mpz_get_si(elem);
elemi = elemli;
Cell* cell = _cellIndices[dim][i];
- chain[cell] = elemi;
- }
+ if(cell == NULL){
+ //printf("Warning: Missing cell in %d-chain %d. \n", dim, chainNumber);
+ continue;
+ }
+ 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++){
+ Cell* subCell = (*it).second;
+ int coeff = (*it).first;
+ chain[subCell] = coeff*elemi;
+ }
+ }
mpz_clear(elem);
}
diff --git a/Geo/ChainComplex.h b/Geo/ChainComplex.h
index fce3e0a9502cef9f506bad4509547d6552dbb518..4c5ebece479808e1c142a8aa6a1d0ad65ad9289f 100644
--- a/Geo/ChainComplex.h
+++ b/Geo/ChainComplex.h
@@ -18,7 +18,6 @@
#include <set>
#include <queue>
#include "CellComplex.h"
-#include "OS.h"
#if defined(HAVE_GMP)
#include "gmp.h"
@@ -73,19 +72,10 @@ class ChainComplex{
public:
- ChainComplex(CellComplex* cellComplex, int domain=0);
-
- ChainComplex(){
- for(int i = 0; i < 5; i++){
- _HMatrix[i] = create_gmp_matrix_zero(1,1);
- _kerH[i] = NULL;
- _codH[i] = NULL;
- _JMatrix[i] = NULL;
- _QMatrix[i] = NULL;
- _Hbasis[i] = NULL;
- }
- _dim = 0;
- }
+ // domain = 0 : relative chain space
+ // domain = 1 : absolute chain space of all cells in cellComplex
+ // domain = 2 : absolute chain space of cells in subdomain
+ ChainComplex(CellComplex* cellComplex, int domain=0);
~ChainComplex(){
for(int i = 0; i < 5; i++){
destroy_gmp_matrix(_HMatrix[i]);
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index 85f3a9d58c69800e9cefe1e544c696b6252fd77b..929e3e838268b5087b7d4b2a2058918b97ed0f88 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -153,10 +153,9 @@ void Homology::findGenerators()
std::string name = "H" + dimension + domainString + generator;
std::set<Cell*, Less_Cell> cells;
cellComplex->getCells(cells, j);
- Chain* chain = new Chain(cells,
- chains->getCoeffVector(j,i),
- cellComplex, _model, name,
- chains->getTorsion(j,i));
+ Chain* chain = new Chain(cells,
+ chains->getCoeffVector(j,i), cellComplex,
+ _model, name, chains->getTorsion(j,i));
t1 = Cpu();
int start = chain->getSize();
chain->smoothenChain();
@@ -170,7 +169,8 @@ void Homology::findGenerators()
j, i, chain->getTorsion());
}
}
- _generators[chain->createPGroup()] = chain;
+ _generators.push_back(chain->createPGroup());
+ delete chain;
}
if(j == cellComplex->getDim() && cellComplex->getNumOmitted() > 0){
for(int i = 0; i < cellComplex->getNumOmitted(); i++){
@@ -181,7 +181,8 @@ void Homology::findGenerators()
Chain* chain = new Chain(cellComplex->getOmitted(i), coeffs,
cellComplex, _model, name, 1);
if(chain->getSize() != 0) HRank[j] = HRank[j] + 1;
- _generators[chain->createPGroup()] = chain;
+ _generators.push_back(chain->createPGroup());
+ delete chain;
}
}
}
@@ -211,8 +212,16 @@ void Homology::findDualGenerators()
Msg::Info("Reducing Cell Complex...");
Msg::StatusBar(1, false, "Reducing...");
double t1 = Cpu();
+ printf("Cell Complex: \n %d volumes, %d faces, %d edges and %d vertices. \n",
+ cellComplex->getSize(3), cellComplex->getSize(2),
+ cellComplex->getSize(1), cellComplex->getSize(0));
+
int omitted = cellComplex->coreduceComplex();
+ printf(" %d volumes, %d faces, %d edges and %d vertices. \n",
+ cellComplex->getSize(3), cellComplex->getSize(2),
+ cellComplex->getSize(1), cellComplex->getSize(0));
+
cellComplex->cocombine(0);
cellComplex->coreduction(1);
cellComplex->cocombine(1);
@@ -220,6 +229,12 @@ void Homology::findDualGenerators()
cellComplex->cocombine(2);
cellComplex->coreduction(3);
cellComplex->checkCoherence();
+
+ printf(" %d volumes, %d faces, %d edges and %d vertices. \n",
+ cellComplex->getSize(3), cellComplex->getSize(2),
+ cellComplex->getSize(1), cellComplex->getSize(0));
+
+
double t2 = Cpu();
Msg::Info("Cell Complex reduction complete (%g s).", t2 - t1);
Msg::Info("%d volumes, %d faces, %d edges and %d vertices.",
@@ -258,7 +273,8 @@ void Homology::findDualGenerators()
Chain* chain = new Chain(cells,
chains->getCoeffVector(j,i), cellComplex,
_model, name, chains->getTorsion(j,i));
- _generators[chain->createPGroup()] = chain;
+ _generators.push_back(chain->createPGroup());
+ delete chain;
if(chain->getSize() != 0){
HRank[dim-j] = HRank[dim-j] + 1;
if(chain->getTorsion() != 1){
@@ -279,8 +295,9 @@ void Homology::findDualGenerators()
std::vector<int> coeffs (cellComplex->getOmitted(i).size(),1);
Chain* chain = new Chain(cellComplex->getOmitted(i), coeffs,
cellComplex, _model, name, 1);
- _generators[chain->createPGroup()] = chain;
if(chain->getSize() != 0) HRank[dim-j] = HRank[dim-j] + 1;
+ _generators.push_back(chain->createPGroup());
+ delete chain;
}
}
}
@@ -302,34 +319,6 @@ void Homology::findDualGenerators()
HRank[0], HRank[1], HRank[2], HRank[3]);
}
-void Homology::computeBettiNumbers()
-{
- CellComplex* cellComplex = createCellComplex(_domainEntities,
- _subdomainEntities);
-
- Msg::Info("Running coreduction...");
- Msg::StatusBar(1, false, "Computing...");
- double t1 = Cpu();
- cellComplex->computeBettiNumbers();
- double t2 = Cpu();
- Msg::Info("Betti number computation complete (%g s).", t2- t1);
-
- Msg::Info("H0 = %d", cellComplex->getBettiNumber(0));
- Msg::Info("H1 = %d", cellComplex->getBettiNumber(1));
- Msg::Info("H2 = %d", cellComplex->getBettiNumber(2));
- 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));
-
- if(_fileName != "") cellComplex->writeBettiNumbers(_fileName);
- delete cellComplex;
-}
-
void Homology::findHomSequence(){
CellComplex* cellComplex = new CellComplex(_domainEntities,
@@ -354,8 +343,6 @@ void Homology::findHomSequence(){
cellComplex->reduction(1);
cellComplex->combine(1);
- cellComplex->storeCells(1);
-
printf(" %d volumes, %d faces, %d edges and %d vertices. \n",
cellComplex->getSize(3), cellComplex->getSize(2),
cellComplex->getSize(1), cellComplex->getSize(0));
@@ -416,7 +403,8 @@ void Homology::findHomSequence(){
j, i, chain->getTorsion());
}
}
- _generators[chain->createPGroup()] = chain;
+ _generators.push_back(chain->createPGroup());
+ delete chain;
}
}
@@ -444,12 +432,6 @@ void Homology::findHomSequence(){
delete cellComplex;
}
-/*void Homology::restoreHomology()
-{
- cellComplexy->restoreComplex();
- for(int i = 0; i < 4; i++) _generators[i].clear();
- }*/
-
std::string Homology::getDomainString(const std::vector<int>& domain,
const std::vector<int>& subdomain)
{
@@ -508,7 +490,6 @@ Chain::Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs,
Cell* subCell = (*it).second;
int coeff = (*it).first;
_cells.insert( std::make_pair(subCell, coeffs.at(i)*coeff));
- subCell->setDeleteWithCellComplex(false);
}
}
i++;
@@ -522,17 +503,55 @@ Chain::Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs,
}
+void Homology::storeCells(CellComplex* cellComplex, int dim)
+{
+ std::vector<MElement*> elements;
+
+ for(CellComplex::citer cit = cellComplex->firstCell(dim);
+ cit != cellComplex->lastCell(dim); cit++){
+ Cell* cell = *cit;
+
+ std::list< std::pair<int, Cell*> > cells;
+ cell->getCells(cells);
+ for(std::list< std::pair<int, Cell*> >::iterator it = cells.begin();
+ it != cells.end(); it++){
+ Cell* subCell = it->second;
+
+ MElement* e = subCell->getImageMElement();
+ subCell->setDeleteImage(false);
+ elements.push_back(e);
+ }
+ }
+
+ int max[4];
+ for(int i = 0; i < 4; i++) max[i] = _model->getMaxElementaryNumber(i);
+ int entityNum = *std::max_element(max,max+4) + 1;
+ for(int i = 0; i < 4; i++) max[i] = _model->getMaxPhysicalNumber(i);
+ int physicalNum = *std::max_element(max,max+4) + 1;
+
+ 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]="Cell Complex";
+ physicalMap[entityNum] = physicalInfo;
+
+ _model->storeChain(dim, entityMap, physicalMap);
+ _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();
}
bool Chain::deform(std::map<Cell*, int, Less_Cell>& cellsInChain,
@@ -788,7 +807,6 @@ 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;
- cell->setDeleteWithCellComplex(false);
}
else if (!insert.second && (*insert.first).second != 0){
Msg::Debug("Error: invalid chain smoothening add! \n");
@@ -825,10 +843,7 @@ void Chain::eraseNullCells()
if((*cit).second == 0) toRemove.push_back((*cit).first);
}
}
- for(unsigned int i = 0; i < toRemove.size(); i++){
- _cells.erase(toRemove[i]);
- toRemove[i]->setDeleteWithCellComplex(true);
- }
+ for(unsigned int i = 0; i < toRemove.size(); i++) _cells.erase(toRemove[i]);
return;
}
diff --git a/Geo/Homology.h b/Geo/Homology.h
index 49f80aa0aa993769c617f815dc3df86112b76121..88a630ea2ccf56b1a961f765b5b92d2076a5cf5e 100644
--- a/Geo/Homology.h
+++ b/Geo/Homology.h
@@ -41,11 +41,9 @@ class Homology
std::vector<GEntity*> _domainEntities;
std::vector<GEntity*> _subdomainEntities;
- // generator chains
- //std::vector<Chain*> _generators[4];
- //std::vector<int> _pGroups;
-
- std::map<int, Chain*> _generators;
+ // generator chains, and their physical group number
+ //std::map<int, Chain*> _generators;
+ std::vector<int> _generators;
std::string _fileName;
@@ -64,14 +62,10 @@ class Homology
// or just compute the ranks of homology spaces
void findGenerators();
void findDualGenerators();
- void computeBettiNumbers();
+ void computeRanks() {}
void findHomSequence();
- //bool swapSubdomain() { return _cellComplex->swapSubdomain(); }
-
- // Restore the cell complex to its original state before cell reductions
- //void restoreHomology();
// Create a string describing the generator
std::string getDomainString(const std::vector<int>& domain,
@@ -79,6 +73,10 @@ class Homology
// write the generators to a file
bool writeGeneratorsMSH(bool binary=false);
+ // store dim-dimensional cells of cellComplex as a pshysical group
+ // in _model, for debugging
+ void storeCells(CellComplex* cellComplex, int dim);
+
};
// A class representing a chain.
diff --git a/Parser/Gmsh.tab.cpp b/Parser/Gmsh.tab.cpp
index efd62cfb978bd138affef7191bea0924ac99080e..7268617aa70159badaf6be2c5fdca08475b9cfaa 100644
--- a/Parser/Gmsh.tab.cpp
+++ b/Parser/Gmsh.tab.cpp
@@ -7937,7 +7937,7 @@ yyreduce:
#if defined(HAVE_KBIPACK)
Homology* homology = new Homology(GModel::current(), domain, subdomain);
homology->setFileName(fileName);
- homology->computeBettiNumbers();
+ homology->computeRanks();
delete homology;
#else
yymsg(0, "Gmsh needs to be configured with option Kbipack to use homology computation.");
diff --git a/Parser/Gmsh.y b/Parser/Gmsh.y
index 955f3f4f6a07d378be4ce46f197f8f847cb8787e..edc45ae33c21ef529939738ee2b1e5a021cd753b 100644
--- a/Parser/Gmsh.y
+++ b/Parser/Gmsh.y
@@ -3333,7 +3333,7 @@ Homology :
#if defined(HAVE_KBIPACK)
Homology* homology = new Homology(GModel::current(), domain, subdomain);
homology->setFileName(fileName);
- homology->computeBettiNumbers();
+ homology->computeRanks();
delete homology;
#else
yymsg(0, "Gmsh needs to be configured with option Kbipack to use homology computation.");
diff --git a/Plugin/HomologyComputation.cpp b/Plugin/HomologyComputation.cpp
index 7027f294c8174101095937e323c2eb4d0449b8f7..33739fc361737f6d5767478a4c93252d8dac8a81 100644
--- a/Plugin/HomologyComputation.cpp
+++ b/Plugin/HomologyComputation.cpp
@@ -108,7 +108,7 @@ PView *GMSH_HomologyComputationPlugin::execute(PView *v)
if(rank != 0){
Homology* homology = new Homology(m, domain, subdomain);
homology->setFileName(fileName);
- homology->computeBettiNumbers();
+ homology->computeRanks();
delete homology;
}