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

Simplified omitting.

parent d71514e7
Branches
Tags
No related merge requests found
Show whitespace changes
Inline Side-by-side
Geo/Cell.cpp
+ 18
− 4
View file @ 364bf6fc
...@@ -43,10 +43,10 @@ Cell::~Cell() ...@@ -43,10 +43,10 @@ Cell::~Cell()
if(_delimage) delete _image; if(_delimage) delete _image;
} }
bool Cell::findBoundaryCells(std::vector<Cell*>& bdCells) bool Cell::findBoundaryElements(std::vector<MElement*>& bdElements)
{ {
if(_combined) return false; if(_combined) return false;
bdCells.clear(); bdElements.clear();
MElementFactory factory; MElementFactory factory;
for(int i = 0; i < getNumFacets(); i++){ for(int i = 0; i < getNumFacets(); i++){
std::vector<MVertex*> vertices; std::vector<MVertex*> vertices;
...@@ -69,8 +69,7 @@ bool Cell::findBoundaryCells(std::vector<Cell*>& bdCells) ...@@ -69,8 +69,7 @@ bool Cell::findBoundaryCells(std::vector<Cell*>& bdCells)
} }
MElement* element = factory.create(newtype, vertices, 0, MElement* element = factory.create(newtype, vertices, 0,
_image->getPartition()); _image->getPartition());
Cell* cell = new Cell(element); bdElements.push_back(element);
bdCells.push_back(cell);
} }
return true; return true;
} }
...@@ -347,4 +346,19 @@ CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell() ...@@ -347,4 +346,19 @@ CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell()
} }
CombinedCell::CombinedCell(std::vector<Cell*>& cells) : Cell()
{
_num = ++_globalNum;
_index = cells.at(0)->getIndex();
_subdomain = cells.at(0)->inSubdomain();
_combined = true;
// cells
for(unsigned int i = 0; i < cells.size(); i++){
Cell* c = cells.at(i);
_cells[c] = 1;
}
}
#endif #endif
Geo/Cell.h
+ 2
− 1
View file @ 364bf6fc
...@@ -100,7 +100,7 @@ class Cell ...@@ -100,7 +100,7 @@ class Cell
else return false; } else return false; }
// find the cells on the boundary of this cell // find the cells on the boundary of this cell
bool findBoundaryCells(std::vector<Cell*>& bdCells); bool findBoundaryElements(std::vector<MElement*>& bdElements);
// get boundary cell orientation // get boundary cell orientation
int findBoundaryCellOrientation(Cell* cell); int findBoundaryCellOrientation(Cell* cell);
...@@ -208,6 +208,7 @@ class CombinedCell : public Cell{ ...@@ -208,6 +208,7 @@ class CombinedCell : public Cell{
public: public:
CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co=false); CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co=false);
CombinedCell(std::vector<Cell*>& cells);
~CombinedCell() {} ~CombinedCell() {}
int getDim() const { return _cells.begin()->first->getDim(); } int getDim() const { return _cells.begin()->first->getDim(); }
......
Geo/CellComplex.cpp
+ 57
− 51
View file @ 364bf6fc
...@@ -28,15 +28,12 @@ CellComplex::CellComplex( std::vector<MElement*>& domainElements, ...@@ -28,15 +28,12 @@ CellComplex::CellComplex( std::vector<MElement*>& domainElements,
void CellComplex::panic_exit(){ void CellComplex::panic_exit(){
for(int i = 0; i < 4; i++){ for(int i = 0; i < 4; i++){
_cells[i].clear();
for(citer cit = _ocells[i].begin(); cit != _ocells[i].end(); cit++){ for(citer cit = _ocells[i].begin(); cit != _ocells[i].end(); cit++){
Cell* cell = *cit; Cell* cell = *cit;
delete cell; delete cell;
} }
_ocells[i].clear();
} }
for(unsigned int i = 0; i < _newcells.size(); i++) delete _newcells.at(i); for(unsigned int i = 0; i < _newcells.size(); i++) delete _newcells.at(i);
_newcells.clear();
printf("ERROR: No proper cell complex could be constructed!\n"); printf("ERROR: No proper cell complex could be constructed!\n");
} }
...@@ -57,26 +54,30 @@ bool CellComplex::insert_cells(std::vector<MElement*>& elements, ...@@ -57,26 +54,30 @@ bool CellComplex::insert_cells(std::vector<MElement*>& elements,
printf("ERROR: mesh element %d not implemented! \n", type); printf("ERROR: mesh element %d not implemented! \n", type);
return false; return false;
} }
//Cell* cell = _cellPool.construct(element);
Cell* cell = new Cell(element); Cell* cell = new Cell(element);
cell->setInSubdomain(subdomain); cell->setInSubdomain(subdomain);
std::pair<citer, bool> insertInfo = _cells[cell->getDim()].insert(cell); std::pair<citer, bool> insertInfo = _cells[cell->getDim()].insert(cell);
if(!insertInfo.second) delete cell; if(!insertInfo.second) delete cell;
//if(!insertInfo.second) _cellPool.free(cell);
} }
// add lower dimensional cells recursively // add lower dimensional cells recursively
for (int dim = 3; dim > 0; dim--){ for (int dim = 3; dim > 0; dim--){
for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){ for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
Cell* cell = *cit; Cell* cell = *cit;
std::vector<Cell*> bdCells; std::vector<MElement*> bdElements;
if(!cell->findBoundaryCells(bdCells)) return false; if(!cell->findBoundaryElements(bdElements)) return false;
for(unsigned int i = 0; i < bdCells.size(); i++){ for(unsigned int i = 0; i < bdElements.size(); i++){
Cell* newCell = bdCells.at(i); //Cell* newCell = _cellPool.construct(bdElements.at(i));
Cell* newCell = new Cell(bdElements.at(i));
newCell->setInSubdomain(subdomain); newCell->setInSubdomain(subdomain);
newCell->setDeleteImage(true); newCell->setDeleteImage(true);
std::pair<citer, bool> insertInfo = std::pair<citer, bool> insertInfo =
_cells[newCell->getDim()].insert(newCell); _cells[newCell->getDim()].insert(newCell);
if(!insertInfo.second){ // the cell was already in the cell complex if(!insertInfo.second){ // the cell was already in the cell complex
delete newCell; delete newCell;
//_cellPool.free(newCell);
newCell = *(insertInfo.first); newCell = *(insertInfo.first);
} }
if(!subdomain) { if(!subdomain) {
...@@ -92,17 +93,12 @@ bool CellComplex::insert_cells(std::vector<MElement*>& elements, ...@@ -92,17 +93,12 @@ bool CellComplex::insert_cells(std::vector<MElement*>& elements,
CellComplex::~CellComplex() CellComplex::~CellComplex()
{ {
for(int i = 0; i < 4; i++){ for(int i = 0; i < 4; i++){
_cells[i].clear();
for(citer cit = _ocells[i].begin(); cit != _ocells[i].end(); cit++){ for(citer cit = _ocells[i].begin(); cit != _ocells[i].end(); cit++){
Cell* cell = *cit; Cell* cell = *cit;
delete cell; delete cell;
} }
_ocells[i].clear();
} }
for(unsigned int i = 0; i < _newcells.size(); i++) delete _newcells.at(i); for(unsigned int i = 0; i < _newcells.size(); i++) delete _newcells.at(i);
_newcells.clear();
} }
void CellComplex::insertCell(Cell* cell) void CellComplex::insertCell(Cell* cell)
...@@ -119,6 +115,7 @@ void CellComplex::insertCell(Cell* cell) ...@@ -119,6 +115,7 @@ void CellComplex::insertCell(Cell* cell)
void CellComplex::removeCell(Cell* cell, bool other) void CellComplex::removeCell(Cell* cell, bool other)
{ {
if(!hasCell(cell)) return;
std::map<Cell*, int, Less_Cell > coboundary; std::map<Cell*, int, Less_Cell > coboundary;
cell->getCoboundary(coboundary); cell->getCoboundary(coboundary);
std::map<Cell*, int, Less_Cell > boundary; std::map<Cell*, int, Less_Cell > boundary;
...@@ -158,7 +155,8 @@ void CellComplex::enqueueCells(std::map<Cell*, int, Less_Cell>& cells, ...@@ -158,7 +155,8 @@ void CellComplex::enqueueCells(std::map<Cell*, int, Less_Cell>& cells,
} }
} }
int CellComplex::coreduction(Cell* startCell, int omitted) int CellComplex::coreduction(Cell* startCell, bool omit,
std::vector<Cell*>& omittedCells)
{ {
int coreductions = 0; int coreductions = 0;
...@@ -184,8 +182,8 @@ int CellComplex::coreduction(Cell* startCell, int omitted) ...@@ -184,8 +182,8 @@ int CellComplex::coreduction(Cell* startCell, int omitted)
bd_s.begin()->first->getCoboundary(cbd_c); bd_s.begin()->first->getCoboundary(cbd_c);
enqueueCells(cbd_c, Q, Qset); enqueueCells(cbd_c, Q, Qset);
removeCell(bd_s.begin()->first); removeCell(bd_s.begin()->first);
if(bd_s.begin()->first->getDim() == 0 && omitted > 0){ if(bd_s.begin()->first->getDim() == 0 && omit){
_store.at(omitted-1).insert(bd_s.begin()->first); omittedCells.push_back(bd_s.begin()->first);
} }
coreductions++; coreductions++;
} }
...@@ -197,7 +195,8 @@ int CellComplex::coreduction(Cell* startCell, int omitted) ...@@ -197,7 +195,8 @@ int CellComplex::coreduction(Cell* startCell, int omitted)
return coreductions; return coreductions;
} }
int CellComplex::reduction(int dim, int omitted) int CellComplex::reduction(int dim, bool omit,
std::vector<Cell*>& omittedCells)
{ {
if(dim < 1 || dim > 3) return 0; if(dim < 1 || dim > 3) return 0;
...@@ -213,9 +212,9 @@ int CellComplex::reduction(int dim, int omitted) ...@@ -213,9 +212,9 @@ int CellComplex::reduction(int dim, int omitted)
Cell* cell = *cit; Cell* cell = *cit;
if( cell->getCoboundarySize() == 1 if( cell->getCoboundarySize() == 1
&& inSameDomain(cell, cell->firstCoboundary()->first)){ && inSameDomain(cell, cell->firstCoboundary()->first)){
++cit; cit++;
if(dim == getDim() && omitted > 0){ if(dim == getDim() && omit){
_store.at(omitted-1).insert(cell->firstCoboundary()->first); omittedCells.push_back(cell->firstCoboundary()->first);
} }
removeCell(cell->firstCoboundary()->first); removeCell(cell->firstCoboundary()->first);
removeCell(cell); removeCell(cell);
...@@ -230,7 +229,8 @@ int CellComplex::reduction(int dim, int omitted) ...@@ -230,7 +229,8 @@ int CellComplex::reduction(int dim, int omitted)
return count; return count;
} }
int CellComplex::coreduction(int dim, int omitted) int CellComplex::coreduction(int dim, bool omit,
std::vector<Cell*>& omittedCells)
{ {
if(dim < 1 || dim > 3) return 0; if(dim < 1 || dim > 3) return 0;
...@@ -247,8 +247,8 @@ int CellComplex::coreduction(int dim, int omitted) ...@@ -247,8 +247,8 @@ int CellComplex::coreduction(int dim, int omitted)
if( cell->getBoundarySize() == 1 if( cell->getBoundarySize() == 1
&& inSameDomain(cell, cell->firstBoundary()->first)){ && inSameDomain(cell, cell->firstBoundary()->first)){
++cit; ++cit;
if(dim-1 == 0 && omitted > 0){ if(dim-1 == 0 && omit){
_store.at(omitted-1).insert(cell->firstBoundary()->first); omittedCells.push_back(cell->firstBoundary()->first);
} }
removeCell(cell->firstBoundary()->first); removeCell(cell->firstBoundary()->first);
removeCell(cell); removeCell(cell);
...@@ -269,25 +269,29 @@ int CellComplex::reduceComplex(bool omit) ...@@ -269,25 +269,29 @@ int CellComplex::reduceComplex(bool omit)
getSize(3), getSize(2), getSize(1), getSize(0)); getSize(3), getSize(2), getSize(1), getSize(0));
int count = 0; int count = 0;
for(int i = 3; i > 0; i--) count = count + reduction(i); std::vector<Cell*> empty;
for(int i = 3; i > 0; i--) count = count + reduction(i, false, empty);
if(omit){ if(omit){
int omitted = 0;
_store.clear();
removeSubdomain(); removeSubdomain();
std::vector<Cell*> newCells;
while (getSize(getDim()) != 0){ while (getSize(getDim()) != 0){
citer cit = firstCell(getDim()); citer cit = firstCell(getDim());
Cell* cell = *cit; Cell* cell = *cit;
removeCell(cell, false); removeCell(cell, false);
std::vector<Cell*> omittedCells;
omittedCells.push_back(cell);
omitted++; for(int j = 3; j > 0; j--) reduction(j, true, omittedCells);
std::set< Cell*, Less_Cell > omittedCells; CombinedCell* newcell = new CombinedCell(omittedCells);
_store.push_back(omittedCells); newCells.push_back(newcell);
_store.at(omitted-1).insert(cell); }
for(int j = 3; j > 0; j--) reduction(j, omitted); for(unsigned int i = 0; i < newCells.size(); i++){
insertCell(newCells.at(i));
} }
} }
...@@ -295,9 +299,9 @@ int CellComplex::reduceComplex(bool omit) ...@@ -295,9 +299,9 @@ int CellComplex::reduceComplex(bool omit)
getSize(3), getSize(2), getSize(1), getSize(0)); getSize(3), getSize(2), getSize(1), getSize(0));
combine(3); combine(3);
reduction(2); reduction(2, false, empty);
combine(2); combine(2);
reduction(1); reduction(1, false, empty);
combine(1); combine(1);
printf(" %d volumes, %d faces, %d edges and %d vertices. \n", printf(" %d volumes, %d faces, %d edges and %d vertices. \n",
...@@ -325,31 +329,33 @@ int CellComplex::coreduceComplex(bool omit) ...@@ -325,31 +329,33 @@ int CellComplex::coreduceComplex(bool omit)
int count = 0; int count = 0;
removeSubdomain(); removeSubdomain();
std::vector<Cell*> empty;
for(int dim = 0; dim < 4; dim++){ for(int dim = 0; dim < 4; dim++){
citer cit = firstCell(dim); citer cit = firstCell(dim);
while(cit != lastCell(dim)){ while(cit != lastCell(dim)){
Cell* cell = *cit; Cell* cell = *cit;
count = count + coreduction(cell); count = count + coreduction(cell, false, empty);
if(count != 0) break; if(count != 0) break;
cit++; cit++;
} }
} }
int omitted = 0;
if(omit){ if(omit){
_store.clear(); std::vector<Cell*> newCells;
while (getSize(0) != 0){ while (getSize(0) != 0){
citer cit = firstCell(0); citer cit = firstCell(0);
Cell* cell = *cit; Cell* cell = *cit;
removeCell(cell, false); removeCell(cell, false);
std::set< Cell*, Less_Cell > omittedCells; std::vector<Cell*> omittedCells;
omitted++; omittedCells.push_back(cell);
_store.push_back(omittedCells);
_store.at(omitted-1).insert(cell); coreduction(cell, true, omittedCells);
coreduction(cell, omitted); CombinedCell* newcell = new CombinedCell(omittedCells);
newCells.push_back(newcell);
}
for(unsigned int i = 0; i < newCells.size(); i++){
insertCell(newCells.at(i));
} }
} }
...@@ -357,17 +363,17 @@ int CellComplex::coreduceComplex(bool omit) ...@@ -357,17 +363,17 @@ int CellComplex::coreduceComplex(bool omit)
getSize(3), getSize(2), getSize(1), getSize(0)); getSize(3), getSize(2), getSize(1), getSize(0));
cocombine(0); cocombine(0);
coreduction(1); coreduction(1, false, empty);
cocombine(1); cocombine(1);
coreduction(2); coreduction(2, false, empty);
cocombine(2); cocombine(2);
coreduction(3); coreduction(3, false, empty);
coherent(); coherent();
printf(" %d volumes, %d faces, %d edges and %d vertices. \n", printf(" %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 0;
} }
int CellComplex::cocombine(int dim) int CellComplex::cocombine(int dim)
...@@ -465,6 +471,7 @@ int CellComplex::combine(int dim) ...@@ -465,6 +471,7 @@ int CellComplex::combine(int dim)
enqueueCells(bd_c, Q, Qset); enqueueCells(bd_c, Q, Qset);
CombinedCell* newCell = new CombinedCell(c1, c2, (or1 != or2)); CombinedCell* newCell = new CombinedCell(c1, c2, (or1 != or2));
//CombinedCell* newCell = _ccellPool.construct(c1, c2, (or1 != or2));
removeCell(c1); removeCell(c1);
removeCell(c2); removeCell(c2);
insertCell(newCell); insertCell(newCell);
...@@ -568,7 +575,6 @@ void CellComplex::restoreComplex() ...@@ -568,7 +575,6 @@ void CellComplex::restoreComplex()
delete cell; delete cell;
} }
_newcells.clear(); _newcells.clear();
_store.clear();
} }
void CellComplex::printComplex(int dim) void CellComplex::printComplex(int dim)
......
Geo/CellComplex.h
+ 8
− 13
View file @ 364bf6fc
...@@ -21,6 +21,7 @@ ...@@ -21,6 +21,7 @@
#include "MElement.h" #include "MElement.h"
#include "ChainComplex.h" #include "ChainComplex.h"
//#include "GModel.h" //#include "GModel.h"
#include <boost/pool/object_pool.hpp>
class Cell; class Cell;
...@@ -32,9 +33,6 @@ class CellComplex ...@@ -32,9 +33,6 @@ class CellComplex
// one for each dimension // one for each dimension
std::set<Cell*, Less_Cell> _cells[4]; std::set<Cell*, Less_Cell> _cells[4];
// temporary store for omitted cells (generators of the highest dimension)
std::vector< std::set<Cell*, Less_Cell> > _store;
// original cells of this cell complex // original cells of this cell complex
std::set<Cell*, Less_Cell> _ocells[4]; std::set<Cell*, Less_Cell> _ocells[4];
...@@ -60,8 +58,9 @@ class CellComplex ...@@ -60,8 +58,9 @@ class CellComplex
void removeCell(Cell* cell, bool other=true); void removeCell(Cell* cell, bool other=true);
void insertCell(Cell* cell); void insertCell(Cell* cell);
// queued coreduction presented in Mrozek's paper // queued coreduction
int coreduction(Cell* startCell, int omitted=0); int coreduction(Cell* startCell, bool omit,
std::vector<Cell*>& omittedCells);
public: public:
...@@ -106,8 +105,8 @@ class CellComplex ...@@ -106,8 +105,8 @@ class CellComplex
// (co)reduction of this cell complex // (co)reduction of this cell complex
// removes (co)reduction pairs of cell of dimension dim and dim-1 // removes (co)reduction pairs of cell of dimension dim and dim-1
int reduction(int dim, int omitted=0); int reduction(int dim, bool omit, std::vector<Cell*>& omittedCells);
int coreduction(int dim, int omitted=0); int coreduction(int dim, bool omit, std::vector<Cell*>& omittedCells);
// Cell combining for reduction and coreduction // Cell combining for reduction and coreduction
int combine(int dim); int combine(int dim);
...@@ -129,10 +128,6 @@ class CellComplex ...@@ -129,10 +128,6 @@ class CellComplex
void printEuler(){ void printEuler(){
printf("Euler characteristic: %d. \n", eulerCharacteristic()); } printf("Euler characteristic: %d. \n", eulerCharacteristic()); }
// get cells omitted by (co)reduction
int getNumOmitted() { return _store.size(); }
std::set<Cell*, Less_Cell> getOmitted(int i) { return _store.at(i); }
// restore the cell complex to its original state before (co)reduction // restore the cell complex to its original state before (co)reduction
void restoreComplex(); void restoreComplex();
......
Geo/ChainComplex.cpp
+ 1
− 2
View file @ 364bf6fc
...@@ -96,7 +96,6 @@ ChainComplex::ChainComplex(CellComplex* cellComplex, int domain) ...@@ -96,7 +96,6 @@ ChainComplex::ChainComplex(CellComplex* cellComplex, int domain)
} }
mpz_clear(elem); mpz_clear(elem);
} }
_kerH[dim] = NULL; _kerH[dim] = NULL;
_codH[dim] = NULL; _codH[dim] = NULL;
_JMatrix[dim] = NULL; _JMatrix[dim] = NULL;
...@@ -482,7 +481,7 @@ void ChainComplex::getBasisChain(std::map<Cell*, int, Less_Cell>& chain, ...@@ -482,7 +481,7 @@ void ChainComplex::getBasisChain(std::map<Cell*, int, Less_Cell>& chain,
} }
mpz_clear(elem); mpz_clear(elem);
if(deform) smoothenChain(chain); if(deform && (dim == 1 || dim == 2) ) smoothenChain(chain);
} }
int ChainComplex::getBasisSize(int dim, int basis) int ChainComplex::getBasisSize(int dim, int basis)
......
Geo/Homology.cpp
+ 1
− 38
View file @ 364bf6fc
...@@ -168,23 +168,6 @@ void Homology::findGenerators(CellComplex* cellComplex) ...@@ -168,23 +168,6 @@ void Homology::findGenerators(CellComplex* cellComplex)
_generators.push_back(chain->createPGroup()); _generators.push_back(chain->createPGroup());
delete chain; delete chain;
} }
if(j == cellComplex->getDim() && cellComplex->getNumOmitted() > 0){
for(int i = 0; i < cellComplex->getNumOmitted(); i++){
std::string generator;
convert(i+1, generator);
std::string name = "H" + dimension + domainString + generator;
std::vector<int> coeffs (cellComplex->getOmitted(i).size(),1);
Chain* chain = new Chain(cellComplex->getOmitted(i), coeffs,
cellComplex, _model, name, 1);
if(chain->getSize() == 0){
delete chain;
continue;
}
HRank[j] = HRank[j] + 1;
_generators.push_back(chain->createPGroup());
delete chain;
}
}
} }
if(_fileName != "") writeGeneratorsMSH(); if(_fileName != "") writeGeneratorsMSH();
...@@ -267,26 +250,6 @@ void Homology::findDualGenerators(CellComplex* cellComplex) ...@@ -267,26 +250,6 @@ void Homology::findDualGenerators(CellComplex* cellComplex)
dim-j, i, chain->getTorsion()); dim-j, i, chain->getTorsion());
} }
} }
if(j == 0 && cellComplex->getNumOmitted() > 0){
for(int i = 0; i < cellComplex->getNumOmitted(); i++){
std::string generator;
convert(i+1, generator);
std::string name
= "H" + dimension + "*" +
getDomainString(_domain, _subdomain) + generator;
std::vector<int> coeffs (cellComplex->getOmitted(i).size(),1);
Chain* chain = new Chain(cellComplex->getOmitted(i), coeffs,
cellComplex, _model, name, 1);
if(chain->getSize() == 0) {
delete chain;
continue;
}
HRank[dim-j] = HRank[dim-j] + 1;
_generators.push_back(chain->createPGroup());
delete chain;
}
}
} }
if(_fileName != "") writeGeneratorsMSH(); if(_fileName != "") writeGeneratorsMSH();
...@@ -313,7 +276,7 @@ void Homology::findHomSequence(){ ...@@ -313,7 +276,7 @@ void Homology::findHomSequence(){
Msg::StatusBar(1, false, "Reducing..."); Msg::StatusBar(1, false, "Reducing...");
double t1 = Cpu(); double t1 = Cpu();
cellComplex->reduceComplex(false); cellComplex->reduceComplex();
double t2 = Cpu(); double t2 = Cpu();
Msg::Info("Cell Complex reduction complete (%g s).", t2 - t1); Msg::Info("Cell Complex reduction complete (%g s).", t2 - t1);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment