From a804c7a42fa51174e773bd6ce11f009de0f85a7f Mon Sep 17 00:00:00 2001
From: Matti Pellika <matti.pellikka@tut.fi>
Date: Tue, 5 May 2009 11:35:01 +0000
Subject: [PATCH] Added ability to compute the "scalar potential cuts" for
homology generators using coreduction. Modified utils/misc/celldriver.cpp to
demonstrate this.
Added some more or less useful utility and experimental functions.
Numerious bug fixes.
---
Geo/CellComplex.cpp | 285 ++++++++++++++++++++++++++++++-----
Geo/CellComplex.h | 63 ++++++--
Geo/ChainComplex.cpp | 137 +++++++++++------
Geo/ChainComplex.h | 37 +++--
contrib/kbipack/gmp_matrix.c | 21 ++-
utils/misc/celldriver.cpp | 52 +++++--
6 files changed, 479 insertions(+), 116 deletions(-)
diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index b7da2fc688..24616c6e70 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -11,10 +11,69 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su
_domain = domain;
_subdomain = subdomain;
+ bool duplicate = false;
+ for(unsigned int j=0; j < _domain.size(); j++){
+ if(_domain.at(j)->dim() == 3){
+ std::list<GFace*> faces = _domain.at(j)->faces();
+ for(std::list<GFace*>::iterator it = faces.begin(); it != faces.end(); it++){
+ GFace* face = *it;
+ duplicate = false;
+ for(std::vector<GEntity*>::iterator itb = _boundary.begin(); itb != _boundary.end(); itb++){
+ GEntity* entity = *itb;
+ if(face->tag() == entity->tag()){
+ _boundary.erase(itb);
+ duplicate = true;
+ break;
+ }
+ }
+ if(!duplicate) _boundary.push_back(face);
+ }
+ }
+ else if(_domain.at(j)->dim() == 2){
+ std::list<GEdge*> edges = _domain.at(j)->edges();
+
+ for(std::list<GEdge*>::iterator it = edges.begin(); it != edges.end(); it++){
+ GEdge* edge = *it;
+ duplicate = false;
+ std::vector<GEntity*>::iterator erase;
+ for(std::vector<GEntity*>::iterator itb = _boundary.begin(); itb != _boundary.end(); itb++){
+ 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++){
+ GVertex* vertex = *it;
+ duplicate = false;
+ for(std::vector<GEntity*>::iterator itb = _boundary.begin(); itb != _boundary.end(); itb++){
+ GEntity* entity = *itb;
+ if(vertex->tag() == entity->tag()){
+ _boundary.erase(itb);
+ duplicate = true;
+ break;
+ }
+ }
+
+ if(!duplicate) _boundary.push_back(vertex);
+ }
+ }
+ }
+
// subdomain need to be inserted first!
- insertCells(true);
- insertCells(false);
+ insertCells(true, true);
+ insertCells(false, true);
+ insertCells(false, false);
int tag = 1;
for(int i = 0; i < 4; i++){
@@ -23,14 +82,20 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su
cell->setTag(tag);
tag++;
}
+
+ _originalCells[i] = _cells[i];
}
+
+
+
}
-void CellComplex::insertCells(bool subdomain){
+void CellComplex::insertCells(bool subdomain, bool boundary){
std::vector<GEntity*> domain;
if(subdomain) domain = _subdomain;
+ else if(boundary) domain = _boundary;
else domain = _domain;
std::vector<int> vertices;
@@ -43,10 +108,10 @@ void CellComplex::insertCells(bool subdomain){
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));
+ _cells[0].insert(new ZeroSimplex(vertex->getNum(), subdomain, boundary));
}
if(domain.at(j)->getMeshElement(i)->getDim() == 3){
- _cells[3].insert(new ThreeSimplex(vertices, subdomain)); // Add volumes
+ _cells[3].insert(new ThreeSimplex(vertices, subdomain, boundary)); // Add volumes
}
for(int k=0; k < domain.at(j)->getMeshElement(i)->getNumFaces(); k++){
@@ -55,7 +120,7 @@ void CellComplex::insertCells(bool subdomain){
vertex = domain.at(j)->getMeshElement(i)->getFace(k).getVertex(l)->getNum();
vertices.push_back(vertex);
}
- _cells[2].insert(new TwoSimplex(vertices, subdomain)); // Add faces
+ _cells[2].insert(new TwoSimplex(vertices, subdomain, boundary)); // Add faces
}
for(int k=0; k < domain.at(j)->getMeshElement(i)->getNumEdges(); k++){
@@ -64,7 +129,7 @@ void CellComplex::insertCells(bool subdomain){
vertex = domain.at(j)->getMeshElement(i)->getEdge(k).getVertex(l)->getNum();
vertices.push_back(vertex);
}
- _cells[1].insert(new OneSimplex(vertices, subdomain)); // Add edges
+ _cells[1].insert(new OneSimplex(vertices, subdomain, boundary)); // Add edges
}
}
@@ -91,11 +156,19 @@ int Simplex::kappa(Cell* tau) const{
return value;
}
-std::set<Cell*, Less_Cell>::iterator CellComplex::findCell(int dim, std::vector<int>& vertices){
- 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));
+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));
+ }
}
std::set<Cell*, Less_Cell>::iterator CellComplex::findCell(int dim, int vertex, int dummy){
@@ -197,18 +270,18 @@ 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;
+ //if(deleteCell) delete cell;
}
void CellComplex::removeCellIt(std::set<Cell*, Less_Cell>::iterator cell){
Cell* c = *cell;
int dim = c->getDim();
_cells[dim].erase(cell);
- delete c;
+ //delete c;
}
void CellComplex::removeCellQset(Cell*& cell, std::set<Cell*, Less_Cell>& Qset){
Qset.erase(cell);
- delete cell;
+ //delete cell;
}
void CellComplex::enqueueCellsIt(std::vector< std::set<Cell*, Less_Cell>::iterator >& cells,
@@ -232,6 +305,16 @@ void CellComplex::enqueueCells(std::vector<Cell*>& cells, std::queue<Cell*>& Q,
}
}
+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));
+ if(it == Q.end()){
+ Q.push_back(cells.at(i));
+ }
+ }
+}
+
+
int CellComplex::coreductionMrozek(Cell* generator){
@@ -245,9 +328,7 @@ int CellComplex::coreductionMrozek(Cell* generator){
//removeCell(generator);
std::vector<Cell*> bd_s;
- //std::vector< std::set<Cell*, Less_Cell>::iterator > bd_s;
std::vector<Cell*> cbd_c;
- //std::vector< std::set<Cell*, Less_Cell>::iterator > cbd_c;
Cell* s;
int round = 0;
while( !Q.empty() ){
@@ -256,33 +337,99 @@ int CellComplex::coreductionMrozek(Cell* generator){
s = Q.front();
Q.pop();
removeCellQset(s, Qset);
-
- //bd_s = bdIt(s);
+
bd_s = bd(s);
if( bd_s.size() == 1 && inSameDomain(s, bd_s.at(0)) ){
removeCell(s);
- //cbd_c = cbdIt(*bd_s.at(0));
cbd_c = cbd(bd_s.at(0));
- //enqueueCellsIt(cbd_c, Q, Qset);
enqueueCells(cbd_c, Q, Qset);
- //removeCellIt(bd_s.at(0));
removeCell(bd_s.at(0));
coreductions++;
}
else if(bd_s.empty()){
- //cbd_c = cbdIt(s);
- //enqueueCellsIt(cbd_c, Q, Qset);
cbd_c = cbd(s);
enqueueCells(cbd_c, Q, Qset);
}
}
- printf("Coreduction: %d loops with %d coreductions\n", round, coreductions);
+ //printf("Coreduction: %d loops with %d coreductions\n", round, coreductions);
+ return coreductions;
+}
+int CellComplex::coreductionMrozek2(Cell* generator){
+
+ int coreductions = 0;
+
+ std::list<Cell*> Q;
+
+ Q.push_back(generator);
+
+ std::vector<Cell*> bd_s;
+ std::vector<Cell*> cbd_c;
+ Cell* s;
+ int round = 0;
+ while( !Q.empty() ){
+ round++;
+ //printf("%d ", round);
+ s = Q.front();
+ Q.pop_front();
+ bd_s = bd(s);
+ if( bd_s.size() == 1 && inSameDomain(s, bd_s.at(0)) ){
+ removeCell(s);
+ cbd_c = cbd(bd_s.at(0));
+ enqueueCells(cbd_c, Q);
+ removeCell(bd_s.at(0));
+ coreductions++;
+ }
+ else if(bd_s.empty()){
+ cbd_c = cbd(s);
+ 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, bool deleteCells){
if(dim < 0 || dim > 2) return 0;
@@ -413,13 +560,14 @@ void CellComplex::coreduceComplex(int generatorDim, std::set<Cell*, Less_Cell>&
}
void CellComplex::coreduceComplex(int generatorDim){
-
+
std::set<Cell*, Less_Cell> generatorCells[4];
printf("Cell complex before coreduction: %d volumes, %d faces, %d edges and %d vertices.\n",
getSize(3), getSize(2), getSize(1), getSize(0));
for(int i = 0; i < 4; i++){
while (getSize(i) != 0){
+ //if(generatorDim == i || i == generatorDim+1) break;
citer cit = firstCell(i);
Cell* cell = *cit;
while(!cell->inSubdomain() && cit != lastCell(i)){
@@ -428,9 +576,11 @@ void CellComplex::coreduceComplex(int generatorDim){
}
generatorCells[i].insert(cell);
removeCell(cell, false);
- coreduction(0);
- coreduction(1);
- coreduction(2);
+
+ //coreduction(0);
+ //coreduction(1);
+ //coreduction(2);
+ coreductionMrozek(cell);
}
if(generatorDim == i) break;
@@ -438,7 +588,7 @@ void CellComplex::coreduceComplex(int generatorDim){
for(int i = 0; i < 4; i++){
for(citer cit = generatorCells[i].begin(); cit != generatorCells[i].end(); cit++){
Cell* cell = *cit;
- if(!cell->inSubdomain()) _cells[i].insert(cell);
+ //if(!cell->inSubdomain()) _cells[i].insert(cell);
}
}
printf("Cell complex after coreduction: %d volumes, %d faces, %d edges and %d vertices.\n",
@@ -447,6 +597,55 @@ void CellComplex::coreduceComplex(int generatorDim){
return;
}
+void CellComplex::repairComplex(){
+
+ for(int i = 3; i > 0; i--){
+
+ for(citer cit = firstCell(i); cit != lastCell(i); cit++){
+ Cell* cell = *cit;
+ std::vector<int> vertices = cell->getVertexVector();
+
+ for(int j=0; j < vertices.size(); j++){
+ std::vector<int> bVertices;
+
+ for(int k=0; k < vertices.size(); k++){
+ if (k !=j ) bVertices.push_back(vertices.at(k));
+ }
+ citer cit2 = findCell(i-1, bVertices, true);
+ Cell* cell2 = *cit2;
+ _cells[i-1].insert(cell2);
+
+ }
+ }
+
+ }
+/*
+ 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++;
+ }
+ }
+ */
+ return;
+}
+
+void CellComplex::swapSubdomain(){
+
+ for(int i = 0; i < 4; i++){
+ for(citer cit = firstCell(i); cit != lastCell(i); cit++){
+ Cell* cell = *cit;
+ if(cell->onDomainBoundary() && cell->inSubdomain()) cell->setInSubdomain(false);
+ else if(cell->onDomainBoundary() && !cell->inSubdomain()) cell->setInSubdomain(true);
+ }
+ }
+
+ return;
+}
+
+
void CellComplex::getDomainVertices(std::set<MVertex*, Less_MVertex>& domainVertices, bool subdomain){
std::vector<GEntity*> domain;
@@ -531,25 +730,39 @@ int CellComplex::writeComplexMSH(const std::string &name){
fprintf(fp, "$Elements\n");
fprintf(fp, "%d\n", _cells[0].size() + _cells[1].size() + _cells[2].size() + _cells[3].size());
+
+ int physical = 0;
for(citer cit = firstCell(0); cit != lastCell(0); cit++) {
Cell* vertex = *cit;
- fprintf(fp, "%d %d %d %d %d %d %d\n", vertex->getTag(), 15, 3, 0, 0, 0, vertex->getVertex(0));
+ if(vertex->inSubdomain()) physical = 3;
+ else if(vertex->onDomainBoundary()) physical = 2;
+ else physical = 1;
+ fprintf(fp, "%d %d %d %d %d %d %d\n", vertex->getTag(), 15, 3, 0, 0, physical, vertex->getVertex(0));
}
for(citer cit = firstCell(1); cit != lastCell(1); cit++) {
Cell* edge = *cit;
- fprintf(fp, "%d %d %d %d %d %d %d %d\n", edge->getTag(), 1, 3, 0, 0, 0, edge->getVertex(0), edge->getVertex(1));
+ if(edge->inSubdomain()) physical = 3;
+ else if(edge->onDomainBoundary()) physical = 2;
+ 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));
}
for(citer cit = firstCell(2); cit != lastCell(2); cit++) {
Cell* face = *cit;
- fprintf(fp, "%d %d %d %d %d %d %d %d %d\n", face->getTag(), 2, 3, 0, 0, 0, face->getVertex(0), face->getVertex(1), face->getVertex(2));
+ if(face->inSubdomain()) physical = 3;
+ else if(face->onDomainBoundary()) physical = 2;
+ 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));
}
for(citer cit = firstCell(3); cit != lastCell(3); cit++) {
Cell* volume = *cit;
- fprintf(fp, "%d %d %d %d %d %d %d %d %d %d\n", volume->getTag(), 4, 3, 0, 0, 0, volume->getVertex(0), volume->getVertex(1), volume->getVertex(2), volume->getVertex(3));
+ if(volume->inSubdomain()) physical = 3;
+ else if(volume->onDomainBoundary()) physical = 2;
+ 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));
}
fprintf(fp, "$EndElements\n");
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index ee70c278c8..587bfd68bd 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -40,6 +40,9 @@ class Cell
// used in relative homology computation
bool _inSubdomain;
+ // whether this cell belongs to the boundary of a cell complex
+ bool _onDomainBoundary;
+
int _tag;
public:
@@ -73,6 +76,10 @@ class Cell
virtual void setCbdSize(int size) { _cbdSize = size; }
virtual bool inSubdomain() const { return _inSubdomain; }
+ virtual void setInSubdomain(bool subdomain) { _inSubdomain = subdomain; }
+
+ virtual bool onDomainBoundary() const { return _onDomainBoundary; }
+ virtual void setOnDomainBoundary(bool domainboundary) { _onDomainBoundary = domainboundary; }
// print cell vertices
virtual void printCell() const = 0;
@@ -138,7 +145,7 @@ class ZeroSimplex : public Simplex
double _x, _y, _z;
public:
- ZeroSimplex(int vertex, bool subdomain=false, double x=0, double y=0, double z=0) : Simplex(){
+ ZeroSimplex(int vertex, bool subdomain=false, bool boundary=false, double x=0, double y=0, double z=0) : Simplex(){
_v = vertex;
_dim = 0;
_bdMaxSize = 0;
@@ -147,8 +154,8 @@ class ZeroSimplex : public Simplex
_y = y;
_z = z;
_inSubdomain = subdomain;
+ _onDomainBoundary = boundary;
}
-
virtual ~ZeroSimplex(){}
virtual int getDim() const { return 0; }
@@ -175,7 +182,7 @@ class OneSimplex : public Simplex
int _v[2];
public:
- OneSimplex(std::vector<int> vertices, bool subdomain=false) : Simplex(){
+ OneSimplex(std::vector<int> vertices, bool subdomain=false, bool boundary=false) : Simplex(){
sort(vertices.begin(), vertices.end());
_v[0] = vertices.at(0);
_v[1] = vertices.at(1);
@@ -183,7 +190,10 @@ class OneSimplex : public Simplex
_bdMaxSize = 2;
_cbdMaxSize = 1000;
_inSubdomain = subdomain;
+ _onDomainBoundary = boundary;
}
+
+
// constructor for the dummy one simplex
// used to find another definite one simplex from the cell complex
@@ -215,7 +225,7 @@ class TwoSimplex : public Simplex
int _v[3];
public:
- TwoSimplex(std::vector<int> vertices, bool subdomain=false) : Simplex(){
+ TwoSimplex(std::vector<int> vertices, bool subdomain=false, bool boundary=false) : Simplex(){
sort(vertices.begin(), vertices.end());
_v[0] = vertices.at(0);
_v[1] = vertices.at(1);
@@ -224,6 +234,7 @@ class TwoSimplex : public Simplex
_bdMaxSize = 3;
_cbdMaxSize = 2;
_inSubdomain = subdomain;
+ _onDomainBoundary = boundary;
}
// constructor for the dummy one simplex
TwoSimplex(int vertex, int dummy){
@@ -254,7 +265,7 @@ class ThreeSimplex : public Simplex
int _v[4];
public:
- ThreeSimplex(std::vector<int> vertices, bool subdomain=false) : Simplex(){
+ ThreeSimplex(std::vector<int> vertices, bool subdomain=false, bool boundary=false) : Simplex(){
sort(vertices.begin(), vertices.end());
_v[0] = vertices.at(0);
_v[1] = vertices.at(1);
@@ -264,6 +275,7 @@ class ThreeSimplex : public Simplex
_bdMaxSize = 4;
_cbdMaxSize = 0;
_inSubdomain = subdomain;
+ _onDomainBoundary = boundary;
}
// constructor for the dummy one simplex
ThreeSimplex(int vertex, int dummy){
@@ -311,6 +323,21 @@ class Less_Cell{
}
};
+
+class Equal_Cell{
+ public:
+ bool operator()(const Cell* c1, const Cell* c2) const {
+
+ if(c1->getNumVertices() != c2->getNumVertices()) return false;
+
+ for(int i=0; i < c1->getNumVertices();i++){
+ if(c1->getVertex(i) != c2->getVertex(i)) return false;
+ }
+ return true;
+ }
+};
+
+
// Ordering for the finite element mesh vertices.
class Less_MVertex{
public:
@@ -330,10 +357,14 @@ class CellComplex
// used in relative homology computation, may be empty
std::vector<GEntity*> _subdomain;
+ std::vector<GEntity*> _boundary;
+
// 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];
+
public:
// iterator for the cells of same dimension
typedef std::set<Cell*, Less_Cell>::iterator citer;
@@ -342,6 +373,7 @@ class CellComplex
// enqueue cells in queue if they are not there already
virtual void enqueueCells(std::vector<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 enqueueCellsIt(std::vector< std::set<Cell*, Less_Cell>::iterator >& cells,
std::queue<Cell*>& Q, std::set<Cell*, Less_Cell>& Qset);
@@ -349,7 +381,7 @@ class CellComplex
virtual void removeCellQset(Cell*& cell, std::set<Cell*, Less_Cell>& Qset);
// insert cells into this cell complex
- virtual void insertCells(bool subdomain);
+ virtual void insertCells(bool subdomain, bool boundary);
public:
CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> subdomain, std::set<Cell*, Less_Cell> cells ) {
@@ -376,7 +408,7 @@ class CellComplex
virtual citer lastCell(int dim) {return _cells[dim].end(); }
// find a cell in this cell complex
- virtual std::set<Cell*, Less_Cell>::iterator findCell(int dim, std::vector<int>& vertices);
+ virtual std::set<Cell*, Less_Cell>::iterator findCell(int dim, std::vector<int>& vertices, bool original=false);
virtual std::set<Cell*, Less_Cell>::iterator findCell(int dim, int vertex, int dummy=0);
// kappa for two cells of this cell complex
@@ -391,7 +423,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=true);
+ virtual void removeCell(Cell* cell, bool deleteCell=false);
virtual void removeCellIt(std::set<Cell*, Less_Cell>::iterator cell);
virtual void insertCell(Cell* cell);
@@ -402,12 +434,12 @@ 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=true);
+ virtual int coreduction(int dim, bool deleteCells=false);
// 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=true);
+ virtual int reduction(int dim, bool deleteCells=false);
// useful functions for (co)reduction of cell complex
virtual void reduceComplex();
@@ -421,11 +453,20 @@ class CellComplex
// queued coreduction presented in Mrozek's paper
// slower, but produces cleaner result
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]; }
+
+ // add every volume, face and edge its missing boundary cells
+ virtual void repairComplex();
+ // change non-subdomain cells to be in subdomain, subdomain cells to not to be in subdomain
+ virtual void swapSubdomain();
// print the vertices of cells of certain dimension
virtual void printComplex(int dim, bool subcomplex);
- // write this cell complex in .msh format
+ // write this cell complex in 2.0 MSH ASCII format
virtual int writeComplexMSH(const std::string &name);
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index 7bb70d2e6a..524e6235f8 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -11,6 +11,7 @@
ChainComplex::ChainComplex(CellComplex* cellComplex){
+ _dim = 0;
_HMatrix[0] = NULL;
_kerH[0] = NULL;
_codH[0] = NULL;
@@ -21,7 +22,10 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){
for(int dim = 1; dim < 4; dim++){
unsigned int cols = cellComplex->getSize(dim);
unsigned int rows = cellComplex->getSize(dim-1);
+
+
+ // ignore subdomain cells
for(std::set<Cell*, Less_Cell>::iterator cit = cellComplex->firstCell(dim); cit != cellComplex->lastCell(dim); cit++){
Cell* cell = *cit;
if(cell->inSubdomain()) cols--;
@@ -33,11 +37,15 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){
if( cols == 0 ){
+ //_HMatrix[dim] = create_gmp_matrix_zero(rows, 1);
_HMatrix[dim] = NULL;
}
else if( rows == 0){
_HMatrix[dim] = create_gmp_matrix_zero(1, cols);
+ //_HMatrix[dim] = NULL;
}
+
+
else{
long int elems[rows*cols];
@@ -61,20 +69,24 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){
}
_HMatrix[dim] = create_gmp_matrix_int(rows, cols, elems);
}
-
+
_kerH[dim] = NULL;
_codH[dim] = NULL;
_JMatrix[dim] = NULL;
_QMatrix[dim] = NULL;
_Hbasis[dim] = NULL;
+
+ if(cellComplex->getSize(dim) != 0) _dim = dim;
+
}
+
return;
}
void ChainComplex::KerCod(int dim){
- if(dim < 1 || dim > 3 || _HMatrix[dim] == NULL) return;
+ if(dim < 0 || dim > 3 || _HMatrix[dim] == NULL) return;
gmp_matrix* HMatrix = copy_gmp_matrix(_HMatrix[dim], 1, 1,
gmp_matrix_rows(_HMatrix[dim]), gmp_matrix_cols(_HMatrix[dim]));
@@ -113,30 +125,34 @@ void ChainComplex::KerCod(int dim){
}
//j:B_(k+1)->Z_k
-void ChainComplex::Inclusion(int dim){
+void ChainComplex::Inclusion(int lowDim, int highDim){
+
+ if(getKerHMatrix(lowDim) == NULL || getCodHMatrix(highDim) == NULL || abs(lowDim-highDim) != 1) return;
- if(dim < 1 || dim > 2 || _kerH[dim] == NULL || _codH[dim+1] == NULL) return;
+ gmp_matrix* Zbasis = copy_gmp_matrix(_kerH[lowDim], 1, 1,
+ gmp_matrix_rows(_kerH[lowDim]), gmp_matrix_cols(_kerH[lowDim]));
+ gmp_matrix* Bbasis = copy_gmp_matrix(_codH[highDim], 1, 1,
+ gmp_matrix_rows(_codH[highDim]), gmp_matrix_cols(_codH[highDim]));
- gmp_matrix* Zbasis = copy_gmp_matrix(_kerH[dim], 1, 1,
- gmp_matrix_rows(_kerH[dim]), gmp_matrix_cols(_kerH[dim]));
- gmp_matrix* Bbasis = copy_gmp_matrix(_codH[dim+1], 1, 1,
- gmp_matrix_rows(_codH[dim+1]), gmp_matrix_cols(_codH[dim+1]));
- int rows = gmp_matrix_rows(Zbasis);
- int cols = gmp_matrix_cols(Zbasis);
+ int rows = gmp_matrix_rows(Bbasis);
+ int cols = gmp_matrix_cols(Bbasis);
if(rows < cols) return;
- rows = gmp_matrix_rows(Bbasis);
- cols = gmp_matrix_cols(Bbasis);
+ rows = gmp_matrix_rows(Zbasis);
+ cols = gmp_matrix_cols(Zbasis);
if(rows < cols) return;
+
// A*inv(V) = U*S
gmp_normal_form* normalForm = create_gmp_Smith_normal_form(Zbasis, INVERTED, INVERTED);
mpz_t elem;
mpz_init(elem);
+
for(int i = 1; i <= cols; i++){
+
gmp_matrix_get_elem(elem, i, i, normalForm->canonical);
if(mpz_cmp_si(elem,0) == 0){
destroy_gmp_normal_form(normalForm);
@@ -175,7 +191,7 @@ void ChainComplex::Inclusion(int dim){
gmp_matrix_left_mult(normalForm->right, LB);
- _JMatrix[dim] = LB;
+ setJMatrix(lowDim, LB);
mpz_clear(elem);
mpz_clear(divisor);
@@ -188,7 +204,7 @@ void ChainComplex::Inclusion(int dim){
void ChainComplex::Quotient(int dim){
- if(dim < 1 || dim > 3 || _JMatrix[dim] == NULL) return;
+ if(dim < 0 || dim > 3 || _JMatrix[dim] == NULL) return;
gmp_matrix* JMatrix = copy_gmp_matrix(_JMatrix[dim], 1, 1,
gmp_matrix_rows(_JMatrix[dim]), gmp_matrix_cols(_JMatrix[dim]));
@@ -196,6 +212,11 @@ void ChainComplex::Quotient(int dim){
int cols = gmp_matrix_cols(JMatrix);
gmp_normal_form* normalForm = create_gmp_Smith_normal_form(JMatrix, NOT_INVERTED, NOT_INVERTED);
+
+ //printMatrix(normalForm->left);
+ //printMatrix(normalForm->canonical);
+ //printMatrix(normalForm->right);
+
mpz_t elem;
mpz_init(elem);
@@ -222,28 +243,43 @@ void ChainComplex::Quotient(int dim){
return;
}
-void ChainComplex::computeHomology(){
+void ChainComplex::computeHomology(bool dual){
+
+ int lowDim = 0;
+ int highDim = 0;
for(int i=0; i < 4; i++){
-
- KerCod(i+1);
- // 1) no edges, but zero cells
- if(i == 0 && gmp_matrix_cols(getHMatrix(0)) > 0 && getHMatrix(1) == NULL) {
- _Hbasis[i] = create_gmp_matrix_identity(gmp_matrix_cols(getHMatrix(0)));
+ if(dual){
+ lowDim = getDim()+1-i;
+ highDim = getDim()+1-(i+1);
+ //KerCod(lowDim);
+ }
+ else{
+ lowDim = i;
+ highDim = i+1;
+ //KerCod(highDim);
}
+ printf("Homology computation process: step %d of 4 \n", i+1 );
+
+ KerCod(highDim);
+
+ // 1) no edges, but zero cells
+ if(lowDim == 0 && gmp_matrix_cols(getHMatrix(lowDim)) > 0 && getHMatrix(highDim) == NULL) {
+ setHbasis( lowDim, create_gmp_matrix_identity(gmp_matrix_cols(getHMatrix(lowDim))) );
+ }
// 2) this dimension is empty
- else if(getHMatrix(i) == NULL && getHMatrix(i+1) == NULL){
- _Hbasis[i] = NULL;
+ else if(getHMatrix(lowDim) == NULL && getHMatrix(highDim) == NULL){
+ setHbasis(lowDim, NULL);
}
// 3) No higher dimension cells -> none of the cycles are boundaries
- else if(getHMatrix(i+1) == NULL){
- _Hbasis[i] = copy_gmp_matrix(getKerHMatrix(i), 1, 1,
- gmp_matrix_rows(getKerHMatrix(i)), gmp_matrix_cols(getKerHMatrix(i)));
+ else if(getHMatrix(highDim) == NULL){
+ setHbasis( lowDim, copy_gmp_matrix(getKerHMatrix(lowDim), 1, 1,
+ 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
@@ -251,39 +287,40 @@ void ChainComplex::computeHomology(){
else {
// 4) No lower dimension cells -> all chains are cycles
- if(getHMatrix(i) == NULL){
- _kerH[i] = create_gmp_matrix_identity(gmp_matrix_rows(getHMatrix(i+1)));
+ if(getHMatrix(lowDim) == NULL){
+ setKerHMatrix(lowDim, create_gmp_matrix_identity(gmp_matrix_rows(getHMatrix(highDim))) );
}
+ Inclusion(lowDim, highDim);
+ Quotient(lowDim);
- Inclusion(i);
- Quotient(i);
-
- if(getCodHMatrix(i+1) == NULL){
- _Hbasis[i] = copy_gmp_matrix(getKerHMatrix(i), 1, 1,
- gmp_matrix_rows(getKerHMatrix(i)), gmp_matrix_cols(getKerHMatrix(i)));
+ if(getCodHMatrix(highDim) == NULL){
+ setHbasis(lowDim, copy_gmp_matrix(getKerHMatrix(lowDim), 1, 1,
+ gmp_matrix_rows(getKerHMatrix(lowDim)),
+ gmp_matrix_cols(getKerHMatrix(lowDim))) );
}
- else if(getJMatrix(i) == NULL || getQMatrix(i) == NULL){
- _Hbasis[i] = NULL;
+ else if(getJMatrix(lowDim) == NULL || getQMatrix(lowDim) == NULL){
+ setHbasis(lowDim, NULL);
}
else{
- _Hbasis[i] = copy_gmp_matrix(getKerHMatrix(i), 1, 1,
- gmp_matrix_rows(getKerHMatrix(i)), gmp_matrix_cols(getKerHMatrix(i)));
+ setHbasis(lowDim, copy_gmp_matrix(getKerHMatrix(lowDim), 1, 1,
+ gmp_matrix_rows(getKerHMatrix(lowDim)),
+ gmp_matrix_cols(getKerHMatrix(lowDim))) );
- gmp_matrix_right_mult(_Hbasis[i], getQMatrix(i));
+ gmp_matrix_right_mult(getHbasis(lowDim), getQMatrix(lowDim));
}
}
- destroy_gmp_matrix(_kerH[i]);
- destroy_gmp_matrix(_codH[i]);
- destroy_gmp_matrix(_JMatrix[i]);
- destroy_gmp_matrix(_QMatrix[i]);
+ destroy_gmp_matrix(getKerHMatrix(lowDim));
+ destroy_gmp_matrix(getCodHMatrix(lowDim));
+ destroy_gmp_matrix(getJMatrix(lowDim));
+ destroy_gmp_matrix(getQMatrix(lowDim));
- _kerH[i] = NULL;
- _codH[i] = NULL;
- _JMatrix[i] = NULL;
- _QMatrix[i] = NULL;
+ setKerHMatrix(lowDim, NULL);
+ setCodHMatrix(lowDim, NULL);
+ setJMatrix(lowDim, NULL);
+ setQMatrix(lowDim, NULL);
}
@@ -359,6 +396,8 @@ int Chain::writeChainMSH(const std::string &name){
//_cellComplex->writeComplexMSH(name);
+ if(getSize() == 0) return 0;
+
FILE *fp = fopen(name.c_str(), "a");
if(!fp){
Msg::Error("Unable to open file '%s'", name.c_str());
@@ -378,7 +417,7 @@ int Chain::writeChainMSH(const std::string &name){
fprintf(fp, "%d \n", getSize());
fprintf(fp, "0 \n");
- for(int i = 0; i < _cells.size(); i++){
+ for(int i = 0; i < getSize(); i++){
fprintf(fp, "%d %d \n", getCell(i)->getTag(), getCoeff(i));
diff --git a/Geo/ChainComplex.h b/Geo/ChainComplex.h
index 4e37193107..0ae7756866 100644
--- a/Geo/ChainComplex.h
+++ b/Geo/ChainComplex.h
@@ -49,6 +49,16 @@ class ChainComplex{
// bases for homology groups
gmp_matrix* _Hbasis[4];
+ int _dim;
+
+ // set the matrices
+ virtual void setHMatrix(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 4) _HMatrix[dim] = matrix;}
+ virtual void setKerHMatrix(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 4) _kerH[dim] = matrix;}
+ virtual void setCodHMatrix(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 4) _codH[dim] = matrix;}
+ virtual void setJMatrix(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 4) _JMatrix[dim] = matrix;}
+ virtual void setQMatrix(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 4) _QMatrix[dim] = matrix;}
+ virtual void setHbasis(int dim, gmp_matrix* matrix) { if(dim > -1 && dim < 4) _Hbasis[dim] = matrix;}
+
public:
@@ -64,31 +74,38 @@ class ChainComplex{
_QMatrix[i] = NULL;
_Hbasis[i] = NULL;
}
+ _dim = 0;
}
virtual ~ChainComplex(){}
+ virtual int getDim() { return _dim; }
+
// get the boundary operator matrix dim->dim-1
- virtual gmp_matrix* getHMatrix(int dim) { if(dim > -1 || dim < 4) return _HMatrix[dim]; else return NULL;}
- virtual gmp_matrix* getKerHMatrix(int dim) { if(dim > -1 || dim < 4) return _kerH[dim]; else return NULL;}
- virtual gmp_matrix* getCodHMatrix(int dim) { if(dim > -1 || dim < 4) return _codH[dim]; else return NULL;}
- virtual gmp_matrix* getJMatrix(int dim) { if(dim > -1 || dim < 4) return _JMatrix[dim]; else return NULL;}
- virtual gmp_matrix* getQMatrix(int dim) { if(dim > -1 || dim < 4) return _QMatrix[dim]; else return NULL;}
- virtual gmp_matrix* getHbasis(int dim) { if(dim > -1 || dim < 4) return _Hbasis[dim]; else return NULL;}
+ virtual gmp_matrix* getHMatrix(int dim) { if(dim > -1 && dim < 4) return _HMatrix[dim]; else return NULL;}
+ virtual gmp_matrix* getKerHMatrix(int dim) { if(dim > -1 && dim < 4) return _kerH[dim]; else return NULL;}
+ virtual gmp_matrix* getCodHMatrix(int dim) { if(dim > -1 && dim < 4) return _codH[dim]; else return NULL;}
+ virtual gmp_matrix* getJMatrix(int dim) { if(dim > -1 && dim < 4) return _JMatrix[dim]; else return NULL;}
+ virtual gmp_matrix* getQMatrix(int dim) { if(dim > -1 && dim < 4) return _QMatrix[dim]; else return NULL;}
+ virtual gmp_matrix* getHbasis(int dim) { if(dim > -1 && dim < 4) return _Hbasis[dim]; else return NULL;}
// Compute basis for kernel and codomain of boundary operator matrix of dimension dim (ie. ker(h_dim) and cod(h_dim) )
virtual void KerCod(int dim);
// Compute matrix representation J for inclusion relation from dim-cells who are boundary of dim+1-cells
// to cycles of dim-cells (ie. j: cod(h_(dim+1)) -> ker(h_dim) )
- virtual void Inclusion(int dim);
+ virtual void Inclusion(int lowDim, int highDim);
// Compute quotient problem for given inclusion relation j to find representatives of homology groups
// and possible torsion coeffcients
virtual void Quotient(int dim);
+ // transpose the boundary operator matrices, these are boundary operator matrices for the dual mesh
+ virtual void transposeHMatrices() { for(int i = 0; i < 4; i++) gmp_matrix_transp(_HMatrix[i]); }
+ virtual void transposeHMatrix(int dim) { if(dim > -1 && dim < 4) gmp_matrix_transp(_HMatrix[dim]); }
+
// Compute bases for the homology groups of this chain complex
- virtual void computeHomology();
+ virtual void computeHomology(bool dual=false);
virtual std::vector<int> getCoeffVector(int dim, int chainNumber);
- virtual int getBasisSize(int dim) { if(dim > -1 || dim < 4) return gmp_matrix_cols(_Hbasis[dim]); else return 0; }
+ virtual int getBasisSize(int dim) { if(dim > -1 && dim < 4) return gmp_matrix_cols(_Hbasis[dim]); else return 0; }
virtual int printMatrix(gmp_matrix* matrix){
printf("%d rows and %d columns\n", gmp_matrix_rows(matrix), gmp_matrix_cols(matrix));
@@ -129,7 +146,7 @@ class Chain{
virtual std::string getName() { return _name; }
virtual void setName(std::string name) { _name=name; }
- // append this chain to a 2.0 .msh file as $ElementData
+ // append this chain to a 2.0 MSH ASCII file as $ElementData
virtual int writeChainMSH(const std::string &name);
};
diff --git a/contrib/kbipack/gmp_matrix.c b/contrib/kbipack/gmp_matrix.c
index 47bf32dd82..531466abd1 100755
--- a/contrib/kbipack/gmp_matrix.c
+++ b/contrib/kbipack/gmp_matrix.c
@@ -22,7 +22,7 @@
P.O.Box 692, FIN-33101 Tampere, Finland
saku.suuriniemi@tut.fi
- $Id: gmp_matrix.c,v 1.4 2009-04-21 07:06:22 matti Exp $
+ $Id: gmp_matrix.c,v 1.5 2009-05-05 11:35:01 matti Exp $
*/
@@ -584,6 +584,23 @@ gmp_matrix_transp(gmp_matrix * M)
return EXIT_FAILURE;
}
+ if(rows == 1){
+ for(i = 1; i <= cols; i++)
+ {
+ mpz_init_set(new_storage[i-1],
+ M-> storage[i-1]);
+ mpz_clear(M-> storage[i-1]);
+ }
+ }
+ else if(cols == 1){
+ for(i = 1; i <= rows; i++)
+ {
+ mpz_init_set(new_storage[i-1],
+ M-> storage[i-1]);
+ mpz_clear(M-> storage[i-1]);
+ }
+ }
+ else{
for(i = 1; i <= rows; i++)
{
for(j = 1; j <= cols; j++)
@@ -593,6 +610,8 @@ gmp_matrix_transp(gmp_matrix * M)
mpz_clear(M-> storage[(i-1)+(j-1)*rows]);
}
}
+ }
+
free(M->storage);
M -> storage = new_storage;
diff --git a/utils/misc/celldriver.cpp b/utils/misc/celldriver.cpp
index a515b2c819..d54683dfd4 100755
--- a/utils/misc/celldriver.cpp
+++ b/utils/misc/celldriver.cpp
@@ -6,8 +6,8 @@
// Then compile this driver with "g++ celldriver.cpp -lGmsh -llapack -lblas -lgmp"
//
//
-// This program creates .msh file chains.msh which represents the generators of
-// an two port transformer model's homology groups.
+// This program creates .msh file chains.msh which represents the generators and
+// the cuts of an two port transformer model's homology groups.
#include <stdio.h>
@@ -33,20 +33,19 @@ int main(int argc, char **argv)
GModel *m = new GModel();
m->readGEO("transformer.geo");
m->mesh(3);
- m->writeMSH("transformer.msh");
printf("This model has: %d GRegions, %d GFaces, %d GEdges and %d GVertices. \n" , m->getNumRegions(), m->getNumFaces(), m->getNumEdges(), m->getNumVertices());
std::vector<GEntity*> domain;
std::vector<GEntity*> subdomain;
- // whole model
+ // whole model the domain
for(GModel::riter rit = m->firstRegion(); rit != m->lastRegion(); rit++){
GEntity *r = *rit;
domain.push_back(r);
}
- // the ports
+ // the ports as subdomain
GModel::fiter sub = m->firstFace();
GEntity *s = *sub;
subdomain.push_back(s);
@@ -60,6 +59,9 @@ int main(int argc, char **argv)
// create a cell complex
CellComplex* complex = new CellComplex(domain, subdomain);
+ // write the cell complex to a .msh file
+ complex->writeComplexMSH("chains.msh");
+
printf("Cell complex of this model has: %d volumes, %d faces, %d edges and %d vertices\n",
complex->getSize(3), complex->getSize(2), complex->getSize(1), complex->getSize(0));
@@ -71,9 +73,6 @@ int main(int argc, char **argv)
// compute the homology using the boundary operator matrices
chains->computeHomology();
- // write the reduced cell complex to a .msh file
- complex->writeComplexMSH("chains.msh");
-
// Append the homology generators to the .msh file as post-processing views.
// To visualize the result, open chains.msh with Gmsh GUI and deselect all mesh
// entities from Tools->Options->Mesh->Visibility.
@@ -93,8 +92,43 @@ int main(int argc, char **argv)
}
}
-
delete chains;
+
+ // restore the cell complex to its prereduced state
+ complex->restoreComplex();
+
+ // reduce the complex by coreduction, this removes all vertices, hence 0 as an argument
+ complex->coreduceComplex(0);
+
+ // create a chain complex of a cell complex (construct boundary operator matrices)
+ ChainComplex* dualChains = new ChainComplex(complex);
+
+ // transpose the boundary operator matrices,
+ // these are the boundary operator matrices for the dual chain complex
+ dualChains->transposeHMatrices();
+
+ // compute the homology of the dual complex
+ dualChains->computeHomology(true);
+
+ // Append the duals of the cuts of the homology generators to the .msh file.
+ for(int j = 3; j > -1; j--){
+ for(int i = 1; i <= dualChains->getBasisSize(j); i++){
+
+ std::string generator;
+ std::string dimension;
+ convert(i, generator);
+ 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->writeChainMSH("chains.msh");
+ delete chain;
+ }
+ }
+
+ delete dualChains;
+
+
delete complex;
delete m;
GmshFinalize();
--
GitLab