Skip to content
Snippets Groups Projects
Commit 61a9dafb authored by Matti Pellika's avatar Matti Pellika
Browse files 

Code cleanup, increased verbosity.
parent ef606b92
No related branches found
No related tags found
No related merge requests found
Expand all Show whitespace changes
Inline Side-by-side
Geo/CMakeLists.txt
+ 1
− 1
View file @ 61a9dafb
...@@ -31,7 +31,7 @@ set(SRC ...@@ -31,7 +31,7 @@ set(SRC
MLine.cpp MTriangle.cpp MQuadrangle.cpp MTetrahedron.cpp MLine.cpp MTriangle.cpp MQuadrangle.cpp MTetrahedron.cpp
MHexahedron.cpp MPrism.cpp MPyramid.cpp MElementCut.cpp MHexahedron.cpp MPrism.cpp MPyramid.cpp MElementCut.cpp
MZone.cpp MZoneBoundary.cpp MZone.cpp MZoneBoundary.cpp
CellComplex.cpp ChainComplex.cpp Homology.cpp Cell.cpp CellComplex.cpp ChainComplex.cpp Homology.cpp
) )
file(GLOB HDR RELATIVE ${CMAKE_SOURCE_DIR}/Geo *.h) file(GLOB HDR RELATIVE ${CMAKE_SOURCE_DIR}/Geo *.h)
......
Geo/CellComplex.cpp
+ 109
− 143
View file @ 61a9dafb
...@@ -18,7 +18,7 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su ...@@ -18,7 +18,7 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su
_dim = 0; _dim = 0;
_simplicial = true; _simplicial = true;
// find boundary elements // find boundary entities
find_boundary(domain, subdomain); find_boundary(domain, subdomain);
// insert cells into cell complex // insert cells into cell complex
...@@ -182,7 +182,7 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){ ...@@ -182,7 +182,7 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){
cell = new CPyramid(vertices, tag, subdomain, boundary); cell = new CPyramid(vertices, tag, subdomain, boundary);
_simplicial = false; _simplicial = false;
}*/ }*/
else printf("Error: mesh element %d not implemented yet! \n", type); else Msg::Error("Error: mesh element %d not implemented yet! \n", type);
tag++; tag++;
//if(domain.at(j)->getMeshElement(i)->getNum() == 8196) cell->setImmune(true); //if(domain.at(j)->getMeshElement(i)->getNum() == 8196) cell->setImmune(true);
//else //else
...@@ -203,7 +203,7 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){ ...@@ -203,7 +203,7 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){
if(dim == 3){ if(dim == 3){
if(vertices.size() == 3) newCell = new TwoSimplex(vertices, tag, subdomain, boundary); if(vertices.size() == 3) newCell = new TwoSimplex(vertices, tag, subdomain, boundary);
else if(vertices.size() == 4) newCell = new CQuadrangle(vertices, tag, subdomain, boundary); else if(vertices.size() == 4) newCell = new CQuadrangle(vertices, tag, subdomain, boundary);
else printf("Error: invalid face! \n"); else Msg::Debug("Error: invalid face! \n");
} }
else if(dim == 2) newCell = new OneSimplex(vertices, tag, subdomain, boundary); else if(dim == 2) newCell = new OneSimplex(vertices, tag, subdomain, boundary);
else if(dim == 1) newCell = new ZeroSimplex(vertices, tag, subdomain, boundary); else if(dim == 1) newCell = new ZeroSimplex(vertices, tag, subdomain, boundary);
...@@ -236,98 +236,54 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){ ...@@ -236,98 +236,54 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){
} }
void CellComplex::insertCell(Cell* cell){ CellComplex::~CellComplex(){
_ecells.insert(cell); for(int i = 0; i < 4; i++){
} _cells[i].clear();
/* for(citer cit = _ocells[i].begin(); cit != _ocells[i].end(); cit++){
int Simplex::kappa(Cell* tau) const{ Cell* cell = *cit;
for(int i=0; i < tau->getNumVertices(); i++){ delete cell;
if( !(this->hasVertex(tau->getVertex(i)->getNum())) ) return 0;
} }
_ocells[i].clear();
if(this->getDim() - tau->getDim() != 1) return 0;
int value=1;
for(int i=0; i < tau->getNumVertices(); i++){
if(this->getSortedVertex(i) != tau->getSortedVertex(i)) return value;
value = value*-1;
} }
//emptyTrash();
return value;
} }
*/
/*
int Cell::kappa(Cell* tau) const{
for(int i=0; i < tau->getNumVertices(); i++){
if( !(this->hasVertex(tau->getVertex(i)->getNum())) ) return 0;
}
if(this->getDim() - tau->getDim() != 1) return 0;
int value=1;
for(int i = 0; i < this->getNumFacets(); i++){
std::vector<MVertex*> vTau = tau->getVertexVector();
std::vector<MVertex*> v;
this->getFacetVertices(i, v);
value = -1;
if(v.size() != vTau.size()) printf("Error: invalid facet!");
do { /*
value = value*-1; CellComplex::CellComplex(CellComplex* cellComplex){
if(v == vTau) return value;
}
while (std::next_permutation(vTau.begin(), vTau.end()) );
vTau = tau->getVertexVector();
value = -1;
do {
value = value*-1;
if(v == vTau) return value;
}
while (std::prev_permutation(vTau.begin(), vTau.end()) );
_domain = cellComplex->getDomain();
_subdomain = cellComplex->getSubdomain;
_boundary = cellComplex->getBoundary;
_domainVertices = cellComplex->getDomainVertices;
for(int i = 0; i < 4; i++){
_betti[i] = cellComplex->getBetti(i);
_cells[i] = ocells[i];
_ocells[i] = ocells[i];
} }
return 0; _dim = cellComplex->getDim();
} }
*/
void CellComplex::insertCell(Cell* cell){
int OneSimplex::kappa(Cell* tau) const{ _ecells.insert(cell);
for(int i=0; i < tau->getNumVertices(); i++){
if( !(this->hasVertex(tau->getVertex(i)->getNum())) ) return 0;
} }
if(tau->getDim() != 0) return 0;
if(tau->getVertex(0) == this->getVertex(0)) return -1;
else return 1;
}*/
void CellComplex::removeCell(Cell* cell, bool other){ void CellComplex::removeCell(Cell* cell, bool other){
std::map<Cell*, int, Less_Cell > coboundary = cell->getOrientedCoboundary(); std::map<Cell*, int, Less_Cell > coboundary = cell->getOrientedCoboundary();
std::map<Cell*, int, Less_Cell > boundary = cell->getOrientedBoundary(); std::map<Cell*, int, Less_Cell > boundary = cell->getOrientedBoundary();
//std::list<Cell*> boundary = cell->getBoundary();
//std::list<Cell*> coboundary = cell->getCoboundary();
for(std::map<Cell*, int, Less_Cell>::iterator it = coboundary.begin(); it != coboundary.end(); it++){ for(std::map<Cell*, int, Less_Cell>::iterator it = coboundary.begin(); it != coboundary.end(); it++){
//for(std::list<Cell*>::iterator it = coboundary.begin(); it != coboundary.end(); it++){
Cell* cbdCell = (*it).first; Cell* cbdCell = (*it).first;
//Cell* cbdCell = *it;
cbdCell->removeBoundaryCell(cell, other); cbdCell->removeBoundaryCell(cell, other);
} }
for(std::map<Cell*, int, Less_Cell>::iterator it = boundary.begin(); it != boundary.end(); it++){ for(std::map<Cell*, int, Less_Cell>::iterator it = boundary.begin(); it != boundary.end(); it++){
//for(std::list<Cell*>::iterator it = boundary.begin(); it != boundary.end(); it++){
Cell* bdCell = (*it).first; Cell* bdCell = (*it).first;
//Cell* bdCell = *it;
bdCell->removeCoboundaryCell(cell, other); bdCell->removeCoboundaryCell(cell, other);
} }
...@@ -527,7 +483,7 @@ int CellComplex::reduceComplex(){ ...@@ -527,7 +483,7 @@ int CellComplex::reduceComplex(){
double t1 = Cpu(); double t1 = Cpu();
printf("Cell complex before reduction: %d volumes, %d faces, %d edges and %d vertices.\n", 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; int count = 0;
...@@ -535,23 +491,15 @@ int CellComplex::reduceComplex(){ ...@@ -535,23 +491,15 @@ int CellComplex::reduceComplex(){
int omitted = 0; int omitted = 0;
_store.clear(); _store.clear();
//if(omit > getDim()) omit = getDim();
removeSubdomain();
CellComplex::removeSubdomain();
//std::set<Cell*, Less_Cell> generatorCells;
while (getSize(getDim()) != 0){ while (getSize(getDim()) != 0){
citer cit = firstCell(getDim()); citer cit = firstCell(getDim());
Cell* cell = *cit; Cell* cell = *cit;
//generatorCells.insert(cell);
removeCell(cell); removeCell(cell);
//makeDualComplex();
//coreduction(cell);
//makeDualComplex();
omitted++; omitted++;
std::set< Cell*, Less_Cell > omittedCells; std::set< Cell*, Less_Cell > omittedCells;
_store.push_back(omittedCells); _store.push_back(omittedCells);
...@@ -561,17 +509,9 @@ int CellComplex::reduceComplex(){ ...@@ -561,17 +509,9 @@ int CellComplex::reduceComplex(){
} }
/*
for(citer cit = generatorCells.begin(); cit != generatorCells.end(); cit++){
Cell* cell = *cit;
cell->clearBoundary();
cell->clearCoboundary();
_cells[cell->getDim()].insert(cell);
}*/
double t2 = Cpu(); double t2 = Cpu();
printf("Cell complex after reduction: %d volumes, %d faces, %d edges and %d vertices (%g s).\n", 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); getSize(3), getSize(2), getSize(1), getSize(0), t2 - t1);
/* /*
...@@ -611,7 +551,7 @@ void CellComplex::removeSubdomain(){ ...@@ -611,7 +551,7 @@ void CellComplex::removeSubdomain(){
int CellComplex::coreduceComplex(){ int CellComplex::coreduceComplex(){
printf("Cell complex before coreduction: %d volumes, %d faces, %d edges and %d vertices.\n", Msg::Debug("Cell complex before 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));
int count = 0; int count = 0;
...@@ -631,11 +571,10 @@ int CellComplex::coreduceComplex(){ ...@@ -631,11 +571,10 @@ int CellComplex::coreduceComplex(){
int omitted = 0; int omitted = 0;
_store.clear(); _store.clear();
//std::set<Cell*, Less_Cell> generatorCells;
while (getSize(0) != 0){ while (getSize(0) != 0){
citer cit = firstCell(0); citer cit = firstCell(0);
Cell* cell = *cit; Cell* cell = *cit;
//generatorCells.insert(cell);
removeCell(cell, false); removeCell(cell, false);
std::set< Cell*, Less_Cell > omittedCells; std::set< Cell*, Less_Cell > omittedCells;
omitted++; omitted++;
...@@ -643,6 +582,7 @@ int CellComplex::coreduceComplex(){ ...@@ -643,6 +582,7 @@ int CellComplex::coreduceComplex(){
_store.at(omitted-1).insert(cell); _store.at(omitted-1).insert(cell);
coreduction(cell, omitted); coreduction(cell, omitted);
} }
/* /*
for(citer cit = generatorCells.begin(); cit != generatorCells.end(); cit++){ for(citer cit = generatorCells.begin(); cit != generatorCells.end(); cit++){
Cell* cell = *cit; Cell* cell = *cit;
...@@ -653,7 +593,7 @@ int CellComplex::coreduceComplex(){ ...@@ -653,7 +593,7 @@ int CellComplex::coreduceComplex(){
*/ */
printf("Cell complex after coreduction: %d volumes, %d faces, %d edges and %d vertices.\n", 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; return omitted;
...@@ -663,7 +603,7 @@ void CellComplex::computeBettiNumbers(){ ...@@ -663,7 +603,7 @@ void CellComplex::computeBettiNumbers(){
for(int i = 0; i < 4; i++){ for(int i = 0; i < 4; i++){
_betti[i] = 0; _betti[i] = 0;
printf("Betti number computation process: step %d of 4 \n", i+1); Msg::Debug("Betti number computation process: step %d of 4 \n", i+1);
while (getSize(i) != 0){ while (getSize(i) != 0){
citer cit = firstCell(i); citer cit = firstCell(i);
...@@ -677,7 +617,7 @@ void CellComplex::computeBettiNumbers(){ ...@@ -677,7 +617,7 @@ void CellComplex::computeBettiNumbers(){
coreduction(cell); coreduction(cell);
} }
} }
printf("Cell complex Betti numbers: \nH0 = %d \nH1 = %d \nH2 = %d \nH3 = %d \n", 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; return;
...@@ -687,7 +627,7 @@ void CellComplex::computeBettiNumbers(){ ...@@ -687,7 +627,7 @@ void CellComplex::computeBettiNumbers(){
int CellComplex::cocombine(int dim){ int CellComplex::cocombine(int dim){
double t1 = Cpu(); double t1 = Cpu();
printf("Cell complex before cocombining: %d volumes, %d faces, %d edges and %d vertices.\n", Msg::Debug("Cell complex before cocombining: %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));
if(dim < 0 || dim > 2) return 0; if(dim < 0 || dim > 2) return 0;
...@@ -735,7 +675,7 @@ int CellComplex::cocombine(int dim){ ...@@ -735,7 +675,7 @@ int CellComplex::cocombine(int dim){
removeCell(c1); removeCell(c1);
removeCell(c2); removeCell(c2);
std::pair<citer, bool> insertInfo = _cells[dim].insert(newCell); std::pair<citer, bool> insertInfo = _cells[dim].insert(newCell);
if(!insertInfo.second) printf("Warning: Combined cell not inserted! \n"); if(!insertInfo.second) Msg::Debug("Warning: Combined cell not inserted! \n");
cit = firstCell(dim); cit = firstCell(dim);
count++; count++;
...@@ -746,7 +686,7 @@ int CellComplex::cocombine(int dim){ ...@@ -746,7 +686,7 @@ int CellComplex::cocombine(int dim){
} }
} }
double t2 = Cpu(); double t2 = Cpu();
printf("Cell complex after cocombining: %d volumes, %d faces, %d edges and %d vertices (%g s).\n", 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); getSize(3), getSize(2), getSize(1), getSize(0), t2-t1);
return count; return count;
...@@ -754,7 +694,7 @@ int CellComplex::cocombine(int dim){ ...@@ -754,7 +694,7 @@ int CellComplex::cocombine(int dim){
int CellComplex::combine(int dim){ int CellComplex::combine(int dim){
double t1 = Cpu(); double t1 = Cpu();
printf("Cell complex before combining: %d volumes, %d faces, %d edges and %d vertices.\n", Msg::Debug("Cell complex before combining: %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));
if(dim < 1 || dim > 3) return 0; if(dim < 1 || dim > 3) return 0;
...@@ -802,7 +742,7 @@ int CellComplex::combine(int dim){ ...@@ -802,7 +742,7 @@ int CellComplex::combine(int dim){
removeCell(c1); removeCell(c1);
removeCell(c2); removeCell(c2);
std::pair<citer, bool> insertInfo = _cells[dim].insert(newCell); std::pair<citer, bool> insertInfo = _cells[dim].insert(newCell);
if(!insertInfo.second) printf("Warning: Combined cell not inserted! \n"); if(!insertInfo.second) Msg::Debug("Warning: Combined cell not inserted! \n");
cit = firstCell(dim); cit = firstCell(dim);
count++; count++;
} }
...@@ -812,7 +752,7 @@ int CellComplex::combine(int dim){ ...@@ -812,7 +752,7 @@ int CellComplex::combine(int dim){
} }
} }
double t2 = Cpu(); double t2 = Cpu();
printf("Cell complex after combining: %d volumes, %d faces, %d edges and %d vertices (%g s).\n", 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); getSize(3), getSize(2), getSize(1), getSize(0), t2 - t1);
...@@ -820,7 +760,7 @@ int CellComplex::combine(int dim){ ...@@ -820,7 +760,7 @@ int CellComplex::combine(int dim){
} }
/* /*
int CellComplex::combine(int dim){ int CellComplex::combine(int dim){ // version of combine that doesn't have a queue
double t1 = Cpu(); double t1 = Cpu();
printf("Cell complex before combining: %d volumes, %d faces, %d edges and %d vertices.\n", printf("Cell complex before combining: %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));
...@@ -901,21 +841,16 @@ void CellComplex::swapSubdomain(){ ...@@ -901,21 +841,16 @@ void CellComplex::swapSubdomain(){
int CellComplex::writeComplexMSH(const std::string &name){ int CellComplex::writeComplexMSH(const std::string &name){
FILE *fp = fopen(name.c_str(), "w"); FILE *fp = fopen(name.c_str(), "w");
if(!fp){ if(!fp){
Msg::Error("Unable to open file '%s'", name.c_str()); Msg::Error("Unable to open file '%s'", name.c_str());
printf("Unable to open file."); Msg::Debug("Unable to open file.");
return 0; return 0;
} }
fprintf(fp, "$MeshFormat\n2.1 0 8\n$EndMeshFormat\n"); fprintf(fp, "$MeshFormat\n2.1 0 8\n$EndMeshFormat\n");
fprintf(fp, "$Nodes\n"); fprintf(fp, "$Nodes\n");
fprintf(fp, "%d\n", (int)_domainVertices.size()); fprintf(fp, "%d\n", (int)_domainVertices.size());
for(std::set<MVertex*, Less_MVertex>::iterator vit = _domainVertices.begin(); vit != _domainVertices.end(); vit++){ for(std::set<MVertex*, Less_MVertex>::iterator vit = _domainVertices.begin(); vit != _domainVertices.end(); vit++){
...@@ -1073,14 +1008,14 @@ bool CellComplex::checkCoherence(){ ...@@ -1073,14 +1008,14 @@ bool CellComplex::checkCoherence(){
int ori = (*it).second; int ori = (*it).second;
citer cit = _cells[bdCell->getDim()].find(bdCell); citer cit = _cells[bdCell->getDim()].find(bdCell);
if(cit == lastCell(bdCell->getDim())){ if(cit == lastCell(bdCell->getDim())){
printf("Warning! Boundary cell not in cell complex! Boundary removed. \n"); Msg::Debug("Warning! Boundary cell not in cell complex! Boundary removed. \n");
//printf(" "); cell->printCell(); //printf(" "); cell->printCell();
//printf(" "); bdCell->printCell(); //printf(" "); bdCell->printCell();
cell->removeBoundaryCell(bdCell); cell->removeBoundaryCell(bdCell);
coherent = false; coherent = false;
} }
if(!bdCell->hasCoboundary(cell)){ if(!bdCell->hasCoboundary(cell)){
printf("Warning! Incoherent boundary/coboundary pair! Fixed. \n"); Msg::Debug("Warning! Incoherent boundary/coboundary pair! Fixed. \n");
//printf(" "); cell->printCell(); //printf(" "); cell->printCell();
//printf(" "); cell->printBoundary(); //printf(" "); cell->printBoundary();
//printf(" "); bdCell->printCell(); //printf(" "); bdCell->printCell();
...@@ -1096,14 +1031,14 @@ bool CellComplex::checkCoherence(){ ...@@ -1096,14 +1031,14 @@ bool CellComplex::checkCoherence(){
int ori = (*it).second; int ori = (*it).second;
citer cit = _cells[cbdCell->getDim()].find(cbdCell); citer cit = _cells[cbdCell->getDim()].find(cbdCell);
if(cit == lastCell(cbdCell->getDim())){ if(cit == lastCell(cbdCell->getDim())){
printf("Warning! Coboundary cell not in cell complex! Coboundary removed. \n"); Msg::Debug("Warning! Coboundary cell not in cell complex! Coboundary removed. \n");
//printf(" "); cell->printCell(); //printf(" "); cell->printCell();
//printf(" "); cbdCell->printCell(); //printf(" "); cbdCell->printCell();
cell->removeCoboundaryCell(cbdCell); cell->removeCoboundaryCell(cbdCell);
coherent = false; coherent = false;
} }
if(!cbdCell->hasBoundary(cell)){ if(!cbdCell->hasBoundary(cell)){
printf("Warning! Incoherent coboundary/boundary pair! Fixed. \n"); Msg::Debug("Warning! Incoherent coboundary/boundary pair! Fixed. \n");
//printf(" "); cell->printCell(); //printf(" "); cell->printCell();
//printf(" "); cell->printCoboundary(); //printf(" "); cell->printCoboundary();
//printf(" "); cbdCell->printCell(); //printf(" "); cbdCell->printCell();
...@@ -1119,24 +1054,55 @@ bool CellComplex::checkCoherence(){ ...@@ -1119,24 +1054,55 @@ bool CellComplex::checkCoherence(){
return coherent; return coherent;
} }
void CellComplex::emptyTrash(){
bool Less_Cell::operator()(const Cell* c1, const Cell* c2) const { for(std::list<Cell*>::iterator cit = _trash.begin(); cit != _trash.end(); cit++){
Cell* cell = *cit;
//cells with fever vertices first delete cell;
}
if(c1->getNumVertices() != c2->getNumVertices()){ _trash.clear();
return (c1->getNumVertices() < c2->getNumVertices());
} }
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();
_ecells.clear();
_trash.clear();
}
bool CellComplex::hasCell(Cell* cell, bool org){
if(!org){
citer cit = _cells[cell->getDim()].find(cell);
if( cit == lastCell(cell->getDim()) ) return false;
else return true;
}
else{
citer cit = _ocells[cell->getDim()].find(cell);
if( cit == lastOrgCell(cell->getDim()) ) return false;
else return true;
}
}
//"natural number" -like ordering (the number of a vertice corresponds a digit) 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 < c1->getNumVertices();i++){ for(int i = 0; i < 4; i++){
if(c1->getSortedVertex(i) < c2->getSortedVertex(i)) return true; for(citer cit = firstCell(i); cit != lastCell(i); cit++){
else if (c1->getSortedVertex(i) > c2->getSortedVertex(i)) return false; Cell* cell = *cit;
cell->makeDualCell();
} }
return false;
} }
}
#endif #endif
Geo/CellComplex.h
+ 51
− 396
View file @ 61a9dafb
This diff is collapsed.
Geo/ChainComplex.cpp
+ 69
− 10
View file @ 61a9dafb
...@@ -82,17 +82,15 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){ ...@@ -82,17 +82,15 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){
//printf("cell1: %d, cell2: %d \n", bdCell->getIndex(), cell->getIndex()); //printf("cell1: %d, cell2: %d \n", bdCell->getIndex(), cell->getIndex());
if(bdCell->getIndex() > (int)gmp_matrix_rows( _HMatrix[dim]) || bdCell->getIndex() < 1 if(bdCell->getIndex() > (int)gmp_matrix_rows( _HMatrix[dim]) || bdCell->getIndex() < 1
|| cell->getIndex() > (int)gmp_matrix_cols( _HMatrix[dim]) || cell->getIndex() < 1){ || cell->getIndex() > (int)gmp_matrix_cols( _HMatrix[dim]) || cell->getIndex() < 1){
printf("Warning: Index out of bound! HMatrix: %d. \n", dim); Msg::Debug("Warning: Index out of bound! HMatrix: %d. \n", dim);
} }
else{ else{
gmp_matrix_get_elem(elem, bdCell->getIndex(), cell->getIndex(), _HMatrix[dim]); gmp_matrix_get_elem(elem, bdCell->getIndex(), cell->getIndex(), _HMatrix[dim]);
old_elem = mpz_get_si(elem); old_elem = mpz_get_si(elem);
mpz_set_si(elem, old_elem + (*it).second); mpz_set_si(elem, old_elem + (*it).second);
/*if( (old_elem + (*it).second) > 1 || (old_elem + (*it).second) < -1 ){ if( abs((old_elem + (*it).second)) > 1){
printf("Warning: Invalid incidence index: %d! HMatrix: %d.", (old_elem + (*it).second), dim); Msg::Debug("Incidence index: %d! HMatrix: %d.", (old_elem + (*it).second), dim);
printf(" Set to %d. \n", (old_elem + (*it).second) % 2); }
mpz_set_si(elem, (old_elem + (*it).second) % 2);
}*/
gmp_matrix_set_elem(elem, bdCell->getIndex(), cell->getIndex(), _HMatrix[dim]); gmp_matrix_set_elem(elem, bdCell->getIndex(), cell->getIndex(), _HMatrix[dim]);
} }
} }
...@@ -291,7 +289,7 @@ void ChainComplex::computeHomology(bool dual){ ...@@ -291,7 +289,7 @@ void ChainComplex::computeHomology(bool dual){
//KerCod(highDim); //KerCod(highDim);
} }
printf("Homology computation process: step %d of 4 \n", i+1); Msg::Debug("Homology computation process: step %d of 4 \n", i+1);
KerCod(highDim); KerCod(highDim);
...@@ -560,7 +558,7 @@ void Chain::smoothenChain(){ ...@@ -560,7 +558,7 @@ void Chain::smoothenChain(){
eraseNullCells(); eraseNullCells();
double t2 = Cpu(); double t2 = Cpu();
printf("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; return;
} }
...@@ -574,7 +572,7 @@ int Chain::writeChainMSH(const std::string &name){ ...@@ -574,7 +572,7 @@ int Chain::writeChainMSH(const std::string &name){
FILE *fp = fopen(name.c_str(), "a"); FILE *fp = fopen(name.c_str(), "a");
if(!fp){ if(!fp){
Msg::Error("Unable to open file '%s'", name.c_str()); Msg::Error("Unable to open file '%s'", name.c_str());
printf("Unable to open file."); Msg::Debug("Unable to open file.");
return 0; return 0;
} }
...@@ -674,4 +672,65 @@ void Chain::createPView(){ ...@@ -674,4 +672,65 @@ void Chain::createPView(){
} }
void Chain::removeCell(Cell* cell) {
citer it = _cells.find(cell);
if(it != _cells.end()){
(*it).second = 0;
}
return;
}
void Chain::addCell(Cell* cell, int coeff) {
std::pair<citer,bool> insert = _cells.insert( std::make_pair( cell, coeff));
if(!insert.second && (*insert.first).second == 0) (*insert.first).second = coeff;
else if (!insert.second && (*insert.first).second != 0) Msg::Debug("Error: invalid chain smoothening add! \n");
if(!_cellComplex->hasCell(cell)){
_cellComplex->insertCell(cell);
}
return;
}
bool Chain::hasCell(Cell* c){
citer it = _cells.find(c);
if(it != _cells.end() && (*it).second != 0) return true;
return false;
}
Cell* Chain::findCell(Cell* c){
citer it = _cells.find(c);
if(it != _cells.end() && (*it).second != 0) return (*it).first;
return NULL;
}
int Chain::getCoeff(Cell* c){
citer it = _cells.find(c);
if(it != _cells.end()) return (*it).second;
return 0;
}
void Chain::eraseNullCells(){
for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
if( (*cit).second == 0){
//cit++;
//_cells.erase(--cit);
_cells.erase(cit);
++cit;
}
}
for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
if( (*cit).second == 0){
_cells.erase(cit);
cit = _cells.begin();
}
}
return;
}
void Chain::deImmuneCells(){
for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
Cell* cell = (*cit).first;
cell->setImmune(false);
}
}
#endif #endif
Geo/ChainComplex.h
+ 10
− 70
View file @ 61a9dafb
...@@ -115,7 +115,7 @@ class ChainComplex{ ...@@ -115,7 +115,7 @@ class ChainComplex{
// Compute bases for the homology groups of this chain complex // Compute bases for the homology groups of this chain complex
void computeHomology(bool dual=false); void computeHomology(bool dual=false);
// get coefficient vector for dim-dimensional chain chainNumber // get coefficient vector for dim-dimensional chain chainNumber (columns of _Hbasis[dim])
std::vector<int> getCoeffVector(int dim, int chainNumber); std::vector<int> getCoeffVector(int dim, int chainNumber);
// torsion coefficient for dim-dimensional chain chainNumber // torsion coefficient for dim-dimensional chain chainNumber
int getTorsion(int dim, int chainNumber); int getTorsion(int dim, int chainNumber);
...@@ -171,48 +171,15 @@ class Chain{ ...@@ -171,48 +171,15 @@ class Chain{
typedef std::map<Cell*, int, Less_Cell>::iterator citer; typedef std::map<Cell*, int, Less_Cell>::iterator citer;
// get i:th cell of this chain
//Cell* getCell(int i) { return _cells.at(i).first; }
// get coeffcient of i:th cell of this chain
//int getCoeff(int i) { return _cells.at(i).second; }
// remove a cell from this chain // remove a cell from this chain
void removeCell(Cell* cell) { void removeCell(Cell* cell);
citer it = _cells.find(cell);
if(it != _cells.end()){
(*it).second = 0;
}
return;
}
// add a cell to this chain // add a cell to this chain
void addCell(Cell* cell, int coeff) { void 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) printf("Error: invalid chain smoothening add! \n");
if(!_cellComplex->hasCell(cell)){
_cellComplex->insertCell(cell);
}
return;
}
bool hasCell(Cell* c){ bool hasCell(Cell* c);
citer it = _cells.find(c); Cell* findCell(Cell* c);
if(it != _cells.end() && (*it).second != 0) return true; int getCoeff(Cell* c);
return false;
}
Cell* findCell(Cell* c){
citer it = _cells.find(c);
if(it != _cells.end() && (*it).second != 0) return (*it).first;
return NULL;
}
int getCoeff(Cell* c){
citer it = _cells.find(c);
if(it != _cells.end()) return (*it).second;
return 0;
}
int getTorsion() {return _torsion;} int getTorsion() {return _torsion;}
...@@ -222,40 +189,12 @@ class Chain{ ...@@ -222,40 +189,12 @@ class Chain{
std::map<Cell*, int, Less_Cell> getCells() {return _cells;} std::map<Cell*, int, Less_Cell> getCells() {return _cells;}
// erase cells from the chain with zero coefficient // erase cells from the chain with zero coefficient
void eraseNullCells(){ void eraseNullCells();
for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
if( (*cit).second == 0){
//cit++;
//_cells.erase(--cit);
_cells.erase(cit);
++cit;
}
}
for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
if( (*cit).second == 0){
_cells.erase(cit);
cit = _cells.begin();
}
}
return;
}
void deImmuneCells();
void deImmuneCells(){
for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
Cell* cell = (*cit).first;
cell->setImmune(false);
}
}
// number of cells in this chain // number of cells in this chain
int getSize() { int getSize() { return _cells.size();}
//eraseNullCells();
return _cells.size();
}
int getNumCells() {
return _cells.size();
}
// get/set chain name // get/set chain name
std::string getName() { return _name; } std::string getName() { return _name; }
...@@ -269,6 +208,7 @@ class Chain{ ...@@ -269,6 +208,7 @@ class Chain{
void smoothenChain(); void smoothenChain();
// append this chain to a MSH ASCII file as $ElementData // append this chain to a MSH ASCII file as $ElementData
// for debugging only
int writeChainMSH(const std::string &name); int writeChainMSH(const std::string &name);
// create a PView of this chain. // create a PView of this chain.
......
Geo/Homology.cpp
+ 106
− 18
View file @ 61a9dafb
...@@ -19,6 +19,8 @@ Homology::Homology(GModel* model, std::vector<int> physicalDomain, std::vector<i ...@@ -19,6 +19,8 @@ Homology::Homology(GModel* model, std::vector<int> physicalDomain, std::vector<i
_subdomain = physicalSubdomain; _subdomain = physicalSubdomain;
Msg::Info("Creating a Cell Complex..."); Msg::Info("Creating a Cell Complex...");
Msg::StatusBar(1, false, "Cell Complex...");
Msg::StatusBar(2, false, "");
double t1 = Cpu(); double t1 = Cpu();
std::map<int, std::vector<GEntity*> > groups[4]; std::map<int, std::vector<GEntity*> > groups[4];
...@@ -68,13 +70,25 @@ Homology::Homology(GModel* model, std::vector<int> physicalDomain, std::vector<i ...@@ -68,13 +70,25 @@ Homology::Homology(GModel* model, std::vector<int> physicalDomain, std::vector<i
Msg::Info("Cell Complex complete (%g s).", t2 - t1); Msg::Info("Cell Complex complete (%g s).", t2 - t1);
Msg::Info("%d volumes, %d faces, %d edges and %d vertices.", 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));
} }
Homology::~Homology(){
delete _cellComplex;
for(int i = 0; i < 4; i++) {
for(int j = 0; j < _generators[i].size(); j++){
Chain* chain = _generators[i].at(j);
//_model->deletePhysicalGroup(chain->getDim(), chain->getNum());
delete chain;
}
}
}
void Homology::findGenerators(std::string fileName){ void Homology::findGenerators(std::string fileName){
Msg::Info("Reducing Cell Complex..."); Msg::Info("Reducing the Cell Complex...");
Msg::StatusBar(1, false, "Reducing...");
double t1 = Cpu(); double t1 = Cpu();
//_cellComplex->printEuler(); //_cellComplex->printEuler();
int omitted = _cellComplex->reduceComplex(); int omitted = _cellComplex->reduceComplex();
...@@ -98,10 +112,12 @@ void Homology::findGenerators(std::string fileName){ ...@@ -98,10 +112,12 @@ void Homology::findGenerators(std::string fileName){
Msg::Info("Cell Complex reduction complete (%g s).", t2 - t1); Msg::Info("Cell Complex reduction complete (%g s).", t2 - t1);
Msg::Info("%d volumes, %d faces, %d edges and %d vertices.", 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));
//for(int i = 0; i < 4; i++) { printf("Dim %d: \n", i); _cellComplex->printComplex(i); } //for(int i = 0; i < 4; i++) { printf("Dim %d: \n", i); _cellComplex->printComplex(i); }
Msg::Info("Computing homology groups..."); Msg::Info("Computing homology spaces...");
Msg::StatusBar(1, false, "Computing...");
t1 = Cpu(); t1 = Cpu();
ChainComplex* chains = new ChainComplex(_cellComplex); ChainComplex* chains = new ChainComplex(_cellComplex);
chains->computeHomology(); chains->computeHomology();
...@@ -149,19 +165,22 @@ void Homology::findGenerators(std::string fileName){ ...@@ -149,19 +165,22 @@ void Homology::findGenerators(std::string fileName){
createPViews(); createPViews();
if(fileName != "") writeGeneratorsMSH(fileName); if(fileName != "") writeGeneratorsMSH(fileName);
Msg::Info("Ranks of homology groups for primal cell complex:"); Msg::Info("Ranks of homology spaces for primal cell complex:");
Msg::Info("H0 = %d", HRank[0]); Msg::Info("H0 = %d", HRank[0]);
Msg::Info("H1 = %d", HRank[1]); Msg::Info("H1 = %d", HRank[1]);
Msg::Info("H2 = %d", HRank[2]); Msg::Info("H2 = %d", HRank[2]);
Msg::Info("H3 = %d", HRank[3]); Msg::Info("H3 = %d", HRank[3]);
if(omitted != 0) Msg::Info("The computation of generators in %d highest dimensions was omitted.", omitted); if(omitted != 0) Msg::Info("The computation of generators in the highest dimension was omitted.");
delete chains; delete chains;
printf("H0 = %d \n", HRank[0]); Msg::Debug("H0 = %d \n", HRank[0]);
printf("H1 = %d \n", HRank[1]); Msg::Debug("H1 = %d \n", HRank[1]);
printf("H2 = %d \n", HRank[2]); Msg::Debug("H2 = %d \n", HRank[2]);
printf("H3 = %d \n", HRank[3]); 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]);
return; return;
} }
...@@ -171,6 +190,7 @@ void Homology::findDualGenerators(std::string fileName){ ...@@ -171,6 +190,7 @@ void Homology::findDualGenerators(std::string fileName){
//for(int i = 0; i < 4; i++) { printf("Dim %d: \n", i); _cellComplex->printComplex(i); } //for(int i = 0; i < 4; i++) { printf("Dim %d: \n", i); _cellComplex->printComplex(i); }
Msg::Info("Reducing Cell Complex..."); Msg::Info("Reducing Cell Complex...");
Msg::StatusBar(1, false, "Reducing...");
double t1 = Cpu(); double t1 = Cpu();
int omitted = _cellComplex->coreduceComplex(); int omitted = _cellComplex->coreduceComplex();
//_cellComplex->emptyTrash(); //_cellComplex->emptyTrash();
...@@ -199,12 +219,14 @@ void Homology::findDualGenerators(std::string fileName){ ...@@ -199,12 +219,14 @@ void Homology::findDualGenerators(std::string fileName){
Msg::Info("Cell Complex reduction complete (%g s).", t2 - t1); Msg::Info("Cell Complex reduction complete (%g s).", t2 - t1);
Msg::Info("%d volumes, %d faces, %d edges and %d vertices.", 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->writeComplexMSH(fileName);
_cellComplex->writeComplexMSH(fileName); Msg::Info("Computing homology spaces...");
Msg::StatusBar(1, false, "Computing...");
Msg::Info("Computing homology groups...");
t1 = Cpu(); t1 = Cpu();
ChainComplex* chains = new ChainComplex(_cellComplex); ChainComplex* chains = new ChainComplex(_cellComplex);
chains->transposeHMatrices(); chains->transposeHMatrices();
...@@ -255,7 +277,7 @@ void Homology::findDualGenerators(std::string fileName){ ...@@ -255,7 +277,7 @@ void Homology::findDualGenerators(std::string fileName){
createPViews(); createPViews();
if(fileName != "") writeGeneratorsMSH(fileName); if(fileName != "") writeGeneratorsMSH(fileName);
Msg::Info("Ranks of homology groups for dual cell complex:"); Msg::Info("Ranks of homology spaces for dual cell complex:");
Msg::Info("H0* = %d", HRank[0]); Msg::Info("H0* = %d", HRank[0]);
Msg::Info("H1* = %d", HRank[1]); Msg::Info("H1* = %d", HRank[1]);
Msg::Info("H2* = %d", HRank[2]); Msg::Info("H2* = %d", HRank[2]);
...@@ -264,10 +286,13 @@ void Homology::findDualGenerators(std::string fileName){ ...@@ -264,10 +286,13 @@ void Homology::findDualGenerators(std::string fileName){
delete chains; delete chains;
printf("H0* = %d \n", HRank[0]); Msg::Debug("H0* = %d \n", HRank[0]);
printf("H1* = %d \n", HRank[1]); Msg::Debug("H1* = %d \n", HRank[1]);
printf("H2* = %d \n", HRank[2]); Msg::Debug("H2* = %d \n", HRank[2]);
printf("H3* = %d \n", HRank[3]); 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]);
return; return;
} }
...@@ -275,6 +300,7 @@ void Homology::findDualGenerators(std::string fileName){ ...@@ -275,6 +300,7 @@ void Homology::findDualGenerators(std::string fileName){
void Homology::computeBettiNumbers(){ void Homology::computeBettiNumbers(){
Msg::Info("Running coreduction..."); Msg::Info("Running coreduction...");
Msg::StatusBar(1, false, "Computing...");
double t1 = Cpu(); double t1 = Cpu();
_cellComplex->computeBettiNumbers(); _cellComplex->computeBettiNumbers();
double t2 = Cpu(); double t2 = Cpu();
...@@ -285,7 +311,69 @@ void Homology::computeBettiNumbers(){ ...@@ -285,7 +311,69 @@ void Homology::computeBettiNumbers(){
Msg::Info("H2 = %d", _cellComplex->getBettiNumber(2)); Msg::Info("H2 = %d", _cellComplex->getBettiNumber(2));
Msg::Info("H3 = %d", _cellComplex->getBettiNumber(3)); 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));
return; return;
} }
void Homology::restoreHomology() {
_cellComplex->restoreComplex();
for(int i = 0; i < 4; i++) _generators[i].clear();
}
std::string Homology::getDomainString() {
std::string domainString = "({";
for(unsigned int i = 0; i < _domain.size(); i++){
std::string temp = "";
convert(_domain.at(i),temp);
domainString += temp;
if (_domain.size()-1 > i){
domainString += ", ";
}
}
domainString += "}";
if(!_subdomain.empty()){
domainString += ", {";
for(unsigned int i = 0; i < _subdomain.size(); i++){
std::string temp = "";
convert(_subdomain.at(i),temp);
domainString += temp;
if (_subdomain.size()-1 > i){
domainString += ", ";
}
}
domainString += "}";
}
domainString += ") ";
return domainString;
}
void Homology::createPViews(){
for(int i = 0; i < 4; i++){
for(int j = 0; j < _generators[i].size(); j++){
Chain* chain = _generators[i].at(j);
chain->createPView();
}
}
}
bool Homology::writeGeneratorsMSH(std::string fileName, bool binary){
if(!_model->writeMSH(fileName, 2.0, binary)) return false;
/*
for(int i = 0; i < 4; i++){
for(int j = 0; j < _generators[i].size(); j++){
Chain* chain = _generators[i].at(j);
if(!chain->writeChainMSH(fileName)) 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());
return true;
}
#endif #endif
Geo/Homology.h
+ 6
− 67
View file @ 61a9dafb
...@@ -41,16 +41,7 @@ class Homology ...@@ -41,16 +41,7 @@ class Homology
public: public:
Homology(GModel* model, std::vector<int> physicalDomain, std::vector<int> physicalSubdomain); Homology(GModel* model, std::vector<int> physicalDomain, std::vector<int> physicalSubdomain);
~Homology(){ ~Homology();
delete _cellComplex;
for(int i = 0; i < 4; i++) {
for(int j = 0; j < _generators[i].size(); j++){
Chain* chain = _generators[i].at(j);
//_model->deletePhysicalGroup(chain->getDim(), chain->getNum());
delete chain;
}
}
}
// Find the generators/duals of homology spaces, or just compute the ranks of homology spaces // Find the generators/duals of homology spaces, or just compute the ranks of homology spaces
void findGenerators(std::string fileName); void findGenerators(std::string fileName);
...@@ -61,67 +52,15 @@ class Homology ...@@ -61,67 +52,15 @@ class Homology
//void swapSubdomain() { _cellComplex->swapSubdomain(); } //void swapSubdomain() { _cellComplex->swapSubdomain(); }
// Restore the cell complex to its original state before cell reductions // Restore the cell complex to its original state before cell reductions
void restoreHomology() { void restoreHomology();
_cellComplex->restoreComplex();
for(int i = 0; i < 4; i++) _generators[i].clear();
}
// Create a string describing the generator // Create a string describing the generator
std::string getDomainString() { std::string getDomainString();
std::string domainString = "({";
for(unsigned int i = 0; i < _domain.size(); i++){
std::string temp = "";
convert(_domain.at(i),temp);
domainString += temp;
if (_domain.size()-1 > i){
domainString += ", ";
}
}
domainString += "}";
if(!_subdomain.empty()){
domainString += ", {";
for(unsigned int i = 0; i < _subdomain.size(); i++){
std::string temp = "";
convert(_subdomain.at(i),temp);
domainString += temp;
if (_subdomain.size()-1 > i){
domainString += ", ";
}
}
domainString += "}";
}
domainString += ") ";
return domainString;
}
// create PViews of the generators
void createPViews(){
for(int i = 0; i < 4; i++){
for(int j = 0; j < _generators[i].size(); j++){
Chain* chain = _generators[i].at(j);
chain->createPView();
}
}
}
// create PViews of the generators and save the chain mesh elements to the mesh of the model
void createPViews();
// write the generators to a file // write the generators to a file
bool writeGeneratorsMSH(std::string fileName, bool binary=false){ bool writeGeneratorsMSH(std::string fileName, bool binary=false);
if(!_model->writeMSH(fileName, 2.0, binary)) return false;
/*
for(int i = 0; i < 4; i++){
for(int j = 0; j < _generators[i].size(); j++){
Chain* chain = _generators[i].at(j);
if(!chain->writeChainMSH(fileName)) return false;
}
}*/
Msg::Info("Wrote homology computation results to %s.", fileName.c_str());
printf("Wrote homology computation results to %s. \n", fileName.c_str());
return true;
}
}; };
......
tutorial/t10.geo
+ 2
− 3
View file @ 61a9dafb
...@@ -133,9 +133,9 @@ Mesh 3; ...@@ -133,9 +133,9 @@ Mesh 3;
// Save the generator chains to t10_hom.msh. // Save the generator chains to t10_hom.msh.
HomGen("t10_hom.msh") = {{69}, {70, 71, 72, 73}}; HomGen("t10_hom.msh") = {{69}, {70, 71, 72, 73}};
// Find the corresponding cuts. // Find the corresponding thin cuts.
// Save the cut chains to t10_hom.msh. // Save the cut chains to t10_hom.msh.
HomCut("t10_hom.msh") = {{69}, {70, 71, 72, 73}}; HomGen("t10_hom.msh") = {{69}, {75}};
// Only find and print the ranks of the relative homology spaces (Betti numbers). // Only find and print the ranks of the relative homology spaces (Betti numbers).
HomRank {{69},{70, 71, 72, 73}}; HomRank {{69},{70, 71, 72, 73}};
...@@ -143,6 +143,5 @@ HomRank {{69},{70, 71, 72, 73}}; ...@@ -143,6 +143,5 @@ HomRank {{69},{70, 71, 72, 73}};
// More examples (uncomment): // More examples (uncomment):
// HomGen("t10_homgen.msh_1") = {{69}, {}}; // HomGen("t10_homgen.msh_1") = {{69}, {}};
// HomGen("t10_homgen.msh_2") = {{69}, {74}}; // HomGen("t10_homgen.msh_2") = {{69}, {74}};
// HomGen("t10_homgen.msh_3") = {{69}, {75}};
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment