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

Greatly improved computational complexity of reduction and coreduction, 

about O(n^2) -> O(n*log(n)).

The bottleneck is now in the Smith normal form computation of the reduced
cell complexes, O(n^3).
parent 4280a4e4
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Geo/CellComplex.cpp
+ 320
− 111
View file @ c9f05afb
......@@ -11,6 +11,8 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su
_domain = domain;
_subdomain = 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){
......@@ -40,16 +42,13 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su
GEntity* entity = *itb;
if(edge->tag() == entity->tag()){
_boundary.erase(itb);
//erase = itb;
duplicate = true;
break;
}
}
//if(duplicate) _boundary.erase(erase);
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++){
......@@ -69,28 +68,62 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su
}
}
// insert cells into cell complex
// subdomain need to be inserted first!
insertCells(true, true);
insertCells(false, true);
insertCells(false, false);
insert_cells2(true, true);
insert_cells2(false, true);
insert_cells2(false, false);
// find mesh vertices in the domain
for(unsigned int j=0; j < domain.size(); j++) {
for(unsigned int i=0; i < domain.at(j)->getNumMeshElements(); i++){
for(int k=0; k < domain.at(j)->getMeshElement(i)->getNumVertices(); k++){
MVertex* vertex = domain.at(j)->getMeshElement(i)->getVertex(k);
_domainVertices.insert(vertex);
}
}
}
for(unsigned int j=0; j < subdomain.size(); j++) {
for(unsigned int i=0; i < subdomain.at(j)->getNumMeshElements(); i++){
for(int k=0; k < subdomain.at(j)->getMeshElement(i)->getNumVertices(); k++){
MVertex* vertex = subdomain.at(j)->getMeshElement(i)->getVertex(k);
_domainVertices.insert(vertex);
}
}
}
// attach induvivial tags to cells
int tag = 1;
for(int i = 0; i < 4; i++){
for(citer cit = firstCell(i); cit != lastCell(i); cit++){
Cell* cell = *cit;
cell->setTag(tag);
tag++;
// make sure that boundary and coboundary information is coherent
std::list<Cell*> boundary = cell->getBoundary();
boundary.sort(Less_Cell());
boundary.unique(Equal_Cell());
cell->setBoundary(boundary);
std::list<Cell*> coboundary = cell->getCoboundary();
coboundary.sort(Less_Cell());
coboundary.unique(Equal_Cell());
cell->setCoboundary(coboundary);
}
_originalCells[i] = _cells[i];
}
_simplicial = true;
}
void CellComplex::insertCells(bool subdomain, bool boundary){
void CellComplex::insert_cells(bool subdomain, bool boundary){
std::vector<GEntity*> domain;
......@@ -101,26 +134,34 @@ void CellComplex::insertCells(bool subdomain, bool boundary){
std::vector<int> vertices;
int vertex = 0;
std::pair<citer, bool> insertInfo;
for(unsigned int j=0; j < domain.size(); j++) {
for(unsigned int i=0; i < domain.at(j)->getNumMeshElements(); i++){
vertices.clear();
for(int k=0; k < domain.at(j)->getMeshElement(i)->getNumVertices(); k++){
MVertex* vertex = domain.at(j)->getMeshElement(i)->getVertex(k);
vertices.push_back(vertex->getNum());
//_cells[0].insert(new ZeroSimplex(vertex->getNum(), subdomain, vertex->x(), vertex->y(), vertex->z() )); // Add vertices
_cells[0].insert(new ZeroSimplex(vertex->getNum(), subdomain, boundary));
Cell* cell = new ZeroSimplex(vertex->getNum(), subdomain, boundary);
insertInfo = _cells[0].insert(cell);
if(!insertInfo.second) delete cell;
}
if(domain.at(j)->getMeshElement(i)->getDim() == 3){
_cells[3].insert(new ThreeSimplex(vertices, subdomain, boundary)); // Add volumes
Cell* cell = new ThreeSimplex(vertices, subdomain, boundary); // Add volumes
insertInfo = _cells[3].insert(cell);
if(!insertInfo.second) delete cell;
}
for(int k=0; k < domain.at(j)->getMeshElement(i)->getNumFaces(); k++){
vertices.clear();
for(int l=0; l < domain.at(j)->getMeshElement(i)->getFace(k).getNumVertices(); l++){
vertex = domain.at(j)->getMeshElement(i)->getFace(k).getVertex(l)->getNum();
vertices.push_back(vertex);
}
_cells[2].insert(new TwoSimplex(vertices, subdomain, boundary)); // Add faces
Cell* cell = new TwoSimplex(vertices, subdomain, boundary); // Add faces
insertInfo = _cells[2].insert(cell);
if(!insertInfo.second) delete cell;
}
for(int k=0; k < domain.at(j)->getMeshElement(i)->getNumEdges(); k++){
......@@ -129,7 +170,9 @@ void CellComplex::insertCells(bool subdomain, bool boundary){
vertex = domain.at(j)->getMeshElement(i)->getEdge(k).getVertex(l)->getNum();
vertices.push_back(vertex);
}
_cells[1].insert(new OneSimplex(vertices, subdomain, boundary)); // Add edges
Cell* cell = new OneSimplex(vertices, subdomain, boundary); // Add edges
insertInfo = _cells[1].insert(cell);
if(!insertInfo.second) delete cell;
}
}
......@@ -137,6 +180,71 @@ void CellComplex::insertCells(bool subdomain, bool boundary){
}
void CellComplex::insert_cells2(bool subdomain, bool boundary){
// works only for simplcial meshes at the moment!
std::vector<GEntity*> domain;
if(subdomain) domain = _subdomain;
else if(boundary) domain = _boundary;
else domain = _domain;
std::vector<int> vertices;
int vertex = 0;
std::pair<citer, bool> insertInfo;
// add highest dimensional cells
for(unsigned int j=0; j < domain.size(); j++) {
for(unsigned int i=0; i < domain.at(j)->getNumMeshElements(); i++){
vertices.clear();
for(int k=0; k < domain.at(j)->getMeshElement(i)->getNumVertices(); k++){
MVertex* vertex = domain.at(j)->getMeshElement(i)->getVertex(k);
vertices.push_back(vertex->getNum());
}
int dim = domain.at(j)->getMeshElement(i)->getDim();
int type = domain.at(j)->getMeshElement(i)->getTypeForMSH();
Cell* cell;
if(dim == 3) cell = new ThreeSimplex(vertices, subdomain, boundary);
else if(dim == 2) cell = new TwoSimplex(vertices, subdomain, boundary);
else if(dim == 1) cell = new OneSimplex(vertices, subdomain, boundary);
else cell = new ZeroSimplex(vertices.at(0), subdomain, boundary);
insertInfo = _cells[dim].insert(cell);
if(!insertInfo.second) delete cell;
}
}
// add lower dimensional cells recursively
for (int dim = 3; dim > 0; dim--){
for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
Cell* cell = *cit;
std::vector<int> vertices;
for(int i = 0; i < cell->getNumFacets(); i++){
cell->getFacetVertices(i, vertices);
Cell* newCell;
if(dim == 3) newCell = new TwoSimplex(vertices, subdomain, boundary);
else if(dim == 2) newCell = new OneSimplex(vertices, subdomain, boundary);
else if(dim == 1) newCell = new ZeroSimplex(vertices.at(0), subdomain, boundary);
insertInfo = _cells[dim-1].insert(newCell);
if(!insertInfo.second){
delete newCell;
Cell* oldCell = *(insertInfo.first);
oldCell->addCoboundaryCell(cell);
cell->addBoundaryCell(oldCell);
}
else{
cell->addBoundaryCell(newCell);
newCell->addCoboundaryCell(cell);
}
}
}
}
}
void CellComplex::insertCell(Cell* cell){
_cells[cell->getDim()].insert(cell);
}
......@@ -157,25 +265,32 @@ int Simplex::kappa(Cell* tau) const{
return value;
}
std::set<Cell*, Less_Cell>::iterator CellComplex::findCell(int dim, std::vector<int>& vertices, bool original){
if(!original){
if(dim == 0) return _cells[dim].find(new ZeroSimplex(vertices.at(0)));
if(dim == 1) return _cells[dim].find(new OneSimplex(vertices));
if(dim == 2) return _cells[dim].find(new TwoSimplex(vertices));
return _cells[3].find(new ThreeSimplex(vertices));
}
else{
if(dim == 0) return _originalCells[dim].find(new ZeroSimplex(vertices.at(0)));
if(dim == 1) return _originalCells[dim].find(new OneSimplex(vertices));
if(dim == 2) return _originalCells[dim].find(new TwoSimplex(vertices));
return _originalCells[3].find(new ThreeSimplex(vertices));
}
Cell* cell;
if(dim == 0) cell = new ZeroSimplex(vertices.at(0));
else if(dim == 1) cell = new OneSimplex(vertices);
else if(dim == 2) cell = new TwoSimplex(vertices);
else cell = new ThreeSimplex(vertices);
citer cit;
if(!original) cit = _cells[dim].find(cell);
else cit = _originalCells[dim].find(cell);
delete cell;
return cit;
}
std::set<Cell*, Less_Cell>::iterator CellComplex::findCell(int dim, int vertex, int dummy){
if(dim == 0) return _cells[dim].lower_bound(new ZeroSimplex(vertex));
if(dim == 1) return _cells[dim].lower_bound(new OneSimplex(vertex, dummy));
if(dim == 2) return _cells[dim].lower_bound(new TwoSimplex(vertex, dummy));
return _cells[3].lower_bound(new ThreeSimplex(vertex, dummy));
Cell* cell;
if(dim == 0) cell = new ZeroSimplex(vertex);
else if(dim == 1) cell = new OneSimplex(vertex, dummy);
else if(dim == 2) cell = new TwoSimplex(vertex, dummy);
else cell = new ThreeSimplex(vertex, dummy);
citer cit = _cells[dim].lower_bound(cell);
delete cell;
return cit;
}
......@@ -184,7 +299,27 @@ std::vector<Cell*> CellComplex::bd(Cell* sigma){
int dim = sigma->getDim();
if(dim < 1) return boundary;
if(simplicial()){ // faster
std::vector<int> vertices = sigma->getVertexVector();
for(int j=0; j < vertices.size(); j++){
std::vector<int> bVertices;
// simplicial mesh assumed here!
for(int k=0; k < vertices.size(); k++){
if (k !=j ) bVertices.push_back(vertices.at(k));
}
citer cit2 = findCell(dim-1, bVertices);
if(cit2 != lastCell(dim-1)){
boundary.push_back(*cit2);
if(boundary.size() == sigma->getBdMaxSize()){
return boundary;
}
}
}
}
else{
citer start = findCell(dim-1, sigma->getVertex(0), 0);
if(start == lastCell(dim-1)) return boundary;
......@@ -193,13 +328,15 @@ std::vector<Cell*> CellComplex::bd(Cell* sigma){
//for(citer cit = firstCell(dim-1); cit != lastCell(dim-1); cit++){
for(citer cit = start; cit != end; cit++){
if(kappa(sigma, *cit) != 0){
//if(kappa(sigma, *cit) != 0){
if(sigma->kappa(*cit) != 0){
boundary.push_back(*cit);
if(boundary.size() == sigma->getBdMaxSize()){
return boundary;
}
}
}
}
return boundary;
}
std::vector< std::set<Cell*, Less_Cell>::iterator > CellComplex::bdIt(Cell* sigma){
......@@ -238,15 +375,19 @@ std::vector<Cell*> CellComplex::cbd(Cell* tau){
}
for(citer cit = firstCell(dim+1); cit != lastCell(dim+1); cit++){
if(kappa(*cit, tau) != 0){
//if(kappa(*cit, tau) != 0){
Cell* cell = *cit;
if(cell->kappa(tau) != 0){
coboundary.push_back(*cit);
if(coboundary.size() == tau->getCbdMaxSize()){
return coboundary;
}
}
}
return coboundary;
}
std::vector< std::set<Cell*, Less_Cell>::iterator > CellComplex::cbdIt(Cell* tau){
std::vector< std::set<Cell*, Less_Cell>::iterator > coboundary;
......@@ -268,20 +409,31 @@ std::vector< std::set<Cell*, Less_Cell>::iterator > CellComplex::cbdIt(Cell* tau
}
void CellComplex::removeCell(Cell* cell, bool deleteCell){
void CellComplex::removeCell(Cell* cell){
_cells[cell->getDim()].erase(cell);
//if(deleteCell) delete cell;
std::list<Cell*> coboundary = cell->getCoboundary();
std::list<Cell*> boundary = cell->getBoundary();
for(std::list<Cell*>::iterator it = coboundary.begin(); it != coboundary.end(); it++){
Cell* cbdCell = *it;
cbdCell->removeBoundaryCell(cell);
}
for(std::list<Cell*>::iterator it = boundary.begin(); it != boundary.end(); it++){
Cell* bdCell = *it;
bdCell->removeCoboundaryCell(cell);
}
}
void CellComplex::removeCellIt(std::set<Cell*, Less_Cell>::iterator cell){
Cell* c = *cell;
int dim = c->getDim();
_cells[dim].erase(cell);
//delete c;
}
void CellComplex::removeCellQset(Cell*& cell, std::set<Cell*, Less_Cell>& Qset){
Qset.erase(cell);
//delete cell;
}
void CellComplex::enqueueCellsIt(std::vector< std::set<Cell*, Less_Cell>::iterator >& cells,
......@@ -295,21 +447,23 @@ void CellComplex::enqueueCellsIt(std::vector< std::set<Cell*, Less_Cell>::iterat
}
}
}
void CellComplex::enqueueCells(std::vector<Cell*>& cells, std::queue<Cell*>& Q, std::set<Cell*, Less_Cell>& Qset){
for(unsigned int i=0; i < cells.size(); i++){
citer cit = Qset.find(cells.at(i));
if(cit == Qset.end()){
Qset.insert(cells.at(i));
Q.push(cells.at(i));
void CellComplex::enqueueCells(std::list<Cell*>& cells, std::queue<Cell*>& Q, std::set<Cell*, Less_Cell>& Qset){
for(std::list<Cell*>::iterator cit = cells.begin(); cit != cells.end(); cit++){
Cell* cell = *cit;
citer it = Qset.find(cell);
if(it == Qset.end()){
Qset.insert(cell);
Q.push(cell);
}
}
}
void CellComplex::enqueueCells(std::vector<Cell*>& cells, std::list<Cell*>& Q){
for(unsigned int i=0; i < cells.size(); i++){
std::list<Cell*>::iterator it = std::find(Q.begin(), Q.end(), cells.at(i));
void CellComplex::enqueueCells(std::list<Cell*>& cells, std::list<Cell*>& Q){
for(std::list<Cell*>::iterator cit = cells.begin(); cit != cells.end(); cit++){
Cell* cell = *cit;
std::list<Cell*>::iterator it = std::find(Q.begin(), Q.end(), cell);
if(it == Q.end()){
Q.push_back(cells.at(i));
Q.push_back(cell);
}
}
}
......@@ -327,27 +481,28 @@ int CellComplex::coreductionMrozek(Cell* generator){
Qset.insert(generator);
//removeCell(generator);
std::vector<Cell*> bd_s;
std::vector<Cell*> cbd_c;
std::list<Cell*> bd_s;
std::list<Cell*> cbd_c;
Cell* s;
int round = 0;
while( !Q.empty() ){
round++;
//printf("%d ", round);
s = Q.front();
Q.pop();
removeCellQset(s, Qset);
bd_s = bd(s);
if( bd_s.size() == 1 && inSameDomain(s, bd_s.at(0)) ){
bd_s = s->getBoundary();
if( bd_s.size() == 1 && inSameDomain(s, bd_s.front()) ){
removeCell(s);
cbd_c = cbd(bd_s.at(0));
cbd_c = bd_s.front()->getCoboundary();
enqueueCells(cbd_c, Q, Qset);
removeCell(bd_s.at(0));
removeCell(bd_s.front());
coreductions++;
}
else if(bd_s.empty()){
cbd_c = cbd(s);
cbd_c = s->getCoboundary();
enqueueCells(cbd_c, Q, Qset);
}
......@@ -356,6 +511,7 @@ int CellComplex::coreductionMrozek(Cell* generator){
//printf("Coreduction: %d loops with %d coreductions\n", round, coreductions);
return coreductions;
}
int CellComplex::coreductionMrozek2(Cell* generator){
int coreductions = 0;
......@@ -364,8 +520,8 @@ int CellComplex::coreductionMrozek2(Cell* generator){
Q.push_back(generator);
std::vector<Cell*> bd_s;
std::vector<Cell*> cbd_c;
std::list<Cell*> bd_s;
std::list<Cell*> cbd_c;
Cell* s;
int round = 0;
while( !Q.empty() ){
......@@ -373,16 +529,16 @@ int CellComplex::coreductionMrozek2(Cell* generator){
//printf("%d ", round);
s = Q.front();
Q.pop_front();
bd_s = bd(s);
if( bd_s.size() == 1 && inSameDomain(s, bd_s.at(0)) ){
bd_s = s->getBoundary();
if( bd_s.size() == 1 && inSameDomain(s, bd_s.front()) ){
removeCell(s);
cbd_c = cbd(bd_s.at(0));
cbd_c = bd_s.front()->getCoboundary();
enqueueCells(cbd_c, Q);
removeCell(bd_s.at(0));
removeCell(bd_s.front());
coreductions++;
}
else if(bd_s.empty()){
cbd_c = cbd(s);
cbd_c = s->getCoboundary();
enqueueCells(cbd_c, Q);
}
......@@ -392,6 +548,7 @@ int CellComplex::coreductionMrozek2(Cell* generator){
//printf("Coreduction: %d loops with %d coreductions\n", round, coreductions);
return coreductions;
}
int CellComplex::coreductionMrozek3(Cell* generator){
int coreductions = 0;
......@@ -430,7 +587,7 @@ int CellComplex::coreductionMrozek3(Cell* generator){
//printf("Coreduction: %d loops with %d coreductions\n", round, coreductions);
return coreductions;
}
int CellComplex::coreduction(int dim, bool deleteCells){
int CellComplex::coreduction(int dim){
if(dim < 0 || dim > 2) return 0;
......@@ -445,8 +602,8 @@ int CellComplex::coreduction(int dim, bool deleteCells){
bd_c = bd(cell);
if(bd_c.size() == 1 && inSameDomain(cell, bd_c.at(0)) ){
removeCell(cell, deleteCells);
removeCell(bd_c.at(0), deleteCells);
removeCell(cell);
removeCell(bd_c.at(0));
count++;
coreduced = true;
}
......@@ -475,8 +632,8 @@ int CellComplex::coreduction(int dim, std::set<Cell*, Less_Cell>& removedCells){
if(bd_c.size() == 1 && inSameDomain(cell, bd_c.at(0)) ){
removedCells.insert(cell);
removedCells.insert(bd_c.at(0));
removeCell(cell, false);
removeCell(bd_c.at(0), false);
removeCell(cell);
removeCell(bd_c.at(0));
count++;
coreduced = true;
}
......@@ -488,12 +645,40 @@ int CellComplex::coreduction(int dim, std::set<Cell*, Less_Cell>& removedCells){
}
int CellComplex::coreduction2(int dim){
if(dim < 0 || dim > 2) return 0;
std::list<Cell*> bd_c;
int count = 0;
bool coreduced = true;
while (coreduced){
coreduced = false;
for(citer cit = firstCell(dim+1); cit != lastCell(dim+1); cit++){
Cell* cell = *cit;
bd_c = cell->getBoundary();
if(bd_c.size() == 1 && inSameDomain(cell, bd_c.front()) ){
removeCell(cell);
removeCell(bd_c.front());
count++;
coreduced = true;
}
}
}
return count;
}
int CellComplex::reduction(int dim, bool deleteCells){
int CellComplex::reduction(int dim){
if(dim < 1 || dim > 3) return 0;
std::vector<Cell*> cbd_c;
std::vector<Cell*> bd_c;
int count = 0;
bool reduced = true;
......@@ -503,23 +688,49 @@ int CellComplex::reduction(int dim, bool deleteCells){
Cell* cell = *cit;
cbd_c = cbd(cell);
if(cbd_c.size() == 1 && inSameDomain(cell, cbd_c.at(0)) ){
removeCell(cell, deleteCells);
removeCell(cbd_c.at(0), deleteCells);
removeCell(cell);
removeCell(cbd_c.at(0));
count++;
reduced = true;
}
}
}
return count;
}
int CellComplex::reduction2(int dim){
if(dim < 1 || dim > 3) return 0;
std::list<Cell*> bd_c;
std::list<Cell*> cbd_c;
int count = 0;
bool reduced = true;
while (reduced){
reduced = false;
for(citer cit = firstCell(dim-1); cit != lastCell(dim-1); cit++){
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;
}
}
}
return count;
}
void CellComplex::reduceComplex(){
printf("Cell complex before reduction: %d volumes, %d faces, %d edges and %d vertices.\n",
getSize(3), getSize(2), getSize(1), getSize(0));
reduction(3);
reduction(2);
reduction(1);
reduction2(3);
reduction2(2);
reduction2(1);
printf("Cell complex after reduction: %d volumes, %d faces, %d edges and %d vertices.\n",
getSize(3), getSize(2), getSize(1), getSize(0));
}
......@@ -541,7 +752,7 @@ void CellComplex::coreduceComplex(int generatorDim, std::set<Cell*, Less_Cell>&
}
generatorCells[i].insert(cell);
removedCells.insert(cell);
removeCell(cell, false);
removeCell(cell);
coreduction(0, removedCells);
coreduction(1, removedCells);
coreduction(2, removedCells);
......@@ -575,11 +786,11 @@ void CellComplex::coreduceComplex(int generatorDim){
cit++;
}
generatorCells[i].insert(cell);
removeCell(cell, false);
removeCell(cell);
//coreduction(0);
//coreduction(1);
//coreduction(2);
//coreduction2(0);
//coreduction2(1);
//coreduction2(2);
coreductionMrozek(cell);
}
if(generatorDim == i) break;
......@@ -597,9 +808,12 @@ void CellComplex::coreduceComplex(int generatorDim){
return;
}
void CellComplex::repairComplex(){
void CellComplex::repairComplex(int i){
printf("Cell complex before repair: %d volumes, %d faces, %d edges and %d vertices.\n",
getSize(3), getSize(2), getSize(1), getSize(0));
for(int i = 3; i > 0; i--){
while(i > 0){
for(citer cit = firstCell(i); cit != lastCell(i); cit++){
Cell* cell = *cit;
......@@ -617,18 +831,13 @@ void CellComplex::repairComplex(){
}
}
i--;
}
/*
int tag = 1;
for(int i = 0; i < 4; i++){
for(citer cit = firstCell(i); cit != lastCell(i); cit++){
Cell* cell = *cit;
cell->setTag(tag);
tag++;
}
}
*/
printf("Cell complex after repair: %d volumes, %d faces, %d edges and %d vertices.\n",
getSize(3), getSize(2), getSize(1), getSize(0));
return;
}
......@@ -714,13 +923,13 @@ int CellComplex::writeComplexMSH(const std::string &name){
fprintf(fp, "$Nodes\n");
std::set<MVertex*, Less_MVertex> domainVertices;
getDomainVertices(domainVertices, true);
getDomainVertices(domainVertices, false);
//std::set<MVertex*, Less_MVertex> domainVertices;
//getDomainVertices(domainVertices, true);
//getDomainVertices(domainVertices, false);
fprintf(fp, "%d\n", domainVertices.size());
fprintf(fp, "%d\n", _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++){
MVertex* vertex = *vit;
fprintf(fp, "%d %.16g %.16g %.16g\n", vertex->getNum(), vertex->x(), vertex->y(), vertex->z() );
}
......
Geo/CellComplex.h
+ 137
− 49
View file @ c9f05afb
......@@ -33,9 +33,6 @@ class Cell
// maximum number of higher dimensional cells on the coboundary of this cell
int _cbdMaxSize;
int _bdSize;
int _cbdSize;
// whether this cell belongs to a subdomain
// used in relative homology computation
bool _inSubdomain;
......@@ -43,8 +40,13 @@ class Cell
// whether this cell belongs to the boundary of a cell complex
bool _onDomainBoundary;
// unique tag for each cell
int _tag;
// cells on the boundary and on the coboundary of thhis cell
std::list<Cell*> _boundary;
std::list<Cell*> _coboundary;
public:
Cell(){}
virtual ~Cell(){}
......@@ -70,10 +72,31 @@ class Cell
virtual unsigned int getBdMaxSize() const { return _bdMaxSize; };
virtual unsigned int getCbdMaxSize() const { return _cbdMaxSize; };
virtual int getBdSize() { return _bdSize; }
virtual void setBdSize(int size) { _bdSize = size; }
virtual int getCbdSize() { return _cbdSize; }
virtual void setCbdSize(int size) { _cbdSize = size; }
virtual int getBoundarySize() { return _boundary.size(); }
virtual int getCoboundarySize() { return _coboundary.size(); }
virtual void setBoundary(std::list<Cell*> boundary){ _boundary = boundary; }
virtual void setBoundary(std::vector<Cell*> boundary){
for(int i = 0; i < boundary.size(); i++) _boundary.push_back(boundary.at(i)); }
virtual void setCoboundary(std::list<Cell*> coboundary){ _coboundary = coboundary; }
virtual void setCoboundary(std::vector<Cell*> coboundary){
for(int i = 0; i < coboundary.size(); i++) _coboundary.push_back(coboundary.at(i)); }
virtual std::list<Cell*> getBoundary() { return _boundary; }
virtual std::list<Cell*> getCoboundary() { return _coboundary; }
virtual void addBoundaryCell(Cell* cell) { _boundary.push_back(cell); }
virtual void addCoboundaryCell(Cell* cell) { _coboundary.push_back(cell); }
virtual void removeBoundaryCell(Cell* cell) {
for(std::list<Cell*>::iterator it = _boundary.begin(); it != _boundary.end(); it++){
Cell* cell2 = *it;
if(cell2->getTag() == cell->getTag()) { _boundary.erase(it); break; }
}
}
virtual void removeCoboundaryCell(Cell* cell) {
for(std::list<Cell*>::iterator it = _coboundary.begin(); it != _coboundary.end(); it++){
Cell* cell2 = *it;
if(cell2->getTag() == cell->getTag()) { _coboundary.erase(it); break; }
}
}
virtual bool inSubdomain() const { return _inSubdomain; }
virtual void setInSubdomain(bool subdomain) { _inSubdomain = subdomain; }
......@@ -89,6 +112,9 @@ class Cell
virtual inline double y() const { return 0; }
virtual inline double z() const { return 0; }
virtual int getNumFacets() const { return 0; }
virtual void getFacetVertices(const int num, std::vector<int> &v) const {};
bool operator==(const Cell* c2){
if(this->getNumVertices() != c2->getNumVertices()){
return false;
......@@ -106,15 +132,13 @@ class Cell
return (this->getNumVertices() < c2->getNumVertices());
}
for(int i=0; i < this->getNumVertices();i++){
if(this->getVertex(i) < c2->getVertex(i)){
return true;
}
if(this->getVertex(i) < c2->getVertex(i)) return true;
else if (this->getVertex(i) > c2->getVertex(i)) return false;
}
return false;
}
};
// Simplicial cell type.
......@@ -123,8 +147,7 @@ class Simplex : public Cell
public:
Simplex() : Cell() {
_cbdSize = 1000; // big number
_bdSize = 1000;
}
virtual ~Simplex(){}
......@@ -203,17 +226,24 @@ class OneSimplex : public Simplex
_v[1] = dummy;
}
virtual ~OneSimplex(){}
~OneSimplex(){}
virtual int getDim() const { return 1; }
virtual int getNumVertices() const { return 2; }
virtual int getVertex(int vertex) const {return _v[vertex]; }
virtual bool hasVertex(int vertex) const {return (_v[0] == vertex || _v[1] == vertex); }
int getDim() const { return 1; }
int getNumVertices() const { return 2; }
int getNumFacets() const { return 2; }
int getVertex(int vertex) const {return _v[vertex]; }
bool hasVertex(int vertex) const {return (_v[0] == vertex || _v[1] == vertex); }
virtual std::vector<int> getVertexVector() const {
std::vector<int> getVertexVector() const {
return std::vector<int>(_v, _v + sizeof(_v)/sizeof(int) ); }
virtual void printCell() const { printf("Vertices: %d %d\n", _v[0], _v[1]); }
void getFacetVertices(const int num, std::vector<int> &v) const {
v.resize(1);
v[0] = _v[num];
}
void printCell() const { printf("Vertices: %d %d\n", _v[0], _v[1]); }
};
// Two simplex cell type.
......@@ -223,6 +253,16 @@ class TwoSimplex : public Simplex
// numbers of the vertices of this simplex
// same as the corresponding vertices in the finite element mesh
int _v[3];
int edges_tri(const int edge, const int vert) const{
static const int e[3][2] = {
{0, 1},
{1, 2},
{2, 0}
};
return e[edge][vert];
}
public:
TwoSimplex(std::vector<int> vertices, bool subdomain=false, bool boundary=false) : Simplex(){
......@@ -243,17 +283,24 @@ class TwoSimplex : public Simplex
_v[2] = dummy;
}
virtual ~TwoSimplex(){}
~TwoSimplex(){}
virtual int getDim() const { return 2; }
virtual int getNumVertices() const { return 3; }
virtual int getVertex(int vertex) const {return _v[vertex]; }
virtual bool hasVertex(int vertex) const {return
int getDim() const { return 2; }
int getNumVertices() const { return 3; }
int getNumFacets() const { return 3; }
int getVertex(int vertex) const {return _v[vertex]; }
bool hasVertex(int vertex) const {return
(_v[0] == vertex || _v[1] == vertex || _v[2] == vertex); }
virtual std::vector<int> getVertexVector() const {
std::vector<int> getVertexVector() const {
return std::vector<int>(_v, _v + sizeof(_v)/sizeof(int) ); }
virtual void printCell() const { printf("Vertices: %d %d %d\n", _v[0], _v[1], _v[2]); }
void getFacetVertices(const int num, std::vector<int> &v) const {
v.resize(2);
v[0] = _v[edges_tri(num, 0)];
v[1] = _v[edges_tri(num, 1)];
}
void printCell() const { printf("Vertices: %d %d %d\n", _v[0], _v[1], _v[2]); }
};
// Three simplex cell type.
......@@ -263,6 +310,17 @@ class ThreeSimplex : public Simplex
// numbers of the vertices of this simplex
// same as the corresponding vertices in the finite element mesh
int _v[4];
int faces_tetra(const int face, const int vert) const{
static const int f[4][3] = {
{0, 2, 1},
{0, 1, 3},
{0, 3, 2},
{3, 1, 2}
};
return f[face][vert];
}
public:
ThreeSimplex(std::vector<int> vertices, bool subdomain=false, bool boundary=false) : Simplex(){
......@@ -285,16 +343,24 @@ class ThreeSimplex : public Simplex
_v[3] = dummy;
}
virtual ~ThreeSimplex(){}
~ThreeSimplex(){}
virtual int getDim() const { return 3; }
virtual int getNumVertices() const { return 4; }
virtual int getVertex(int vertex) const {return _v[vertex]; }
virtual bool hasVertex(int vertex) const {return
int getDim() const { return 3; }
int getNumVertices() const { return 4; }
int getNumFacets() const { return 4; }
int getVertex(int vertex) const {return _v[vertex]; }
bool hasVertex(int vertex) const {return
(_v[0] == vertex || _v[1] == vertex || _v[2] == vertex || _v[3] == vertex); }
virtual std::vector<int> getVertexVector() const {
std::vector<int> getVertexVector() const {
return std::vector<int>(_v, _v + sizeof(_v)/sizeof(int) ); }
void getFacetVertices(const int num, std::vector<int> &v) const {
v.resize(3);
v[0] = _v[faces_tetra(num, 0)];
v[1] = _v[faces_tetra(num, 1)];
v[2] = _v[faces_tetra(num, 2)];
}
virtual void printCell() const { printf("Vertices: %d %d %d %d\n", _v[0], _v[1], _v[2], _v[3]); }
};
......@@ -323,15 +389,16 @@ class Less_Cell{
}
};
class Equal_Cell{
public:
bool operator()(const Cell* c1, const Cell* c2) const {
if(c1->getNumVertices() != c2->getNumVertices()) return false;
bool operator ()(const Cell* c1, const Cell* c2){
if(c1->getNumVertices() != c2->getNumVertices()){
return false;
}
for(int i=0; i < c1->getNumVertices();i++){
if(c1->getVertex(i) != c2->getVertex(i)) return false;
if(c1->getVertex(i) != c2->getVertex(i)){
return false;
}
}
return true;
}
......@@ -359,21 +426,25 @@ class CellComplex
std::vector<GEntity*> _boundary;
std::set<MVertex*, Less_MVertex> _domainVertices;
// sorted containers of unique cells in this cell complex
// one for each dimension
std::set<Cell*, Less_Cell> _cells[4];
std::set<Cell*, Less_Cell> _originalCells[4];
bool _simplicial;
public:
// iterator for the cells of same dimension
typedef std::set<Cell*, Less_Cell>::iterator citer;
protected:
// enqueue cells in queue if they are not there already
virtual void enqueueCells(std::vector<Cell*>& cells,
virtual void enqueueCells(std::list<Cell*>& cells,
std::queue<Cell*>& Q, std::set<Cell*, Less_Cell>& Qset);
virtual void enqueueCells(std::vector<Cell*>& cells, std::list<Cell*>& Q);
virtual void enqueueCells(std::list<Cell*>& cells, std::list<Cell*>& Q);
virtual void enqueueCellsIt(std::vector< std::set<Cell*, Less_Cell>::iterator >& cells,
std::queue<Cell*>& Q, std::set<Cell*, Less_Cell>& Qset);
......@@ -381,7 +452,8 @@ class CellComplex
virtual void removeCellQset(Cell*& cell, std::set<Cell*, Less_Cell>& Qset);
// insert cells into this cell complex
virtual void insertCells(bool subdomain, bool boundary);
virtual void insert_cells(bool subdomain, bool boundary);
virtual void insert_cells2(bool subdomain, bool boundary);
public:
CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> subdomain, std::set<Cell*, Less_Cell> cells ) {
......@@ -395,6 +467,8 @@ class CellComplex
CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> subdomain );
virtual ~CellComplex(){}
virtual bool simplicial() { return _simplicial; }
// get the number of certain dimensional cells
virtual int getSize(int dim){ return _cells[dim].size(); }
......@@ -413,7 +487,7 @@ class CellComplex
// kappa for two cells of this cell complex
// implementation will vary depending on cell type
virtual int kappa(Cell* sigma, Cell* tau) const { return sigma->kappa(tau); }
virtual inline int kappa(Cell* sigma, Cell* tau) const { return sigma->kappa(tau); }
// cells on the boundary of a cell
virtual std::vector<Cell*> bd(Cell* sigma);
......@@ -423,7 +497,7 @@ class CellComplex
virtual std::vector< std::set<Cell*, Less_Cell>::iterator > cbdIt(Cell* tau);
// remove a cell from this cell complex
virtual void removeCell(Cell* cell, bool deleteCell=false);
virtual void removeCell(Cell* cell);
virtual void removeCellIt(std::set<Cell*, Less_Cell>::iterator cell);
virtual void insertCell(Cell* cell);
......@@ -434,12 +508,14 @@ class CellComplex
// coreduction of this cell complex
// removes corection pairs of cells of dimension dim and dim+1
virtual int coreduction(int dim, bool deleteCells=false);
virtual int coreduction(int dim);
virtual int coreduction2(int dim);
// stores removed cells
virtual int coreduction(int dim, std::set<Cell*, Less_Cell>& removedCells);
// reduction of this cell complex
// removes reduction pairs of cell of dimension dim and dim-1
virtual int reduction(int dim, bool deleteCells=false);
virtual int reduction(int dim);
virtual int reduction2(int dim);
// useful functions for (co)reduction of cell complex
virtual void reduceComplex();
......@@ -456,10 +532,22 @@ class CellComplex
virtual int coreductionMrozek2(Cell* generator);
virtual int coreductionMrozek3(Cell* generator);
virtual void restoreComplex(){ for(int i = 0; i < 4; i++) _cells[i] = _originalCells[i]; }
// never use, O(n^2) complexity
virtual void restoreComplex(){
for(int i = 0; i < 4; i++) _cells[i] = _originalCells[i];
for(int i = 0; i < 4; i++){
for(citer cit = firstCell(i); cit != lastCell(i); cit++){
Cell* cell = *cit;
cell->setBoundary(bd(cell));
cell->setCoboundary(cbd(cell));
}
}
}
// add every volume, face and edge its missing boundary cells
virtual void repairComplex();
virtual void repairComplex(int i=3);
// change non-subdomain cells to be in subdomain, subdomain cells to not to be in subdomain
virtual void swapSubdomain();
......
Geo/ChainComplex.cpp
+ 2
− 2
View file @ c9f05afb
Geo/Makefile
+ 1
− 1
View file @ c9f05afb
utils/misc/celldriver.cpp
+ 7
− 6
View file @ c9f05afb
......@@ -93,15 +93,16 @@ int main(int argc, char **argv)
}
delete chains;
delete complex;
// restore the cell complex to its prereduced state
complex->restoreComplex();
// create new complex to compute the cuts
CellComplex* complex2 = new CellComplex(domain, subdomain);
// reduce the complex by coreduction, this removes all vertices, hence 0 as an argument
complex->coreduceComplex(0);
complex2->coreduceComplex(0);
// create a chain complex of a cell complex (construct boundary operator matrices)
ChainComplex* dualChains = new ChainComplex(complex);
ChainComplex* dualChains = new ChainComplex(complex2);
// transpose the boundary operator matrices,
// these are the boundary operator matrices for the dual chain complex
......@@ -120,7 +121,7 @@ int main(int argc, char **argv)
convert(3-(j-1), dimension);
std::string name = dimension + "D Dual cut " + generator;
Chain* chain = new Chain(complex->getCells(j-1), dualChains->getCoeffVector(j,i), complex, name);
Chain* chain = new Chain(complex2->getCells(j-1), dualChains->getCoeffVector(j,i), complex2, name);
chain->writeChainMSH("chains.msh");
delete chain;
}
......@@ -129,7 +130,7 @@ int main(int argc, char **argv)
delete dualChains;
delete complex;
delete complex2;
delete m;
GmshFinalize();
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment