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

Code cleanup. 

Added cell combining methods (about O(n^2)) to further reduce the size of
the cell complex before applying Smith normal form computation.

(Don't know why there seems to be some benchmarks/ files modified too,
Somebody else committed between my update and commit?)
parent 6f56a324
No related branches found
No related tags found
No related merge requests found
Expand all Show whitespace changes
Inline Side-by-side
Geo/CellComplex.cpp
+ 334
− 500
View file @ 6091a8d3
...@@ -70,9 +70,9 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su ...@@ -70,9 +70,9 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su
// insert cells into cell complex // insert cells into cell complex
// subdomain need to be inserted first! // subdomain need to be inserted first!
insert_cells2(true, true); insert_cells(true, true);
insert_cells2(false, true); insert_cells(false, true);
insert_cells2(false, false); insert_cells(false, false);
// find mesh vertices in the domain // find mesh vertices in the domain
for(unsigned int j=0; j < domain.size(); j++) { for(unsigned int j=0; j < domain.size(); j++) {
...@@ -93,94 +93,20 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su ...@@ -93,94 +93,20 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su
} }
// attach induvivial tags to cells // attach individual tags to cells
int tag = 1; int tag = 1;
for(int i = 0; i < 4; i++){ for(int i = 0; i < 4; i++){
for(citer cit = firstCell(i); cit != lastCell(i); cit++){ for(citer cit = firstCell(i); cit != lastCell(i); cit++){
Cell* cell = *cit; Cell* cell = *cit;
cell->setTag(tag); cell->setTag(tag);
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]; _originalCells[i] = _cells[i];
} }
_simplicial = true;
}
void CellComplex::insert_cells(bool subdomain, bool boundary){
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;
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
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){
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);
}
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++){
vertices.clear();
for(int l=0; l < domain.at(j)->getMeshElement(i)->getEdge(k).getNumVertices(); l++){
vertex = domain.at(j)->getMeshElement(i)->getEdge(k).getVertex(l)->getNum();
vertices.push_back(vertex);
}
Cell* cell = new OneSimplex(vertices, subdomain, boundary); // Add edges
insertInfo = _cells[1].insert(cell);
if(!insertInfo.second) delete cell;
}
}
} }
} void CellComplex::insert_cells(bool subdomain, bool boundary){
void CellComplex::insert_cells2(bool subdomain, bool boundary){
// works only for simplcial meshes at the moment! // works only for simplcial meshes at the moment!
...@@ -208,10 +134,10 @@ void CellComplex::insert_cells2(bool subdomain, bool boundary){ ...@@ -208,10 +134,10 @@ void CellComplex::insert_cells2(bool subdomain, bool boundary){
int dim = domain.at(j)->getMeshElement(i)->getDim(); int dim = domain.at(j)->getMeshElement(i)->getDim();
int type = domain.at(j)->getMeshElement(i)->getTypeForMSH(); int type = domain.at(j)->getMeshElement(i)->getTypeForMSH();
Cell* cell; Cell* cell;
if(dim == 3) cell = new ThreeSimplex(vertices, subdomain, boundary); if(dim == 3) cell = new ThreeSimplex(vertices, 0, subdomain, boundary);
else if(dim == 2) cell = new TwoSimplex(vertices, subdomain, boundary); else if(dim == 2) cell = new TwoSimplex(vertices, 0, subdomain, boundary);
else if(dim == 1) cell = new OneSimplex(vertices, subdomain, boundary); else if(dim == 1) cell = new OneSimplex(vertices, 0, subdomain, boundary);
else cell = new ZeroSimplex(vertices.at(0), subdomain, boundary); else cell = new ZeroSimplex(vertices.at(0), 0, subdomain, boundary);
insertInfo = _cells[dim].insert(cell); insertInfo = _cells[dim].insert(cell);
if(!insertInfo.second) delete cell; if(!insertInfo.second) delete cell;
} }
...@@ -225,19 +151,23 @@ void CellComplex::insert_cells2(bool subdomain, bool boundary){ ...@@ -225,19 +151,23 @@ void CellComplex::insert_cells2(bool subdomain, bool boundary){
for(int i = 0; i < cell->getNumFacets(); i++){ for(int i = 0; i < cell->getNumFacets(); i++){
cell->getFacetVertices(i, vertices); cell->getFacetVertices(i, vertices);
Cell* newCell; Cell* newCell;
if(dim == 3) newCell = new TwoSimplex(vertices, subdomain, boundary); if(dim == 3) newCell = new TwoSimplex(vertices, 0, subdomain, boundary);
else if(dim == 2) newCell = new OneSimplex(vertices, subdomain, boundary); else if(dim == 2) newCell = new OneSimplex(vertices, 0, subdomain, boundary);
else if(dim == 1) newCell = new ZeroSimplex(vertices.at(0), subdomain, boundary); else if(dim == 1) newCell = new ZeroSimplex(vertices.at(0), 0, subdomain, boundary);
insertInfo = _cells[dim-1].insert(newCell); insertInfo = _cells[dim-1].insert(newCell);
if(!insertInfo.second){ if(!insertInfo.second){
delete newCell; delete newCell;
Cell* oldCell = *(insertInfo.first); Cell* oldCell = *(insertInfo.first);
oldCell->addCoboundaryCell(cell); if(!subdomain && !boundary){
cell->addBoundaryCell(oldCell); int ori = cell->kappa(oldCell);
oldCell->addCoboundaryCell( ori, cell );
cell->addBoundaryCell( ori, oldCell);
} }
else{ }
cell->addBoundaryCell(newCell); else if(!subdomain && !boundary) {
newCell->addCoboundaryCell(cell); int ori = cell->kappa(newCell);
cell->addBoundaryCell( ori, newCell );
newCell->addCoboundaryCell( ori, cell);
} }
} }
} }
...@@ -249,6 +179,7 @@ void CellComplex::insertCell(Cell* cell){ ...@@ -249,6 +179,7 @@ void CellComplex::insertCell(Cell* cell){
_cells[cell->getDim()].insert(cell); _cells[cell->getDim()].insert(cell);
} }
int Simplex::kappa(Cell* tau) const{ int Simplex::kappa(Cell* tau) const{
for(int i=0; i < tau->getNumVertices(); i++){ for(int i=0; i < tau->getNumVertices(); i++){
if( !(this->hasVertex(tau->getVertex(i))) ) return 0; if( !(this->hasVertex(tau->getVertex(i))) ) return 0;
...@@ -258,158 +189,120 @@ int Simplex::kappa(Cell* tau) const{ ...@@ -258,158 +189,120 @@ int Simplex::kappa(Cell* tau) const{
int value=1; int value=1;
for(int i=0; i < tau->getNumVertices(); i++){ for(int i=0; i < tau->getNumVertices(); i++){
if(this->getVertex(i) != tau->getVertex(i)) return value; if(this->getSortedVertex(i) != tau->getSortedVertex(i)) return value;
value = value*-1; value = value*-1;
} }
return value; return value;
} }
std::set<Cell*, Less_Cell>::iterator CellComplex::findCell(int dim, std::vector<int>& vertices, bool original){
Cell* cell;
if(dim == 0) cell = new ZeroSimplex(vertices.at(0)); /*
else if(dim == 1) cell = new OneSimplex(vertices); int Simplex::kappa(Cell* tau) const{
else if(dim == 2) cell = new TwoSimplex(vertices); for(int i=0; i < tau->getNumVertices(); i++){
else cell = new ThreeSimplex(vertices); if( !(this->hasVertex(tau->getVertex(i))) ) return 0;
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){
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;
} }
if(this->getDim() - tau->getDim() != 1) return 0;
std::vector<Cell*> CellComplex::bd(Cell* sigma){ int value=1;
std::vector<Cell*> boundary;
int dim = sigma->getDim();
if(dim < 1) return boundary;
if(simplicial()){ // faster for(int i = 0; i < this->getNumFacets(); i++){
std::vector<int> vertices = sigma->getVertexVector(); std::vector<int> vTau = tau->getVertexVector();
std::vector<int> v;
this->getFacetVertices(i, v);
value = -1;
for(int j=0; j < vertices.size(); j++){ do {
std::vector<int> bVertices; value = value*-1;
if(v == vTau) return value;
// 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;
}
}
} }
while (std::next_permutation(vTau.begin(), vTau.end()) );
vTau = tau->getVertexVector();
value = -1;
do {
value = value*-1;
if(v == vTau) return value;
} }
else{ while (std::prev_permutation(vTau.begin(), vTau.end()) );
citer start = findCell(dim-1, sigma->getVertex(0), 0);
if(start == lastCell(dim-1)) return boundary;
citer end = findCell(dim-1, sigma->getVertex(sigma->getNumVertices()-1), sigma->getVertex(sigma->getNumVertices()-1));
if(end != lastCell(dim-1)) end++;
//for(citer cit = firstCell(dim-1); cit != lastCell(dim-1); cit++){
for(citer cit = start; cit != end; cit++){
//if(kappa(sigma, *cit) != 0){
if(sigma->kappa(*cit) != 0){
boundary.push_back(*cit);
if(boundary.size() == sigma->getBdMaxSize()){
return boundary;
}
}
}
} }
return boundary;
return 0;
} }
std::vector< std::set<Cell*, Less_Cell>::iterator > CellComplex::bdIt(Cell* sigma){
int dim = sigma->getDim(); int OneSimplex::kappa(Cell* tau) const{
std::vector< std::set<Cell*, Less_Cell>::iterator > boundary; for(int i=0; i < tau->getNumVertices(); i++){
if(dim < 1){ if( !(this->hasVertex(tau->getVertex(i))) ) return 0;
return boundary;
} }
citer start = findCell(dim-1, sigma->getVertex(0), 0); if(tau->getDim() != 0) return 0;
if(start == lastCell(dim-1)) return boundary;
citer end = findCell(dim-1, sigma->getVertex(sigma->getNumVertices()-1), sigma->getVertex(sigma->getNumVertices()-1)); if(tau->getVertex(0) == this->getVertex(0)) return -1;
if(end != lastCell(dim-1)) end++; else return 1;
for(citer cit = start; cit != end; cit++){
if(kappa(sigma, *cit) != 0){
boundary.push_back(cit);
if(boundary.size() == sigma->getBdMaxSize()){
return boundary;
}
}
} }
*/
return boundary; int Quadrangle::kappa(Cell* tau) const{
for(int i=0; i < tau->getNumVertices(); i++){
if( !(this->hasVertex(tau->getVertex(i))) ) return 0;
} }
if(tau->getDim() != 1) return 0;
int value=1;
std::vector<Cell*> CellComplex::cbd(Cell* tau){ std::vector<int> vTau = tau->getVertexVector();
std::vector<Cell*> coboundary;
int dim = tau->getDim();
if(dim > 2){
return coboundary;
}
for(citer cit = firstCell(dim+1); cit != lastCell(dim+1); cit++){ for(int i = 0; i < this->getNumFacets(); i++){
//if(kappa(*cit, tau) != 0){ std::vector<int> v;
Cell* cell = *cit; this->getFacetVertices(i, v);
if(cell->kappa(tau) != 0){ value = -1;
coboundary.push_back(*cit); do {
if(coboundary.size() == tau->getCbdMaxSize()){ value = value*-1;
return coboundary; if(v == vTau) return value;
}
} }
while (std::next_permutation(v.begin(), v.end()) );
} }
return coboundary; return value;
} }
std::vector< std::set<Cell*, Less_Cell>::iterator > CellComplex::cbdIt(Cell* tau){
std::vector< std::set<Cell*, Less_Cell>::iterator > coboundary; std::set<Cell*, Less_Cell>::iterator CellComplex::findCell(int dim, std::vector<int>& vertices, bool original){
int dim = tau->getDim(); Cell* cell;
if(dim > 2){
return coboundary;
}
for(citer cit = firstCell(dim+1); cit != lastCell(dim+1); cit++){ if(dim == 0) cell = new ZeroSimplex(vertices.at(0));
if(kappa(*cit, tau) != 0){ else if(dim == 1) cell = new OneSimplex(vertices);
coboundary.push_back(cit); else if(dim == 2) cell = new TwoSimplex(vertices);
if(coboundary.size() == tau->getCbdMaxSize()){ else cell = new ThreeSimplex(vertices);
return coboundary;
} citer cit;
} if(!original) cit = _cells[dim].find(cell);
else cit = _originalCells[dim].find(cell);
delete cell;
return cit;
} }
return coboundary; std::set<Cell*, Less_Cell>::iterator CellComplex::findCell(int dim, int vertex, int 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;
} }
void CellComplex::removeCell(Cell* cell){ void CellComplex::removeCell(Cell* cell){
_cells[cell->getDim()].erase(cell); _cells[cell->getDim()].erase(cell);
std::list<Cell*> coboundary = cell->getCoboundary(); std::list<Cell*> coboundary = cell->getCoboundary();
...@@ -419,34 +312,18 @@ void CellComplex::removeCell(Cell* cell){ ...@@ -419,34 +312,18 @@ void CellComplex::removeCell(Cell* cell){
Cell* cbdCell = *it; Cell* cbdCell = *it;
cbdCell->removeBoundaryCell(cell); cbdCell->removeBoundaryCell(cell);
} }
for(std::list<Cell*>::iterator it = boundary.begin(); it != boundary.end(); it++){ for(std::list<Cell*>::iterator it = boundary.begin(); it != boundary.end(); it++){
Cell* bdCell = *it; Cell* bdCell = *it;
bdCell->removeCoboundaryCell(cell); bdCell->removeCoboundaryCell(cell);
} }
}
void CellComplex::removeCellIt(std::set<Cell*, Less_Cell>::iterator cell){
Cell* c = *cell;
int dim = c->getDim();
_cells[dim].erase(cell);
} }
void CellComplex::removeCellQset(Cell*& cell, std::set<Cell*, Less_Cell>& Qset){ void CellComplex::removeCellQset(Cell*& cell, std::set<Cell*, Less_Cell>& Qset){
Qset.erase(cell); Qset.erase(cell);
} }
void CellComplex::enqueueCellsIt(std::vector< std::set<Cell*, Less_Cell>::iterator >& cells,
std::queue<Cell*>& Q, std::set<Cell*, Less_Cell>& Qset){
for(unsigned int i=0; i < cells.size(); i++){
Cell* cell = *cells.at(i);
citer cit = Qset.find(cell);
if(cit == Qset.end()){
Qset.insert(cell);
Q.push(cell);
}
}
}
void CellComplex::enqueueCells(std::list<Cell*>& cells, std::queue<Cell*>& Q, std::set<Cell*, Less_Cell>& Qset){ 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++){ for(std::list<Cell*>::iterator cit = cells.begin(); cit != cells.end(); cit++){
Cell* cell = *cit; Cell* cell = *cit;
...@@ -458,18 +335,6 @@ void CellComplex::enqueueCells(std::list<Cell*>& cells, std::queue<Cell*>& Q, st ...@@ -458,18 +335,6 @@ void CellComplex::enqueueCells(std::list<Cell*>& cells, std::queue<Cell*>& Q, st
} }
} }
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(cell);
}
}
}
int CellComplex::coreductionMrozek(Cell* generator){ int CellComplex::coreductionMrozek(Cell* generator){
int coreductions = 0; int coreductions = 0;
...@@ -482,6 +347,7 @@ int CellComplex::coreductionMrozek(Cell* generator){ ...@@ -482,6 +347,7 @@ int CellComplex::coreductionMrozek(Cell* generator){
//removeCell(generator); //removeCell(generator);
std::list<Cell*> bd_s; std::list<Cell*> bd_s;
//std::list< std::pair<int, Cell*> >bd_s;
std::list<Cell*> cbd_c; std::list<Cell*> cbd_c;
Cell* s; Cell* s;
int round = 0; int round = 0;
...@@ -512,145 +378,11 @@ int CellComplex::coreductionMrozek(Cell* generator){ ...@@ -512,145 +378,11 @@ int CellComplex::coreductionMrozek(Cell* generator){
return coreductions; return coreductions;
} }
int CellComplex::coreductionMrozek2(Cell* generator){
int coreductions = 0;
std::list<Cell*> Q;
Q.push_back(generator);
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_front();
bd_s = s->getBoundary();
if( bd_s.size() == 1 && inSameDomain(s, bd_s.front()) ){
removeCell(s);
cbd_c = bd_s.front()->getCoboundary();
enqueueCells(cbd_c, Q);
removeCell(bd_s.front());
coreductions++;
}
else if(bd_s.empty()){
cbd_c = s->getCoboundary();
enqueueCells(cbd_c, Q);
}
//Q.unique(Equal_Cell());
}
//printf("Coreduction: %d loops with %d coreductions\n", round, coreductions);
return coreductions;
}
int CellComplex::coreductionMrozek3(Cell* generator){
int coreductions = 0;
std::queue<Cell*> Q;
std::set<Cell*, Less_Cell> Qset;
Q.push(generator);
Qset.insert(generator);
std::vector< std::set<Cell*, Less_Cell>::iterator > bd_s;
std::vector< std::set<Cell*, Less_Cell>::iterator > cbd_c;
Cell* s;
int round = 0;
while( !Q.empty() ){
round++;
//printf("%d ", round);
s = Q.front();
Q.pop();
removeCellQset(s, Qset);
bd_s = bdIt(s);
if( bd_s.size() == 1 && inSameDomain(s, *bd_s.at(0)) ){
removeCell(s);
cbd_c = cbdIt(*bd_s.at(0));
enqueueCellsIt(cbd_c, Q, Qset);
removeCellIt(bd_s.at(0));
coreductions++;
}
else if(bd_s.empty()){
cbd_c = cbdIt(s);
enqueueCellsIt(cbd_c, Q, Qset);
}
}
//printf("Coreduction: %d loops with %d coreductions\n", round, coreductions);
return coreductions;
}
int CellComplex::coreduction(int dim){ int CellComplex::coreduction(int dim){
if(dim < 0 || dim > 2) return 0; if(dim < 0 || dim > 2) return 0;
std::vector<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 = bd(cell);
if(bd_c.size() == 1 && inSameDomain(cell, bd_c.at(0)) ){
removeCell(cell);
removeCell(bd_c.at(0));
count++;
coreduced = true;
}
}
}
return count;
}
int CellComplex::coreduction(int dim, std::set<Cell*, Less_Cell>& removedCells){
if(dim < 0 || dim > 2) return 0;
std::vector<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 = bd(cell);
if(bd_c.size() == 1 && inSameDomain(cell, bd_c.at(0)) ){
removedCells.insert(cell);
removedCells.insert(bd_c.at(0));
removeCell(cell);
removeCell(bd_c.at(0));
count++;
coreduced = true;
}
}
}
return count;
}
int CellComplex::coreduction2(int dim){
if(dim < 0 || dim > 2) return 0;
std::list<Cell*> bd_c; std::list<Cell*> bd_c;
int count = 0; int count = 0;
bool coreduced = true; bool coreduced = true;
...@@ -669,40 +401,14 @@ int CellComplex::coreduction2(int dim){ ...@@ -669,40 +401,14 @@ int CellComplex::coreduction2(int dim){
} }
} }
return count;
}
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;
while (reduced){
reduced = false;
for(citer cit = firstCell(dim-1); cit != lastCell(dim-1); cit++){
Cell* cell = *cit;
cbd_c = cbd(cell);
if(cbd_c.size() == 1 && inSameDomain(cell, cbd_c.at(0)) ){
removeCell(cell);
removeCell(cbd_c.at(0));
count++;
reduced = true;
}
}
}
return count; return count;
}
}
int CellComplex::reduction2(int dim){ int CellComplex::reduction(int dim){
if(dim < 1 || dim > 3) return 0; if(dim < 1 || dim > 3) return 0;
std::list<Cell*> bd_c;
std::list<Cell*> cbd_c; std::list<Cell*> cbd_c;
int count = 0; int count = 0;
...@@ -728,22 +434,22 @@ int CellComplex::reduction2(int dim){ ...@@ -728,22 +434,22 @@ int CellComplex::reduction2(int dim){
void CellComplex::reduceComplex(){ void CellComplex::reduceComplex(){
printf("Cell complex before reduction: %d volumes, %d faces, %d edges and %d vertices.\n", printf("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));
reduction2(3); reduction(3);
reduction2(2); reduction(2);
reduction2(1); reduction(1);
printf("Cell complex after reduction: %d volumes, %d faces, %d edges and %d vertices.\n", printf("Cell complex after 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));
} }
void CellComplex::coreduceComplex(int generatorDim, std::set<Cell*, Less_Cell>& removedCells){ void CellComplex::coreduceComplex(int generatorDim){
std::set<Cell*, Less_Cell> generatorCells[4]; std::set<Cell*, Less_Cell> generatorCells[4];
printf("Cell complex before coreduction: %d volumes, %d faces, %d edges and %d vertices.\n", printf("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));
for(int i = 0; i < 4; i++){
int i = generatorDim;
while (getSize(i) != 0){ while (getSize(i) != 0){
//if(generatorDim == i || i == generatorDim+1) break;
citer cit = firstCell(i); citer cit = firstCell(i);
Cell* cell = *cit; Cell* cell = *cit;
while(!cell->inSubdomain() && cit != lastCell(i)){ while(!cell->inSubdomain() && cit != lastCell(i)){
...@@ -751,63 +457,212 @@ void CellComplex::coreduceComplex(int generatorDim, std::set<Cell*, Less_Cell>& ...@@ -751,63 +457,212 @@ void CellComplex::coreduceComplex(int generatorDim, std::set<Cell*, Less_Cell>&
cit++; cit++;
} }
generatorCells[i].insert(cell); generatorCells[i].insert(cell);
removedCells.insert(cell);
removeCell(cell); removeCell(cell);
coreduction(0, removedCells);
coreduction(1, removedCells);
coreduction(2, removedCells);
}
//coreduction(0);
//coreduction(1);
//coreduction(2);
coreductionMrozek(cell);
}
if(generatorDim == i) break;
}
/*
for(int i = 0; i < 4; i++){
for(citer cit = generatorCells[i].begin(); cit != generatorCells[i].end(); cit++){ for(citer cit = generatorCells[i].begin(); cit != generatorCells[i].end(); cit++){
Cell* cell = *cit; Cell* cell = *cit;
if(!cell->inSubdomain()) _cells[i].insert(cell); if(!cell->inSubdomain()) _cells[i].insert(cell);
if(!cell->inSubdomain()) removedCells.insert(cell);
} }
}
*/
printf("Cell complex after coreduction: %d volumes, %d faces, %d edges and %d vertices.\n", printf("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; return;
}
void CellComplex::replaceCells(Cell* c1, Cell* c2, Cell* newCell, bool orMatch, bool co){
int dim = c1->getDim();
std::list< std::pair<int, Cell*> > coboundary1 = c1->getOrientedCoboundary();
std::list< std::pair<int, Cell*> > coboundary2 = c2->getOrientedCoboundary();
std::list< std::pair<int, Cell*> > boundary1 = c1->getOrientedBoundary();
std::list< std::pair<int, Cell*> > boundary2 = c2->getOrientedBoundary();
for(std::list< std::pair<int, Cell*> >::iterator it = coboundary1.begin(); it != coboundary1.end(); it++){
Cell* cbdCell = (*it).second;
int ori = (*it).first;
cbdCell->removeBoundaryCell(c1);
cbdCell->addBoundaryCell(ori, newCell);
}
for(std::list< std::pair<int, Cell*> >::iterator it = coboundary2.begin(); it != coboundary2.end(); it++){
Cell* cbdCell = (*it).second;
int ori = (*it).first;
if(!orMatch) ori = -1*ori;
cbdCell->removeBoundaryCell(c2);
if(!co){
bool old = false;
for(std::list< std::pair<int, Cell* > >::iterator it2 = coboundary1.begin(); it2 != coboundary1.end(); it2++){
Cell* cell2 = (*it2).second;
if(*cell2 == *cbdCell) old = true;
}
if(!old) cbdCell->addBoundaryCell(ori, newCell);
}
else cbdCell->addBoundaryCell(ori, newCell);
} }
void CellComplex::coreduceComplex(int generatorDim){
std::set<Cell*, Less_Cell> generatorCells[4]; 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);
}
for(std::list< std::pair<int, Cell* > >::iterator it = boundary2.begin(); it != boundary2.end(); it++){
Cell* bdCell = (*it).second;
int ori = (*it).first;
if(!orMatch) ori = -1*ori;
bdCell->removeCoboundaryCell(c2);
if(co){
bool old = false;
for(std::list< std::pair<int, Cell* > >::iterator it2 = boundary1.begin(); it2 != boundary1.end(); it2++){
Cell* cell2 = (*it2).second;
if(*cell2 == *bdCell) old = true;
}
if(!old) bdCell->addCoboundaryCell(ori, newCell);
}
else bdCell->addCoboundaryCell(ori, newCell);
printf("Cell complex before coreduction: %d volumes, %d faces, %d edges and %d vertices.\n", }
_cells[dim].erase(c1);
_cells[dim].erase(c2);
_cells[dim].insert(newCell);
return;
}
int CellComplex::cocombine(int dim){
printf("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));
for(int i = 0; i < 4; i++){
while (getSize(i) != 0){ if(dim < 1 || dim > 3) return 0;
//if(generatorDim == i || i == generatorDim+1) break;
citer cit = firstCell(i); std::queue<Cell*> Q;
std::set<Cell*, Less_Cell> Qset;
std::list<Cell*> cbd_c;
std::list< std::pair< int, Cell*> > bd_c;
int count = 0;
for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
Cell* cell = *cit; Cell* cell = *cit;
while(!cell->inSubdomain() && cit != lastCell(i)){ cbd_c = cell->getCoboundary();
cell = *cit; enqueueCells(cbd_c, Q, Qset);
cit++; while(Q.size() != 0){
Cell* s = Q.front();
Q.pop();
bd_c = s->getOrientedBoundary();
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) ){
Cell* c1 = bd_c.front().second;
Cell* c2 = bd_c.back().second;
cbd_c = c1->getCoboundary();
enqueueCells(cbd_c, Q, Qset);
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++;
} }
generatorCells[i].insert(cell); removeCellQset(s, Qset);
removeCell(cell);
//coreduction2(0);
//coreduction2(1);
//coreduction2(2);
coreductionMrozek(cell);
} }
if(generatorDim == i) break; }
printf("Cell complex after cocombining: %d volumes, %d faces, %d edges and %d vertices.\n",
getSize(3), getSize(2), getSize(1), getSize(0));
return count;
} }
for(int i = 0; i < 4; i++){
for(citer cit = generatorCells[i].begin(); cit != generatorCells[i].end(); cit++){ int CellComplex::combine(int dim){
printf("Cell complex before combining: %d volumes, %d faces, %d edges and %d vertices.\n",
getSize(3), getSize(2), getSize(1), getSize(0));
if(dim < 1 || dim > 3) return 0;
std::queue<Cell*> Q;
std::set<Cell*, Less_Cell> Qset;
std::list< std::pair<int, Cell*> > cbd_c;
std::list<Cell*> bd_c;
int count = 0;
for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
Cell* cell = *cit; Cell* cell = *cit;
//if(!cell->inSubdomain()) _cells[i].insert(cell); bd_c = cell->getBoundary();
enqueueCells(bd_c, Q, Qset);
while(Q.size() != 0){
Cell* s = Q.front();
Q.pop();
cbd_c = s->getOrientedCoboundary();
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) ){
Cell* c1 = cbd_c.front().second;
Cell* c2 = cbd_c.back().second;
bd_c = c1->getBoundary();
enqueueCells(bd_c, Q, Qset);
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);
count++;
} }
removeCellQset(s, Qset);
} }
printf("Cell complex after coreduction: %d volumes, %d faces, %d edges and %d vertices.\n", }
printf("Cell complex after 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));
return;
return count;
} }
void CellComplex::repairComplex(int i){ void CellComplex::repairComplex(int i){
printf("Cell complex before repair: %d volumes, %d faces, %d edges and %d vertices.\n", printf("Cell complex before repair: %d volumes, %d faces, %d edges and %d vertices.\n",
...@@ -851,61 +706,24 @@ void CellComplex::swapSubdomain(){ ...@@ -851,61 +706,24 @@ void CellComplex::swapSubdomain(){
} }
} }
return; // make subdomain closed
} for(int i = 0; i < 4; i++){
for(citer cit = firstCell(i); cit != lastCell(i); cit++){
Cell* cell = *cit;
void CellComplex::getDomainVertices(std::set<MVertex*, Less_MVertex>& domainVertices, bool subdomain){ if(cell->inSubdomain()){
std::list<Cell*> boundary = cell->getBoundary();
std::vector<GEntity*> domain; for(std::list<Cell*>::iterator it = boundary.begin(); it != boundary.end(); it++){
if(subdomain) domain = _subdomain; Cell* bdCell = *it;
else domain = _domain; bdCell->setInSubdomain(true);
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 < _domain.size(); j++) {
for(unsigned int i=0; i < _domain.at(j)->getNumMeshVertices(); i++){
MVertex* vertex = _domain.at(j)->getMeshVertex(i);
domainVertices.insert(vertex);
}
std::list<GFace*> faces = _domain.at(j)->faces();
for(std::list<GFace*>::iterator fit = faces.begin(); fit != faces.end(); fit++){
GFace* face = *fit;
for(unsigned int i=0; i < face->getNumMeshVertices(); i++){
MVertex* vertex = face->getMeshVertex(i);
domainVertices.insert(vertex);
}
}
std::list<GEdge*> edges = _domain.at(j)->edges();
for(std::list<GEdge*>::iterator eit = edges.begin(); eit != edges.end(); eit++){
GEdge* edge = *eit;
for(unsigned int i=0; i < edge->getNumMeshVertices(); i++){
MVertex* vertex = edge->getMeshVertex(i);
domainVertices.insert(vertex);
} }
} }
std::list<GVertex*> vertices = _domain.at(j)->vertices();
for(std::list<GVertex*>::iterator vit = vertices.begin(); vit != vertices.end(); vit++){
GVertex* vertex = *vit;
for(unsigned int i=0; i < vertex->getNumMeshVertices(); i++){
MVertex* mvertex = vertex->getMeshVertex(i);
domainVertices.insert(mvertex);
} }
} }
}
*/
return; return;
} }
int CellComplex::writeComplexMSH(const std::string &name){ int CellComplex::writeComplexMSH(const std::string &name){
...@@ -923,10 +741,6 @@ int CellComplex::writeComplexMSH(const std::string &name){ ...@@ -923,10 +741,6 @@ int CellComplex::writeComplexMSH(const std::string &name){
fprintf(fp, "$Nodes\n"); fprintf(fp, "$Nodes\n");
//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++){
...@@ -938,7 +752,15 @@ int CellComplex::writeComplexMSH(const std::string &name){ ...@@ -938,7 +752,15 @@ int CellComplex::writeComplexMSH(const std::string &name){
fprintf(fp, "$EndNodes\n"); fprintf(fp, "$EndNodes\n");
fprintf(fp, "$Elements\n"); fprintf(fp, "$Elements\n");
fprintf(fp, "%d\n", _cells[0].size() + _cells[1].size() + _cells[2].size() + _cells[3].size()); int count = 0;
for(int i = 0; i < 4; i++){
for(citer cit = firstCell(i); cit != lastCell(i); cit++){
Cell* cell = *cit;
count = count + cell->getNumCells();
}
}
fprintf(fp, "%d\n", count);
int physical = 0; int physical = 0;
...@@ -950,13 +772,17 @@ int CellComplex::writeComplexMSH(const std::string &name){ ...@@ -950,13 +772,17 @@ int CellComplex::writeComplexMSH(const std::string &name){
fprintf(fp, "%d %d %d %d %d %d %d\n", vertex->getTag(), 15, 3, 0, 0, physical, vertex->getVertex(0)); fprintf(fp, "%d %d %d %d %d %d %d\n", vertex->getTag(), 15, 3, 0, 0, physical, vertex->getVertex(0));
} }
std::list< std::pair<int, Cell*> > cells;
for(citer cit = firstCell(1); cit != lastCell(1); cit++) { for(citer cit = firstCell(1); cit != lastCell(1); cit++) {
Cell* edge = *cit; Cell* edge = *cit;
if(edge->inSubdomain()) physical = 3; if(edge->inSubdomain()) physical = 3;
else if(edge->onDomainBoundary()) physical = 2; else if(edge->onDomainBoundary()) physical = 2;
else physical = 1; else physical = 1;
fprintf(fp, "%d %d %d %d %d %d %d %d\n", edge->getTag(), 1, 3, 0, 0, physical, edge->getVertex(0), edge->getVertex(1)); cells = edge->getCells();
for(std::list< std::pair<int, Cell*> >::iterator it = cells.begin(); it != cells.end(); it++){
Cell* cell = (*it).second;
fprintf(fp, "%d %d %d %d %d %d %d %d\n", cell->getTag(), 1, 3, 0, 0, physical, cell->getVertex(0), cell->getVertex(1));
}
} }
for(citer cit = firstCell(2); cit != lastCell(2); cit++) { for(citer cit = firstCell(2); cit != lastCell(2); cit++) {
...@@ -964,14 +790,22 @@ int CellComplex::writeComplexMSH(const std::string &name){ ...@@ -964,14 +790,22 @@ int CellComplex::writeComplexMSH(const std::string &name){
if(face->inSubdomain()) physical = 3; if(face->inSubdomain()) physical = 3;
else if(face->onDomainBoundary()) physical = 2; else if(face->onDomainBoundary()) physical = 2;
else physical = 1; else physical = 1;
fprintf(fp, "%d %d %d %d %d %d %d %d %d\n", face->getTag(), 2, 3, 0, 0, physical, face->getVertex(0), face->getVertex(1), face->getVertex(2)); cells = face->getCells();
for(std::list< std::pair<int, Cell*> >::iterator it = cells.begin(); it != cells.end(); it++){
Cell* cell = (*it).second;
fprintf(fp, "%d %d %d %d %d %d %d %d %d\n", cell->getTag(), 2, 3, 0, 0, physical, cell->getVertex(0), cell->getVertex(1), cell->getVertex(2));
}
} }
for(citer cit = firstCell(3); cit != lastCell(3); cit++) { for(citer cit = firstCell(3); cit != lastCell(3); cit++) {
Cell* volume = *cit; Cell* volume = *cit;
if(volume->inSubdomain()) physical = 3; if(volume->inSubdomain()) physical = 3;
else if(volume->onDomainBoundary()) physical = 2; else if(volume->onDomainBoundary()) physical = 2;
else physical = 1; else physical = 1;
fprintf(fp, "%d %d %d %d %d %d %d %d %d %d\n", volume->getTag(), 4, 3, 0, 0, physical, volume->getVertex(0), volume->getVertex(1), volume->getVertex(2), volume->getVertex(3)); cells = volume->getCells();
for(std::list< std::pair<int, Cell*> >::iterator it = cells.begin(); it != cells.end(); it++){
Cell* cell = (*it).second;
fprintf(fp, "%d %d %d %d %d %d %d %d %d %d\n", cell->getTag(), 4, 3, 0, 0, physical, cell->getVertex(0), cell->getVertex(1), cell->getVertex(2), cell->getVertex(3));
}
} }
fprintf(fp, "$EndElements\n"); fprintf(fp, "$EndElements\n");
...@@ -982,11 +816,11 @@ int CellComplex::writeComplexMSH(const std::string &name){ ...@@ -982,11 +816,11 @@ int CellComplex::writeComplexMSH(const std::string &name){
} }
void CellComplex::printComplex(int dim, bool subcomplex){ void CellComplex::printComplex(int 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;
for(int i = 0; i < cell->getNumVertices(); i++){ for(int i = 0; i < cell->getNumVertices(); i++){
if(!subcomplex && !cell->inSubdomain()) printf("%d ", cell->getVertex(i)); printf("%d ", cell->getVertex(i));
} }
printf("\n"); printf("\n");
} }
......
Geo/CellComplex.h
+ 311
− 112
View file @ 6091a8d3
This diff is collapsed.
Geo/ChainComplex.cpp
+ 63
− 9
View file @ 6091a8d3
...@@ -24,15 +24,26 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){ ...@@ -24,15 +24,26 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){
unsigned int rows = cellComplex->getSize(dim-1); unsigned int rows = cellComplex->getSize(dim-1);
int index = 1;
// ignore subdomain cells // ignore subdomain cells
for(std::set<Cell*, Less_Cell>::iterator cit = cellComplex->firstCell(dim); cit != cellComplex->lastCell(dim); cit++){ for(std::set<Cell*, Less_Cell>::iterator cit = cellComplex->firstCell(dim); cit != cellComplex->lastCell(dim); cit++){
Cell* cell = *cit; Cell* cell = *cit;
if(cell->inSubdomain()) cols--; cell->setIndex(index);
index++;
if(cell->inSubdomain()){
index--;
cols--;
}
} }
index = 1;
for(std::set<Cell*, Less_Cell>::iterator cit = cellComplex->firstCell(dim-1); cit != cellComplex->lastCell(dim-1); cit++){ for(std::set<Cell*, Less_Cell>::iterator cit = cellComplex->firstCell(dim-1); cit != cellComplex->lastCell(dim-1); cit++){
Cell* cell = *cit; Cell* cell = *cit;
if(cell->inSubdomain()) rows--; cell->setIndex(index);
index++;
if(cell->inSubdomain()){
index--;
rows--;
}
} }
...@@ -45,7 +56,34 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){ ...@@ -45,7 +56,34 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){
//_HMatrix[dim] = NULL; //_HMatrix[dim] = NULL;
} }
else{
mpz_t elem;
mpz_init(elem);
_HMatrix[dim] = create_gmp_matrix_zero(rows, cols);
for( std::set<Cell*, Less_Cell>::iterator cit = cellComplex->firstCell(dim); cit != cellComplex->lastCell(dim); cit++){
Cell* cell = *cit;
if(!cell->inSubdomain()){
std::list< std::pair<int,Cell*> >bdCell = cell->getOrientedBoundary();
for(std::list< std::pair<int, Cell*> >::iterator it = bdCell.begin(); it != bdCell.end(); it++){
Cell* bdCell = (*it).second;
if(!bdCell->inSubdomain()){
int old_elem = 0;
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);
gmp_matrix_set_elem(elem, bdCell->getIndex(), cell->getIndex(), _HMatrix[dim]);
}
}
}
}
mpz_clear(elem);
}
/*
else{ else{
long int elems[rows*cols]; long int elems[rows*cols];
...@@ -59,6 +97,16 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){ ...@@ -59,6 +97,16 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){
Cell* lowcell = *low; Cell* lowcell = *low;
Cell* highcell = *high; Cell* highcell = *high;
if(!(highcell->inSubdomain() || lowcell->inSubdomain())){ if(!(highcell->inSubdomain() || lowcell->inSubdomain())){
std::list< std::pair<int, Cell*> >bdHigh = highcell->getBoundary();
for(std::list< std::pair<int, Cell*> >::iterator it = bdHigh.begin(); it != bdHigh.end(); it++){
Cell* bdCell = (*it).second;
if(bdCell->getTag() == lowcell->getTag()) elems[i] = (*it).first;
else elems[i] = 0;
}
elems[i] = cellComplex->kappa(*high, *low); elems[i] = cellComplex->kappa(*high, *low);
i++; i++;
} }
...@@ -69,6 +117,8 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){ ...@@ -69,6 +117,8 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){
} }
_HMatrix[dim] = create_gmp_matrix_int(rows, cols, elems); _HMatrix[dim] = create_gmp_matrix_int(rows, cols, elems);
} }
*/
_kerH[dim] = NULL; _kerH[dim] = NULL;
_codH[dim] = NULL; _codH[dim] = NULL;
...@@ -414,13 +464,17 @@ int Chain::writeChainMSH(const std::string &name){ ...@@ -414,13 +464,17 @@ int Chain::writeChainMSH(const std::string &name){
fprintf(fp, "4 \n"); fprintf(fp, "4 \n");
fprintf(fp, "0 \n"); fprintf(fp, "0 \n");
fprintf(fp, "1 \n"); fprintf(fp, "1 \n");
fprintf(fp, "%d \n", getSize()); fprintf(fp, "%d \n", getNumCells());
fprintf(fp, "0 \n"); fprintf(fp, "0 \n");
std::list< std::pair<int, Cell*> > cells;
for(int i = 0; i < getSize(); i++){ for(int i = 0; i < getSize(); i++){
Cell* cell = getCell(i);
fprintf(fp, "%d %d \n", getCell(i)->getTag(), getCoeff(i)); cells = cell->getCells();
for(std::list< std::pair<int, Cell*> >::iterator it = cells.begin(); it != cells.end(); it++){
Cell* cell2 = (*it).second;
fprintf(fp, "%d %d \n", cell2->getTag(), getCoeff(i)*(*it).first );
}
} }
fprintf(fp, "$EndElementData\n"); fprintf(fp, "$EndElementData\n");
......
Geo/ChainComplex.h
+ 10
− 0
View file @ 6091a8d3
...@@ -142,6 +142,16 @@ class Chain{ ...@@ -142,6 +142,16 @@ class Chain{
// number of cells in this chain // number of cells in this chain
virtual int getSize() { return _cells.size(); } virtual int getSize() { return _cells.size(); }
virtual int getNumCells() {
int count = 0;
for(std::vector< std::pair <Cell*, int> >::iterator it = _cells.begin(); it != _cells.end(); it++){
Cell* cell = (*it).first;
count = count + cell->getNumCells();
}
return count;
}
// get/set chain name // get/set chain name
virtual std::string getName() { return _name; } virtual std::string getName() { return _name; }
virtual void setName(std::string name) { _name=name; } virtual void setName(std::string name) { _name=name; }
......
utils/misc/celldriver.cpp
+ 9
− 0
View file @ 6091a8d3
...@@ -68,6 +68,11 @@ int main(int argc, char **argv) ...@@ -68,6 +68,11 @@ int main(int argc, char **argv)
// reduce the complex in order to ease homology computation // reduce the complex in order to ease homology computation
complex->reduceComplex(); complex->reduceComplex();
// perform series of cell combinatios in order to further ease homology computation
complex->combine(3);
complex->combine(2);
complex->combine(1);
// create a chain complex of a cell complex (construct boundary operator matrices) // create a chain complex of a cell complex (construct boundary operator matrices)
ChainComplex* chains = new ChainComplex(complex); ChainComplex* chains = new ChainComplex(complex);
// compute the homology using the boundary operator matrices // compute the homology using the boundary operator matrices
...@@ -101,6 +106,10 @@ int main(int argc, char **argv) ...@@ -101,6 +106,10 @@ int main(int argc, char **argv)
// reduce the complex by coreduction, this removes all vertices, hence 0 as an argument // reduce the complex by coreduction, this removes all vertices, hence 0 as an argument
complex2->coreduceComplex(0); complex2->coreduceComplex(0);
// perform series of "dual complex" cell combinations in order to ease homology computation
complex2->cocombine(1);
complex2->cocombine(2);
// create a chain complex of a cell complex (construct boundary operator matrices) // create a chain complex of a cell complex (construct boundary operator matrices)
ChainComplex* dualChains = new ChainComplex(complex2); ChainComplex* dualChains = new ChainComplex(complex2);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment