diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 24616c6e7095da78c3ea33e2cf18a0b2e578a7cc..1d6337ad85f06ba81cd9a6c05fed2b745a7197cb 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -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,14 +170,81 @@ 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;
}
-
+
}
}
}
+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,19 +299,41 @@ std::vector<Cell*> CellComplex::bd(Cell* sigma){
int dim = sigma->getDim();
if(dim < 1) return boundary;
-
- 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){
- boundary.push_back(*cit);
- if(boundary.size() == sigma->getBdMaxSize()){
- 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;
+
+ 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;
+ }
}
}
}
@@ -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){
- _cells[cell->getDim()].erase(cell);
- //if(deleteCell) delete cell;
+void CellComplex::removeCell(Cell* cell){
+ _cells[cell->getDim()].erase(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,35 +481,37 @@ 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);
}
-
+
}
//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){
- for(int i = 3; i > 0; i--){
+ printf("Cell complex before repair: %d volumes, %d faces, %d edges and %d vertices.\n",
+ getSize(3), getSize(2), getSize(1), getSize(0));
+
+ 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() );
}
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index 587bfd68bde24bf0136a6a15ad526e16445b1685..ebeafaedec417de6a7c084a95c0cab501b4e62bf 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -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();
@@ -455,11 +531,23 @@ class CellComplex
virtual int coreductionMrozek(Cell* generator);
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();
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index 524e6235f87c050b1e15f96ea77a11f607791693..c0eae87ce57ebaf0f4487d95a5b1532a97c6b1cf 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -279,7 +279,7 @@ void ChainComplex::computeHomology(bool dual){
gmp_matrix_rows(getKerHMatrix(lowDim)),
gmp_matrix_cols(getKerHMatrix(lowDim))) );
}
-
+
// 5) General case:
// 1) Find the bases of boundaries B and cycles Z
// 2) find j: B -> Z and
@@ -321,7 +321,7 @@ void ChainComplex::computeHomology(bool dual){
setCodHMatrix(lowDim, NULL);
setJMatrix(lowDim, NULL);
setQMatrix(lowDim, NULL);
-
+
}
diff --git a/Geo/Makefile b/Geo/Makefile
index b834b47e646fbf89bffbe341fe5df24453129bd1..248d3aa02dad6f2ab6526cb5a2ae96503a53eddc 100644
--- a/Geo/Makefile
+++ b/Geo/Makefile
@@ -47,7 +47,7 @@ OBJ = ${SRC:.cpp=${OBJEXT}}
.SUFFIXES: ${OBJEXT} .cpp
-${LIB}: ${OBJ}
+${LIB}: ${OBJ}
${AR} ${ARFLAGS}${LIB} ${OBJ}
${RANLIB} ${LIB}
diff --git a/utils/misc/celldriver.cpp b/utils/misc/celldriver.cpp
index d54683dfd4800f9f1bf64da09b3fc7fcda52f1aa..14040cae5b07d901a12649106f1a098448791f7f 100755
--- a/utils/misc/celldriver.cpp
+++ b/utils/misc/celldriver.cpp
@@ -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();