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

Bugfixes to cell combining. I can still find examples where it fails.

parent 03e776f6
No related branches found
No related tags found
No related merge requests found
......@@ -105,11 +105,12 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su
cell->setTag(tag);
tag++;
}
//_originalCells[i] = _cells[i];
_cells2[i] = _cells[i];
_betti[i] = 0;
if(getSize(i) > _dim) _dim = i;
}
}
void CellComplex::insert_cells(bool subdomain, bool boundary){
......@@ -254,6 +255,8 @@ int OneSimplex::kappa(Cell* tau) const{
void CellComplex::removeCell(Cell* cell){
//_trash.insert(cell);
_cells[cell->getDim()].erase(cell);
std::list<Cell*> coboundary = cell->getCoboundary();
......@@ -269,7 +272,6 @@ void CellComplex::removeCell(Cell* cell){
bdCell->removeCoboundaryCell(cell);
}
//_trash.push_back(cell);
}
......@@ -332,7 +334,6 @@ int CellComplex::coreduction(Cell* generator){
int CellComplex::reduction(int dim){
if(dim < 1 || dim > 3) return 0;
std::list<Cell*> cbd_c;
int count = 0;
......@@ -343,7 +344,9 @@ int CellComplex::reduction(int dim){
Cell* cell = *cit;
cbd_c = cell->getCoboundary();
if(cbd_c.size() == 1 && inSameDomain(cell, cbd_c.front()) ){
removeCell(cell);
removeCell(cbd_c.front());
count++;
reduced = true;
......@@ -362,21 +365,30 @@ int CellComplex::reduceComplex(bool omitHighdim){
int count = 0;
for(int i = 3; i > 0; i--) count = count + reduction(i);
int omitted = 0;
if(count == 0 && omitHighdim){
removeSubdomain();
CellComplex::removeSubdomain();
CellComplex::removeSubdomain();
std::set<Cell*, Less_Cell> generatorCells;
while (getSize(getDim()) != 0){
citer cit = firstCell(getDim());
Cell* cell = *cit;
generatorCells.insert(cell);
removeCell(cell);
reduction(3);
reduction(2);
reduction(1);
omitted++;
}
......@@ -390,7 +402,6 @@ int CellComplex::reduceComplex(bool omitHighdim){
}
printf("Cell complex after reduction: %d volumes, %d faces, %d edges and %d vertices.\n",
getSize(3), getSize(2), getSize(1), getSize(0));
......@@ -422,7 +433,7 @@ int CellComplex::coreduceComplex(bool omitLowdim){
int count = 0;
CellComplex::removeSubdomain();
CellComplex::removeSubdomain();
for(int dim = 0; dim < 4; dim++){
citer cit = firstCell(dim);
......@@ -506,7 +517,7 @@ void CellComplex::replaceCells(Cell* c1, Cell* c2, Cell* newCell, bool orMatch,
Cell* cbdCell = (*it).second;
int ori = (*it).first;
cbdCell->removeBoundaryCell(c1);
cbdCell->addBoundaryCell(ori, newCell);
cbdCell->addBoundaryCell(ori, newCell, true);
}
for(std::list< std::pair<int, Cell*> >::iterator it = coboundary2.begin(); it != coboundary2.end(); it++){
Cell* cbdCell = (*it).second;
......@@ -519,16 +530,16 @@ void CellComplex::replaceCells(Cell* c1, Cell* c2, Cell* newCell, bool orMatch,
Cell* cell2 = (*it2).second;
if(*cell2 == *cbdCell) old = true;
}
if(!old) cbdCell->addBoundaryCell(ori, newCell);
if(!old) cbdCell->addBoundaryCell(ori, newCell, true);
}
else cbdCell->addBoundaryCell(ori, newCell);
else cbdCell->addBoundaryCell(ori, newCell, true);
}
for(std::list< std::pair<int, Cell* > >::iterator it = boundary1.begin(); it != boundary1.end(); it++){
Cell* bdCell = (*it).second;
int ori = (*it).first;
bdCell->removeCoboundaryCell(c1);
bdCell->addCoboundaryCell(ori, newCell);
bdCell->addCoboundaryCell(ori, newCell, true);
}
for(std::list< std::pair<int, Cell* > >::iterator it = boundary2.begin(); it != boundary2.end(); it++){
Cell* bdCell = (*it).second;
......@@ -541,14 +552,16 @@ void CellComplex::replaceCells(Cell* c1, Cell* c2, Cell* newCell, bool orMatch,
Cell* cell2 = (*it2).second;
if(*cell2 == *bdCell) old = true;
}
if(!old) bdCell->addCoboundaryCell(ori, newCell);
if(!old) bdCell->addCoboundaryCell(ori, newCell, true);
}
else bdCell->addCoboundaryCell(ori, newCell);
else bdCell->addCoboundaryCell(ori, newCell, true);
}
_cells[dim].erase(c1);
_cells[dim].erase(c2);
//_trash.insert(c1);
//_trash.insert(c2);
_cells[dim].insert(newCell);
......@@ -584,6 +597,8 @@ int CellComplex::cocombine(int dim){
if(bd_c.size() == 2 && !(*(bd_c.front().second) == *(bd_c.back().second))
&& inSameDomain(s, bd_c.front().second) && inSameDomain(s, bd_c.back().second) ){
int or1 = bd_c.front().first;
int or2 = bd_c.back().first;
Cell* c1 = bd_c.front().second;
Cell* c2 = bd_c.back().second;
......@@ -592,15 +607,13 @@ int CellComplex::cocombine(int dim){
cbd_c = c2->getCoboundary();
enqueueCells(cbd_c, Q, Qset);
int or1 = bd_c.front().first;
int or2 = bd_c.back().first;
removeCell(s);
CombinedCell* newCell = new CombinedCell(c1, c2, (or1 != or2), true );
replaceCells(c1, c2, newCell, (or1 != or2), true);
cit = firstCell(dim);
count++;
}
removeCellQset(s, Qset);
......@@ -641,7 +654,8 @@ int CellComplex::combine(int dim){
if(cbd_c.size() == 2 && !(*(cbd_c.front().second) == *(cbd_c.back().second))
&& inSameDomain(s, cbd_c.front().second) && inSameDomain(s, cbd_c.back().second) ){
int or1 = cbd_c.front().first;
int or2 = cbd_c.back().first;
Cell* c1 = cbd_c.front().second;
Cell* c2 = cbd_c.back().second;
......@@ -650,16 +664,13 @@ int CellComplex::combine(int dim){
bd_c = c2->getBoundary();
enqueueCells(bd_c, Q, Qset);
int or1 = cbd_c.front().first;
int or2 = cbd_c.back().first;
removeCell(s);
CombinedCell* newCell = new CombinedCell(c1, c2, (or1 != or2) );
replaceCells(c1, c2, newCell, (or1 != or2));
cit = firstCell(dim);
//cit++;
count++;
}
removeCellQset(s, Qset);
......
......@@ -88,6 +88,7 @@ class Cell
for( std::list< std::pair<int, Cell*> >::iterator it= _boundary.begin();it!= _boundary.end();it++){
Cell* cell = (*it).second;
boundary.push_back(cell);
if((*it).first == 0) boundary.push_back(cell);
}
return boundary;
}
......@@ -97,24 +98,48 @@ class Cell
for( std::list< std::pair<int, Cell*> >::iterator it= _coboundary.begin();it!= _coboundary.end();it++){
Cell* cell = (*it).second;
coboundary.push_back(cell);
if((*it).first == 0) coboundary.push_back(cell);
}
return coboundary;
}
virtual void addBoundaryCell(int orientation, Cell* cell) { _boundary.push_back( std::make_pair(orientation, cell)); }
virtual void addCoboundaryCell(int orientation, Cell* cell) { _coboundary.push_back( std::make_pair(orientation, cell)); }
virtual void addBoundaryCell(int orientation, Cell* cell, bool duplicates=false) {
if(!duplicates){
for(std::list< std::pair<int, Cell*> >::iterator it = _boundary.begin(); it != _boundary.end(); it++){
Cell* cell2 = (*it).second;
int ori2 = (*it).first;
if(*cell2 == *cell && !duplicates) return;
}
}
_boundary.push_back( std::make_pair(orientation, cell));
virtual void removeBoundaryCell(Cell* cell) {
}
virtual void addCoboundaryCell(int orientation, Cell* cell, bool duplicates=false) {
if(!duplicates){
for(std::list< std::pair<int, Cell*> >::iterator it = _coboundary.begin(); it != _coboundary.end(); it++){
Cell* cell2 = (*it).second;
int ori2 = (*it).first;
if(*cell2 == *cell && !duplicates) return;
}
}
_coboundary.push_back( std::make_pair(orientation, cell));
}
virtual int removeBoundaryCell(Cell* cell) {
int count = 0;
for(std::list< std::pair<int, Cell*> >::iterator it = _boundary.begin(); it != _boundary.end(); it++){
Cell* cell2 = (*it).second;
if(*cell2 == *cell) { _boundary.erase(it); break; }
if(*cell2 == *cell) { _boundary.erase(it); count++; it = _boundary.begin(); }
}
return count;
}
virtual void removeCoboundaryCell(Cell* cell) {
virtual int removeCoboundaryCell(Cell* cell) {
int count = 0;
for(std::list< std::pair<int, Cell*> >::iterator it = _coboundary.begin(); it != _coboundary.end(); it++){
Cell* cell2 = (*it).second;
if(*cell2 == *cell) { _coboundary.erase(it); break; }
if(*cell2 == *cell) { _coboundary.erase(it); count++; it = _coboundary.begin(); }
}
return count;
}
virtual void clearBoundary() { _boundary.clear(); }
......@@ -187,6 +212,7 @@ class Simplex : public Cell
public:
Simplex() : Cell() {
_combined = false;
_index = 0;
}
~Simplex(){}
......@@ -413,6 +439,7 @@ class CombinedCell : public Cell{
public:
CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co=false) : Cell(){
_index = c1->getIndex();
_tag = c1->getTag();
_dim = c1->getDim();
_inSubdomain = c1->inSubdomain();
......@@ -518,7 +545,9 @@ class CellComplex
// one for each dimension
std::set<Cell*, Less_Cell> _cells[4];
std::vector<Cell*> _trash;
// trashbin for cell pointers removed from cell complex
std::set<Cell*, Less_Cell> _cells2[4];
std::set<Cell*, Less_Cell> _trash;
//std::set<Cell*, Less_Cell> _originalCells[4];
......@@ -531,7 +560,7 @@ class CellComplex
// iterator for the cells of same dimension
typedef std::set<Cell*, Less_Cell>::iterator citer;
protected:
private:
// enqueue cells in queue if they are not there already
void enqueueCells(std::list<Cell*>& cells,
std::queue<Cell*>& Q, std::set<Cell*, Less_Cell>& Qset);
......@@ -542,6 +571,19 @@ class CellComplex
//virtual void insert_cells(bool subdomain, bool boundary);
void insert_cells(bool subdomain, bool boundary);
// remove a cell from this cell complex
void removeCell(Cell* cell);
void replaceCells(Cell* c1, Cell* c2, Cell* newCell, bool orMatch, bool co=false);
void insertCell(Cell* cell);
// reduction of this cell complex
// removes reduction pairs of cell of dimension dim and dim-1
int reduction(int dim);
// queued coreduction presented in Mrozek's paper
int coreduction(Cell* generator);
public:
CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> subdomain, std::set<Cell*, Less_Cell> cells ) {
_domain = domain;
......@@ -578,12 +620,21 @@ class CellComplex
CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> subdomain );
CellComplex(){}
~CellComplex(){
for(int i = 0; i < _trash.size(); i++){
Cell* cell = _trash.at(i);
//for(std::set<Cell*, Less_Cell>::iterator it = _trash.begin(); it != _trash.end(); it++){
// Cell* cell = *it;
// delete cell;
//}
for(int i = 0; i < 4; i++){
_cells[i].clear();
for(citer cit = _cells2[i].begin(); cit != _cells2[i].end(); cit++){
Cell* cell = *cit;
delete cell;
}
}
}
// get the number of certain dimensional cells
int getSize(int dim){ return _cells[dim].size(); }
......@@ -600,29 +651,16 @@ class CellComplex
// implementation will vary depending on cell type
inline int kappa(Cell* sigma, Cell* tau) const { return sigma->kappa(tau); }
// remove a cell from this cell complex
void removeCell(Cell* cell);
void replaceCells(Cell* c1, Cell* c2, Cell* newCell, bool orMatch, bool co=false);
void insertCell(Cell* cell);
// 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()) ); }
// reduction of this cell complex
// removes reduction pairs of cell of dimension dim and dim-1
int reduction(int dim);
// useful functions for (co)reduction of cell complex
int reduceComplex(bool omitHighdim = false);
int coreduceComplex(bool omitlowDim = false);
// queued coreduction presented in Mrozek's paper
int coreduction(Cell* generator);
// add every volume, face and edge its missing boundary cells
// void repairComplex(int i=3);
......@@ -640,6 +678,13 @@ class CellComplex
int combine(int dim);
int cocombine(int dim);
void emptyTrash() {
for(std::set<Cell*, Less_Cell>::iterator it = _trash.begin(); it != _trash.end(); it++){
Cell* cell = *it;
delete cell;
}
}
void computeBettiNumbers();
int getBettiNumber(int i) { if(i > -1 && i < 4) return _betti[i]; else return 0; }
......@@ -659,6 +704,27 @@ class CellComplex
}
}
void checkCoherence(){
for(int i = 0; i < 4; i++){
for(citer cit = firstCell(i); cit != lastCell(i); cit++){
Cell* cell = *cit;
std::list< std::pair<int, Cell*> > boundary;
for(std::list< std::pair<int, Cell* > >::iterator it = boundary.begin(); it != boundary.end(); it++){
Cell* bdCell = (*it).second;
citer cit = _cells[bdCell->getDim()].find(bdCell);
if(cit == lastCell(bdCell->getDim())) printf("Warning! Incoherent boundary!");
}
std::list< std::pair<int, Cell*> > coboundary;
for(std::list< std::pair<int, Cell* > >::iterator it = coboundary.begin(); it != coboundary.end(); it++){
Cell* cbdCell = (*it).second;
citer cit = _cells[cbdCell->getDim()].find(cbdCell);
if(cit == lastCell(cbdCell->getDim())) printf("Warning! Incoherent coboundary!");
}
}
}
}
};
#endif
......
......@@ -76,15 +76,25 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){
Cell* bdCell = (*it).second;
if(!bdCell->inSubdomain()){
int old_elem = 0;
//printf("kaka1: %d, kaka2: %d \n", bdCell->getIndex(), cell->getIndex());
//printf("cell1: %d, cell2: %d \n", bdCell->getIndex(), cell->getIndex());
if(bdCell->getIndex() > gmp_matrix_rows( _HMatrix[dim]) || bdCell->getIndex() < 1
|| cell->getIndex() > gmp_matrix_cols( _HMatrix[dim]) || cell->getIndex() < 1){
printf("Warning: Index out of bound! HMatrix: %d. \n", dim);
}
else{
gmp_matrix_get_elem(elem, bdCell->getIndex(), cell->getIndex(), _HMatrix[dim]);
old_elem = mpz_get_si(elem);
mpz_set_si(elem, old_elem + (*it).first);
if( (old_elem + (*it).first) > 1 || (old_elem + (*it).first) < -1 ){
printf("Warning: Invalid incidence index: %d! HMatrix: %d. \n", (old_elem + (*it).first), dim);
mpz_set_si(elem, 0);
}
gmp_matrix_set_elem(elem, bdCell->getIndex(), cell->getIndex(), _HMatrix[dim]);
}
}
}
}
}
mpz_clear(elem);
......
......@@ -77,7 +77,16 @@ class ChainComplex{
}
_dim = 0;
}
~ChainComplex(){}
~ChainComplex(){
for(int i = 0; i < 5; i++){
destroy_gmp_matrix(_HMatrix[i]);
destroy_gmp_matrix(_kerH[i]);
destroy_gmp_matrix(_codH[i]);
destroy_gmp_matrix(_JMatrix[i]);
destroy_gmp_matrix(_QMatrix[i]);
destroy_gmp_matrix(_Hbasis[i]);
}
}
int getDim() { return _dim; }
......
......@@ -14,6 +14,7 @@
Homology::Homology(GModel* model, std::vector<int> physicalDomain, std::vector<int> physicalSubdomain){
_model = model;
_combine = true;
Msg::Info("Creating a Cell Complex...");
double t1 = Cpu();
......@@ -74,14 +75,22 @@ void Homology::findGenerators(std::string fileName){
Msg::Info("Reducing Cell Complex...");
double t1 = Cpu();
int omitted = _cellComplex->reduceComplex(true);
//_cellComplex->checkCoherence();
if(getCombine()){
_cellComplex->combine(3);
_cellComplex->reduceComplex(false);
_cellComplex->combine(2);
_cellComplex->reduceComplex(false);
_cellComplex->combine(1);
_cellComplex->reduceComplex(false);
}
//_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->writeComplexMSH(fileName);
Msg::Info("Computing homology groups...");
......@@ -119,29 +128,37 @@ void Homology::findGenerators(std::string fileName){
Msg::Info("H3 = %d", HRank[3]);
if(omitted != 0) Msg::Info("Computation of %dD generators was omitted.", _cellComplex->getDim());
Msg::Info("Wrote results to %s.", fileName.c_str());
delete chains;
printf("H0 = %d \n", HRank[0]);
printf("H1 = %d \n", HRank[1]);
printf("H2 = %d \n", HRank[2]);
printf("H3 = %d \n", HRank[3]);
return;
}
void Homology::findThickCuts(std::string fileName){
_cellComplex->removeSubdomain();
Msg::Info("Reducing Cell Complex...");
double t1 = Cpu();
int omitted = _cellComplex->coreduceComplex(true);
//_cellComplex->checkCoherence();
if(getCombine()){
_cellComplex->cocombine(0);
_cellComplex->cocombine(1);
_cellComplex->cocombine(2);
}
//_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->writeComplexMSH(fileName);
Msg::Info("Computing homology groups...");
......@@ -154,6 +171,8 @@ void Homology::findThickCuts(std::string fileName){
int dim = _cellComplex->getDim();
int HRank[4];
for(int i = 0; i < 4; i++) HRank[i] = 0;
for(int j = 3; j > -1; j--){
......@@ -183,10 +202,15 @@ void Homology::findThickCuts(std::string fileName){
Msg::Info("H3 = %d", HRank[3]);
if(omitted != 0) Msg::Info("Computation of %dD thick cuts was omitted.", dim);
Msg::Info("Wrote results to %s.", fileName.c_str());
delete chains;
printf("H0 = %d \n", HRank[0]);
printf("H1 = %d \n", HRank[1]);
printf("H2 = %d \n", HRank[2]);
printf("H3 = %d \n", HRank[3]);
return;
}
#endif
......@@ -28,11 +28,16 @@ class Homology
CellComplex* _cellComplex;
GModel* _model;
bool _combine;
public:
Homology(GModel* model, std::vector<int> physicalDomain, std::vector<int> physicalSubdomain);
~Homology(){ delete _cellComplex; }
bool getCombine() { return _combine; }
bool setCombine(bool combine) { _combine = combine; }
void findGenerators(std::string fileName);
void findThickCuts(std::string fileName);
......
......@@ -15,11 +15,14 @@
#if defined(HAVE_KBIPACK)
StringXNumber HomologyComputationOptions_Number[] = {
{GMSH_FULLRC, "Physical group for domain", NULL, 0.},
{GMSH_FULLRC, "Physical group for subdomain", NULL, 0.},
{GMSH_FULLRC, "1. Physical group for domain", NULL, 0.},
{GMSH_FULLRC, "2. Physical group for domain", NULL, 0.},
{GMSH_FULLRC, "1. Physical group for subdomain", NULL, 0.},
{GMSH_FULLRC, "2. Physical group for subdomain", NULL, 0.},
{GMSH_FULLRC, "Compute generators", NULL, 1.},
{GMSH_FULLRC, "Compute thick cuts", NULL, 0.},
{GMSH_FULLRC, "Swap subdomain", NULL, 0.},
{GMSH_FULLRC, "Combine cells", NULL, 1.}
};
StringXString HomologyComputationOptions_String[] = {
......@@ -89,22 +92,27 @@ PView *GMSH_HomologyComputationPlugin::execute(PView *v)
std::vector<int> subdomain;
domain.push_back((int)HomologyComputationOptions_Number[0].def);
subdomain.push_back((int)HomologyComputationOptions_Number[1].def);
domain.push_back((int)HomologyComputationOptions_Number[1].def);
subdomain.push_back((int)HomologyComputationOptions_Number[2].def);
subdomain.push_back((int)HomologyComputationOptions_Number[3].def);
int gen = (int)HomologyComputationOptions_Number[2].def;
int cuts = (int)HomologyComputationOptions_Number[3].def;
int swap = (int)HomologyComputationOptions_Number[4].def;
int gens = (int)HomologyComputationOptions_Number[4].def;
int cuts = (int)HomologyComputationOptions_Number[5].def;
int swap = (int)HomologyComputationOptions_Number[6].def;
int combine = (int)HomologyComputationOptions_Number[7].def;
GModel *m = GModel::current();
Homology* homology = new Homology(m, domain, subdomain);
if(combine == 0) homology->setCombine(false);
if(swap == 1) homology->swapSubdomain();
if(gen == 1 && cuts != 1) {
if(gens == 1 && cuts != 1) {
homology->findGenerators(fileName);
GmshMergeFile(fileName);
}
else if(cuts == 1 && gen != 1) {
else if(cuts == 1 && gens != 1) {
homology->findThickCuts(fileName);
GmshMergeFile(fileName);
}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment