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

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.
parent 14ca4b5b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Geo/CellComplex.cpp
+ 249
− 36
View file @ a804c7a4
......@@ -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");
......
Geo/CellComplex.h
+ 52
− 11
View file @ a804c7a4
......@@ -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);
......
Geo/ChainComplex.cpp
+ 88
− 49
View file @ a804c7a4
......@@ -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));
......
Geo/ChainComplex.h
+ 27
− 10
View file @ a804c7a4
......@@ -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);
};
......
contrib/kbipack/gmp_matrix.c
+ 20
− 1
View file @ a804c7a4
......@@ -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;
......
utils/misc/celldriver.cpp
+ 43
− 9
View file @ a804c7a4
......@@ -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();
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment