From 71295301f346f0b1048c1b8e6efca144d1dc0d73 Mon Sep 17 00:00:00 2001
From: Matti Pellika <matti.pellikka@tut.fi>
Date: Mon, 24 Aug 2009 12:46:56 +0000
Subject: [PATCH] Improved reduction algorithm when generators of highest
dimension are present. Added Betti number computation using coreduction to
the UI. Added support to quadrangles (not yet to prisms, hexaherdons or
pyramids due to unresolved issues). Bugfixes.
---
Geo/CellComplex.cpp | 169 ++++++++++++++--------
Geo/CellComplex.h | 257 ++++++++++++++++++++++++++++++---
Geo/Homology.cpp | 38 ++++-
Geo/Homology.h | 5 +-
Plugin/HomologyComputation.cpp | 21 ++-
5 files changed, 394 insertions(+), 96 deletions(-)
diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 46eb4b36c2..2eb31a86ac 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -103,10 +103,10 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su
for(int i = 0; i < 4; i++){
for(citer cit = firstCell(i); cit != lastCell(i); cit++){
Cell* cell = *cit;
- cell->setTag(tag);
+ //cell->setTag(tag);
tag++;
}
- _cells2[i] = _cells[i];
+ //_cells2[i] = _cells[i];
_betti[i] = 0;
if(getSize(i) > _dim) _dim = i;
}
@@ -127,7 +127,7 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){
std::vector<MVertex*> vertices;
std::pair<citer, bool> insertInfo;
-
+ int tag = 1;
// add highest dimensional cells
for(unsigned int j=0; j < domain.size(); j++) {
@@ -140,11 +140,34 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){
}
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, 0, subdomain, boundary);
- else if(dim == 2) cell = new TwoSimplex(vertices, 0, subdomain, boundary);
- else if(dim == 1) cell = new OneSimplex(vertices, 0, subdomain, boundary);
- else cell = new ZeroSimplex(vertices, 0, subdomain, boundary);
+ // simplex types
+ if(type == MSH_LIN_2 || type == MSH_TRI_3 || type == MSH_TET_4 || type == MSH_LIN_3 || type == MSH_TRI_6 ||
+ type == MSH_TET_10 || type == MSH_PNT || type == MSH_TRI_9 || type == MSH_TRI_10 || type == MSH_TRI_12 ||
+ type == MSH_TRI_15 || type == MSH_TRI_15I || type == MSH_TRI_21 || type == MSH_LIN_4 ||
+ type == MSH_LIN_5 || type == MSH_LIN_6 || type == MSH_TET_20 || type == MSH_TET_35 || type == MSH_TET_56
+ || type == MSH_TET_34 || type == MSH_TET_52 ){
+ if(dim == 3) cell = new ThreeSimplex(vertices, tag, subdomain, boundary);
+ else if(dim == 2) cell = new TwoSimplex(vertices, tag, subdomain, boundary);
+ else if(dim == 1) cell = new OneSimplex(vertices, tag, subdomain, boundary);
+ else cell = new ZeroSimplex(vertices, tag, subdomain, boundary);
+ }
+ else if(type == MSH_QUA_4 || type == MSH_QUA_8 || type == MSH_QUA_9){
+ cell = new CQuadrangle(vertices, tag, subdomain, boundary);
+ }
+ /*
+ else if(type == MSH_HEX_8 || type == MSH_HEX_27 || type == MSH_HEX_20){
+ cell = new CHexahedron(vertices, tag, subdomain, boundary);
+ }
+ else if(type == MSH_PRI_6 || type == MSH_PRI_18 || type == MSH_PRI_15){
+ cell = new CPrism(vertices, tag, subdomain, boundary);
+ }
+ else if(type == MSH_PYR_5 || type == MSH_PYR_14 || type == MSH_PYR_13){
+ cell = new CPyramid(vertices, tag, subdomain, boundary);
+ }*/
+ else printf("Error: mesh element %d not implemented yet! \n", type);
+ tag++;
insertInfo = _cells[dim].insert(cell);
if(!insertInfo.second) delete cell;
}
@@ -158,9 +181,14 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){
for(int i = 0; i < cell->getNumFacets(); i++){
cell->getFacetVertices(i, vertices);
Cell* newCell;
- if(dim == 3) newCell = new TwoSimplex(vertices, 0, subdomain, boundary);
- else if(dim == 2) newCell = new OneSimplex(vertices, 0, subdomain, boundary);
- else if(dim == 1) newCell = new ZeroSimplex(vertices, 0, subdomain, boundary);
+ if(dim == 3){
+ if(vertices.size() == 3) newCell = new TwoSimplex(vertices, tag, subdomain, boundary);
+ else if(vertices.size() == 4) newCell = new CQuadrangle(vertices, tag, subdomain, boundary);
+ else printf("Error: invalid face! \n");
+ }
+ else if(dim == 2) newCell = new OneSimplex(vertices, tag, subdomain, boundary);
+ else if(dim == 1) newCell = new ZeroSimplex(vertices, tag, subdomain, boundary);
+ tag++;
insertInfo = _cells[dim-1].insert(newCell);
if(!insertInfo.second){
delete newCell;
@@ -186,7 +214,7 @@ void CellComplex::insertCell(Cell* cell){
_cells[cell->getDim()].insert(cell);
}
-
+/*
int Simplex::kappa(Cell* tau) const{
for(int i=0; i < tau->getNumVertices(); i++){
if( !(this->hasVertex(tau->getVertex(i)->getNum())) ) return 0;
@@ -202,11 +230,11 @@ int Simplex::kappa(Cell* tau) const{
return value;
}
+*/
-/*
-int Simplex::kappa(Cell* tau) const{
+int Cell::kappa(Cell* tau) const{
for(int i=0; i < tau->getNumVertices(); i++){
- if( !(this->hasVertex(tau->getVertex(i))) ) return 0;
+ if( !(this->hasVertex(tau->getVertex(i)->getNum())) ) return 0;
}
if(this->getDim() - tau->getDim() != 1) return 0;
@@ -219,6 +247,8 @@ int Simplex::kappa(Cell* tau) const{
this->getFacetVertices(i, v);
value = -1;
+ if(v.size() != vTau.size()) printf("Error: invalid facet!");
+
do {
value = value*-1;
if(v == vTau) return value;
@@ -238,11 +268,12 @@ int Simplex::kappa(Cell* tau) const{
return 0;
}
-*/
-/*
+
+
int OneSimplex::kappa(Cell* tau) const{
+
for(int i=0; i < tau->getNumVertices(); i++){
- if( !(this->hasVertex(tau->getVertex(i))) ) return 0;
+ if( !(this->hasVertex(tau->getVertex(i)->getNum())) ) return 0;
}
if(tau->getDim() != 0) return 0;
@@ -251,9 +282,9 @@ int OneSimplex::kappa(Cell* tau) const{
else return 1;
}
-*/
-void CellComplex::removeCell(Cell* cell){
+
+void CellComplex::removeCell(Cell* cell, bool other){
std::list< std::pair< int, Cell*> > coboundary = cell->getOrientedCoboundary();
std::list< std::pair< int, Cell*> > boundary = cell->getOrientedBoundary();
@@ -264,7 +295,7 @@ void CellComplex::removeCell(Cell* cell){
//for(std::list<Cell*>::iterator it = coboundary.begin(); it != coboundary.end(); it++){
Cell* cbdCell = (*it).second;
//Cell* cbdCell = *it;
- cbdCell->removeBoundaryCell(cell);
+ cbdCell->removeBoundaryCell(cell, other);
}
@@ -272,7 +303,7 @@ void CellComplex::removeCell(Cell* cell){
//for(std::list<Cell*>::iterator it = boundary.begin(); it != boundary.end(); it++){
Cell* bdCell = (*it).second;
//Cell* bdCell = *it;
- bdCell->removeCoboundaryCell(cell);
+ bdCell->removeCoboundaryCell(cell, other);
}
_cells[cell->getDim()].erase(cell);
@@ -317,14 +348,18 @@ int CellComplex::coreduction(Cell* generator){
s = Q.front();
Q.pop();
removeCellQset(s, Qset);
-
bd_s = s->getBoundary();
+
if( bd_s.size() == 1 && inSameDomain(s, bd_s.front()) ){
removeCell(s);
cbd_c = bd_s.front()->getCoboundary();
enqueueCells(cbd_c, Q, Qset);
removeCell(bd_s.front());
coreductions++;
+
+ _trash.push_back(s);
+ _trash.push_back(bd_s.front());
+
}
else if(bd_s.empty()){
cbd_c = s->getCoboundary();
@@ -344,24 +379,32 @@ int CellComplex::reduction(int dim){
bool reduced = true;
while (reduced){
+
reduced = false;
- for(citer cit = firstCell(dim-1); cit != lastCell(dim-1); cit++){
+ citer cit = firstCell(dim-1);
+ while(cit != lastCell(dim-1)){
+
Cell* cell = *cit;
cbd_c = cell->getCoboundary();
if( cbd_c.size() == 1 && inSameDomain(cell, cbd_c.front()) ){
-
+
++cit;
removeCell(cbd_c.front());
removeCell(cell);
+ _trash.push_back(cell);
+ _trash.push_back(cbd_c.front());
+
count++;
reduced = true;
}
+ if(getSize(dim) == 0 || getSize(dim-1) == 0) break;
+ cit++;
}
-
+
}
-
+
return count;
}
/*
@@ -455,16 +498,14 @@ int CellComplex::reduceComplex(int omit){
printf("Cell complex before reduction: %d volumes, %d faces, %d edges and %d vertices.\n",
getSize(3), getSize(2), getSize(1), getSize(0));
-
+
int count = 0;
for(int i = 3; i > 0; i--) count = count + reduction(i);
-
-
int omitted = 0;
if(omit > getDim()) omit = getDim();
- if(count == 0 && omit > 0){
+ //if(count == 0 && omit > 0){
CellComplex::removeSubdomain();
@@ -472,7 +513,7 @@ int CellComplex::reduceComplex(int omit){
std::set<Cell*, Less_Cell> generatorCells;
while (getSize(getDim()) != 0){
-
+
citer cit = firstCell(getDim());
Cell* cell = *cit;
@@ -494,14 +535,13 @@ int CellComplex::reduceComplex(int omit){
cell->clearCoboundary();
_cells[cell->getDim()].insert(cell);
}
-
- }
+ //}
double t2 = Cpu();
printf("Cell complex after reduction: %d volumes, %d faces, %d edges and %d vertices (%g s).\n",
getSize(3), getSize(2), getSize(1), getSize(0), t2 - t1);
-
+
return 0;
}
@@ -512,6 +552,7 @@ void CellComplex::removeSubdomain(){
for(citer cit = firstCell(i); cit != lastCell(i); cit++){
Cell* cell = *cit;
if(cell->inSubdomain()) {
+ //++cit;
removeCell(cell);
cit = firstCell(i);
}
@@ -543,7 +584,7 @@ int CellComplex::coreduceComplex(int omit){
}
int omitted = 0;
- if(count == 0 && omit > 0){
+ //if(count == 0 && omit > 0){
//if(omit > getDim()) omit = getDim();
//for(int i = 0; i < omit; i++){
std::set<Cell*, Less_Cell> generatorCells;
@@ -551,9 +592,8 @@ int CellComplex::coreduceComplex(int omit){
citer cit = firstCell(0);
Cell* cell = *cit;
generatorCells.insert(cell);
- removeCell(cell);
- coreduceComplex(0);
- //coreduction(cell);
+ removeCell(cell, false);
+ coreduction(cell);
omitted++;
}
@@ -564,7 +604,7 @@ int CellComplex::coreduceComplex(int omit){
_cells[0].insert(cell);
}
- }
+ //}
printf("Cell complex after coreduction: %d volumes, %d faces, %d edges and %d vertices.\n",
getSize(3), getSize(2), getSize(1), getSize(0));
@@ -575,11 +615,10 @@ int CellComplex::coreduceComplex(int omit){
void CellComplex::computeBettiNumbers(){
for(int i = 0; i < 4; i++){
+ _betti[i] = 0;
printf("Betti number computation process: step %d of 4 \n", i+1);
-
-
+
while (getSize(i) != 0){
-
citer cit = firstCell(i);
Cell* cell = *cit;
while(!cell->inSubdomain() && cit != lastCell(i)){
@@ -587,11 +626,9 @@ void CellComplex::computeBettiNumbers(){
cit++;
}
if(!cell->inSubdomain()) _betti[i] = _betti[i] + 1;
- removeCell(cell);
-
+ removeCell(cell, false);
coreduction(cell);
}
-
}
printf("Cell complex Betti numbers: \n b0 = %d \n b1 = %d \n b2 = %d \n b3 = %d \n",
getBettiNumber(0), getBettiNumber(1), getBettiNumber(2), getBettiNumber(3));
@@ -636,6 +673,7 @@ int CellComplex::cocombine(int dim){
Cell* c2 = bd_c.back().second;
removeCell(s);
+ _trash.push_back(s);
cbd_c = c1->getCoboundary();
enqueueCells(cbd_c, Q, Qset);
@@ -697,7 +735,8 @@ int CellComplex::combine(int dim){
Cell* c2 = cbd_c.back().second;
removeCell(s);
-
+ _trash.push_back(s);
+
bd_c = c1->getBoundary();
enqueueCells(bd_c, Q, Qset);
bd_c = c2->getBoundary();
@@ -840,49 +879,53 @@ int CellComplex::writeComplexMSH(const std::string &name){
fprintf(fp, "%d\n", count);
- int physical = 0;
+ int partition = 0;
for(citer cit = firstCell(0); cit != lastCell(0); cit++) {
Cell* vertex = *cit;
- 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->getNum(), 15, 3, 0, 0, physical, vertex->getVertex(0)->getNum());
+ if(vertex->inSubdomain()) partition = 3;
+ else if(vertex->onDomainBoundary()) partition = 2;
+ else partition = 1;
+ fprintf(fp, "%d %d %d %d %d %d %d\n", vertex->getNum(), 15, 3, 0, 0, partition, vertex->getVertex(0)->getNum());
}
std::list< std::pair<int, Cell*> > cells;
for(citer cit = firstCell(1); cit != lastCell(1); cit++) {
Cell* edge = *cit;
- if(edge->inSubdomain()) physical = 3;
- else if(edge->onDomainBoundary()) physical = 2;
- else physical = 1;
+ if(edge->inSubdomain()) partition = 3;
+ else if(edge->onDomainBoundary()) partition = 2;
+ else partition = 1;
cells = edge->getCells();
for(std::list< std::pair<int, Cell*> >::iterator it = cells.begin(); it != cells.end(); it++){
Cell* cell = (*it).second;
- fprintf(fp, "%d %d %d %d %d %d %d %d\n", cell->getNum(), 1, 3, 0, 0, physical, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum());
+ fprintf(fp, "%d %d %d %d %d %d %d %d\n", cell->getNum(), 1, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum());
}
}
for(citer cit = firstCell(2); cit != lastCell(2); cit++) {
Cell* face = *cit;
- if(face->inSubdomain()) physical = 3;
- else if(face->onDomainBoundary()) physical = 2;
- else physical = 1;
+ if(face->inSubdomain()) partition = 3;
+ else if(face->onDomainBoundary()) partition = 2;
+ else partition = 1;
cells = face->getCells();
for(std::list< std::pair<int, Cell*> >::iterator it = cells.begin(); it != cells.end(); it++){
Cell* cell = (*it).second;
- fprintf(fp, "%d %d %d %d %d %d %d %d %d\n", cell->getNum(), 2, 3, 0, 0, physical, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum());
+ if(cell->getNumVertices() == 3) fprintf(fp, "%d %d %d %d %d %d %d %d %d\n", cell->getNum(), 2, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum());
+ else if (cell->getNumVertices() == 4) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), 3, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum());
}
}
for(citer cit = firstCell(3); cit != lastCell(3); cit++) {
Cell* volume = *cit;
- if(volume->inSubdomain()) physical = 3;
- else if(volume->onDomainBoundary()) physical = 2;
- else physical = 1;
+ if(volume->inSubdomain()) partition = 3;
+ else if(volume->onDomainBoundary()) partition = 2;
+ else partition = 1;
cells = volume->getCells();
for(std::list< std::pair<int, Cell*> >::iterator it = cells.begin(); it != cells.end(); it++){
Cell* cell = (*it).second;
- fprintf(fp, "%d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), 4, 3, 0, 0, physical, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum());
+ if(cell->getNumVertices() == 4) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), 4, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum());
+ if(cell->getNumVertices() == 8) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), 12, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(), cell->getVertex(4)->getNum(), cell->getVertex(5)->getNum(), cell->getVertex(6)->getNum(), cell->getVertex(7)->getNum() );
+ if(cell->getNumVertices() == 6) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), 13, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(),cell->getVertex(4)->getNum(), cell->getVertex(5)->getNum());
+ if(cell->getNumVertices() == 5) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), 14, 3, 0, 0, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(), cell->getVertex(4)->getNum());
}
}
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index ad55b013c9..3f6bdd0983 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -18,6 +18,10 @@
#include "MPoint.h"
#include "MLine.h"
#include "MTriangle.h"
+#include "MQuadrangle.h"
+#include "MHexahedron.h"
+#include "MPrism.h"
+#include "MPyramid.h"
#include "MTetrahedron.h"
#include "GModel.h"
#include "GEntity.h"
@@ -58,9 +62,10 @@ class Cell
public:
Cell(){
-
_bdSize = 0;
_cbdSize = 0;
+ _combined = false;
+ _index = 0;
}
virtual ~Cell(){}
@@ -80,7 +85,7 @@ class Cell
// returns 1 or -1 if lower dimensional cell tau is on the boundary of this cell
// otherwise returns 0
// implementation will vary depending on cell type
- virtual int kappa(Cell* tau) const = 0;
+ //virtual int kappa(Cell* tau) const = 0;
// true if this cell has given vertex
virtual bool hasVertex(int vertex) const = 0;
@@ -233,7 +238,7 @@ class Cell
virtual int getNumCells() {return 1;}
bool operator==(const Cell& c2){
-
+
if(this->getNumVertices() != c2.getNumVertices()){
return false;
}
@@ -243,12 +248,14 @@ class Cell
}
}
return true;
-
+
+ //return (this->getTag() == c2.getTag());
}
bool operator<(const Cell& c2) const {
+
if(this->getNumVertices() != c2.getNumVertices()){
return (this->getNumVertices() < c2.getNumVertices());
}
@@ -257,8 +264,11 @@ class Cell
else if (this->getSortedVertex(i) > c2.getSortedVertex(i)) return false;
}
return false;
+
+ //return (this->getTag() < c2.getTag());
}
+ virtual int kappa(Cell* tau) const;
};
@@ -268,15 +278,12 @@ class Simplex : public Cell
{
public:
- Simplex() : Cell() {
- _combined = false;
- _index = 0;
- }
+ Simplex() : Cell() {}
~Simplex(){}
// kappa for simplices
// checks whether lower dimensional simplex tau is on the boundary of this cell
- virtual int kappa(Cell* tau) const;
+ //virtual int kappa(Cell* tau) const;
};
@@ -348,7 +355,7 @@ class OneSimplex : public Simplex, public MLine
v[0] = _v[num];
}
- //int kappa(Cell* tau) const;
+ int kappa(Cell* tau) const;
void printCell() const { printf("Cell dimension: %d, Vertices: %d %d, in subdomain: %d \n", getDim(), _v[0]->getNum(), _v[1]->getNum(), _inSubdomain); }
};
@@ -397,7 +404,53 @@ class TwoSimplex : public Simplex, public MTriangle
void printCell() const { printf("Cell dimension: %d, Vertices: %d %d %d, in subdomain: %d\n", getDim(), _v[0]->getNum(), _v[1]->getNum(), _v[2]->getNum(), _inSubdomain); }
};
-// Three simplex cell type.
+// Quadrangle cell type.
+class CQuadrangle : public Cell, public MQuadrangle
+{
+ private:
+
+ int _vs[4];
+
+ public:
+
+ CQuadrangle(std::vector<MVertex*> v, int tag=0, bool subdomain=false, bool boundary=false, int num=0, int part=0)
+ : MQuadrangle(v, num, part){
+ _tag = tag;
+ _inSubdomain = subdomain;
+ _onDomainBoundary = boundary;
+ _vs[0] = v.at(0)->getNum();
+ _vs[1] = v.at(1)->getNum();
+ _vs[2] = v.at(2)->getNum();
+ _vs[3] = v.at(3)->getNum();
+ std::sort(_vs, _vs+3);
+ }
+
+ ~CQuadrangle(){}
+
+ int getDim() const { return 2; }
+ int getNum() { return MQuadrangle::getNum(); }
+ int getNumVertices() const { return 4; }
+ int getNumFacets() const { return 4; }
+ MVertex* getVertex(int vertex) const {return _v[vertex]; }
+ int getSortedVertex(int vertex) const {return _vs[vertex]; }
+ bool hasVertex(int vertex) const {return
+ (_v[0]->getNum() == vertex || _v[1]->getNum() == vertex || _v[2]->getNum() == vertex || _v[3]->getNum() == vertex); }
+
+ std::vector<MVertex*> getVertexVector() const {
+ return std::vector<MVertex*>(_v, _v + sizeof(_v)/sizeof(MVertex*) ); }
+ std::vector<int> getSortedVertexVector() const {
+ return std::vector<int>(_vs, _vs + sizeof(_vs)/sizeof(int) );
+ }
+
+ void getFacetVertices(const int num, std::vector<MVertex*> &v) const {
+ MQuadrangle::getEdgeVertices(num, v);
+ }
+
+ void printCell() const { printf("Cell dimension: %d, Vertices: %d %d %d, in subdomain: %d\n", getDim(), _v[0]->getNum(), _v[1]->getNum(), _v[2]->getNum(), _v[3]->getNum(), _inSubdomain); }
+};
+
+
+// ThreeSimplex cell type.
class ThreeSimplex : public Simplex, public MTetrahedron
{
private:
@@ -438,6 +491,142 @@ class ThreeSimplex : public Simplex, public MTetrahedron
virtual void printCell() const { printf("Cell dimension: %d, Vertices: %d %d %d %d, in subdomain: %d \n", getDim(), _v[0]->getNum(), _v[1]->getNum(), _v[2]->getNum(), _v[3]->getNum(), _inSubdomain); }
};
+// Hexahedron cell type.
+class CHexahedron : public Cell, public MHexahedron
+{
+ private:
+ int _vs[8];
+
+ public:
+ CHexahedron(std::vector<MVertex*> v, int tag=0, bool subdomain=false, bool boundary=false, int num=0, int part=0)
+ : MHexahedron(v, num, part){
+ _tag = tag;
+ _inSubdomain = subdomain;
+ _onDomainBoundary = boundary;
+ _vs[0] = v.at(0)->getNum();
+ _vs[1] = v.at(1)->getNum();
+ _vs[2] = v.at(2)->getNum();
+ _vs[3] = v.at(3)->getNum();
+ _vs[4] = v.at(4)->getNum();
+ _vs[5] = v.at(5)->getNum();
+ _vs[6] = v.at(6)->getNum();
+ _vs[7] = v.at(7)->getNum();
+ std::sort(_vs, _vs+8);
+ }
+
+ ~CHexahedron(){}
+
+ int getDim() const { return 3; }
+ int getNum() { return MHexahedron::getNum(); }
+ int getNumVertices() const { return 8; }
+ int getNumFacets() const { return 6; }
+ MVertex* getVertex(int vertex) const {return _v[vertex]; }
+ int getSortedVertex(int vertex) const {return _vs[vertex]; }
+ bool hasVertex(int vertex) const {return
+ (_v[0]->getNum() == vertex || _v[1]->getNum() == vertex || _v[2]->getNum() == vertex || _v[3]->getNum() == vertex ||
+ _v[4]->getNum() == vertex || _v[5]->getNum() == vertex || _v[6]->getNum() == vertex || _v[7]->getNum() == vertex ); }
+ std::vector<MVertex*> getVertexVector() const {
+ return std::vector<MVertex*>(_v, _v + sizeof(_v)/sizeof(MVertex*) ); }
+ std::vector<int> getSortedVertexVector() const {
+ return std::vector<int>(_vs, _vs + sizeof(_vs)/sizeof(int) ); }
+
+ void getFacetVertices(const int num, std::vector<MVertex*> &v) const {
+ MHexahedron::getFaceVertices(num, v);
+ }
+
+ virtual void printCell() const { printf("Cell dimension: %d, Vertices: %d %d %d %d %d %d %d %d, in subdomain: %d \n", getDim(), _v[0]->getNum(), _v[1]->getNum(), _v[2]->getNum(), _v[3]->getNum(), _v[4]->getNum(), _v[5]->getNum(), _v[6]->getNum(), _v[7]->getNum(), _inSubdomain); }
+};
+
+
+// Prism cell type.
+class CPrism : public Cell, public MPrism
+{
+ private:
+ int _vs[6];
+
+ public:
+ CPrism(std::vector<MVertex*> v, int tag=0, bool subdomain=false, bool boundary=false, int num=0, int part=0)
+ : MPrism(v, num, part){
+ _tag = tag;
+ _inSubdomain = subdomain;
+ _onDomainBoundary = boundary;
+ _vs[0] = v.at(0)->getNum();
+ _vs[1] = v.at(1)->getNum();
+ _vs[2] = v.at(2)->getNum();
+ _vs[3] = v.at(3)->getNum();
+ _vs[4] = v.at(4)->getNum();
+ _vs[5] = v.at(5)->getNum();
+ std::sort(_vs, _vs+6);
+ }
+
+ ~CPrism(){}
+
+ int getDim() const { return 3; }
+ int getNum() { return MPrism::getNum(); }
+ int getNumVertices() const { return 6; }
+ int getNumFacets() const { return 5; }
+ MVertex* getVertex(int vertex) const {return _v[vertex]; }
+ int getSortedVertex(int vertex) const {return _vs[vertex]; }
+ bool hasVertex(int vertex) const {return
+ (_v[0]->getNum() == vertex || _v[1]->getNum() == vertex || _v[2]->getNum() == vertex || _v[3]->getNum() == vertex ||
+ _v[4]->getNum() == vertex || _v[5]->getNum() == vertex); }
+ std::vector<MVertex*> getVertexVector() const {
+ return std::vector<MVertex*>(_v, _v + sizeof(_v)/sizeof(MVertex*) ); }
+ std::vector<int> getSortedVertexVector() const {
+ return std::vector<int>(_vs, _vs + sizeof(_vs)/sizeof(int) ); }
+
+ void getFacetVertices(const int num, std::vector<MVertex*> &v) const {
+ MPrism::getFaceVertices(num, v);
+ }
+
+ virtual void printCell() const { printf("Cell dimension: %d, Vertices: %d %d %d %d %d %d, in subdomain: %d \n", getDim(), _v[0]->getNum(), _v[1]->getNum(), _v[2]->getNum(), _v[3]->getNum(), _v[4]->getNum(), _v[5]->getNum(), _inSubdomain); }
+};
+
+
+// Pyramid cell type.
+class CPyramid : public Cell, public MPyramid
+{
+ private:
+ int _vs[5];
+
+ public:
+ CPyramid(std::vector<MVertex*> v, int tag=0, bool subdomain=false, bool boundary=false, int num=0, int part=0)
+ : MPyramid(v, num, part){
+ _tag = tag;
+ _inSubdomain = subdomain;
+ _onDomainBoundary = boundary;
+ _vs[0] = v.at(0)->getNum();
+ _vs[1] = v.at(1)->getNum();
+ _vs[2] = v.at(2)->getNum();
+ _vs[3] = v.at(3)->getNum();
+ _vs[4] = v.at(4)->getNum();
+ std::sort(_vs, _vs+5);
+ }
+
+ ~CPyramid(){}
+
+ int getDim() const { return 3; }
+ int getNum() { return MPyramid::getNum(); }
+ int getNumVertices() const { return 5; }
+ int getNumFacets() const { return 5; }
+ MVertex* getVertex(int vertex) const {return _v[vertex]; }
+ int getSortedVertex(int vertex) const {return _vs[vertex]; }
+ bool hasVertex(int vertex) const {return
+ (_v[0]->getNum() == vertex || _v[1]->getNum() == vertex || _v[2]->getNum() == vertex || _v[3]->getNum() == vertex ||
+ _v[4]->getNum() == vertex); }
+ std::vector<MVertex*> getVertexVector() const {
+ return std::vector<MVertex*>(_v, _v + sizeof(_v)/sizeof(MVertex*) ); }
+ std::vector<int> getSortedVertexVector() const {
+ return std::vector<int>(_vs, _vs + sizeof(_vs)/sizeof(int) ); }
+
+ void getFacetVertices(const int num, std::vector<MVertex*> &v) const {
+ MPyramid::getFaceVertices(num, v);
+ }
+
+ virtual void printCell() const { printf("Cell dimension: %d, Vertices: %d %d %d %d %d, in subdomain: %d \n", getDim(), _v[0]->getNum(), _v[1]->getNum(), _v[2]->getNum(), _v[3]->getNum(), _v[4]->getNum(), _inSubdomain); }
+};
+
+
// Ordering for the cells.
class Less_Cell{
public:
@@ -459,6 +648,10 @@ class Less_Cell{
}
return false;
+
+
+ //return (c1->getTag() < c2->getTag());
+
}
};
@@ -466,6 +659,8 @@ class Less_Cell{
class Equal_Cell{
public:
bool operator ()(const Cell* c1, const Cell* c2){
+
+
if(c1->getNumVertices() != c2->getNumVertices()){
return false;
}
@@ -475,6 +670,9 @@ class Equal_Cell{
}
}
return true;
+
+ //return (c1->getTag() == c2->getTag());
+
}
};
@@ -595,7 +793,13 @@ class CombinedCell : public Cell{
}
- ~CombinedCell(){}
+ ~CombinedCell(){
+ std::list< std::pair<int, Cell*> > cells = getCells();
+ for(std::list< std::pair<int, Cell*> >::iterator it = cells.begin(); it != cells.end(); it++){
+ Cell* cell2 = (*it).second;
+ delete cell2;
+ }
+ }
int getNum() { return _num; }
int getNumVertices() const { return _v.size(); }
@@ -650,7 +854,9 @@ class CellComplex
std::set<Cell*, Less_Cell> _cells[4];
// storage for cell pointers to delete them
- std::set<Cell*, Less_Cell> _cells2[4];
+ //std::set<Cell*, Less_Cell> _cells2[4];
+ std::list<Cell*> _trash;
+
//std::set<Cell*, Less_Cell> _originalCells[4];
@@ -675,7 +881,7 @@ class CellComplex
void insert_cells(bool subdomain, bool boundary);
// remove a cell from this cell complex
- void removeCell(Cell* cell);
+ void removeCell(Cell* cell, bool other=true);
void insertCell(Cell* cell);
public:
@@ -722,23 +928,33 @@ class CellComplex
CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> subdomain );
CellComplex(){}
- ~CellComplex(){
+ ~CellComplex(){
for(int i = 0; i < 4; i++){
for(citer cit = _cells[i].begin(); cit != _cells[i].end(); cit++){
Cell* cell = *cit;
- if(cell->isCombined()) delete cell;
+ delete cell;
+ //if(cell->isCombined()) delete cell;
}
_cells[i].clear();
+ /*
for(citer cit = _cells2[i].begin(); cit != _cells2[i].end(); cit++){
Cell* cell = *cit;
delete cell;
}
_cells2[i].clear();
+ */
}
-
+ emptyTrash();
}
+ void emptyTrash(){
+ for(std::list<Cell*>::iterator cit = _trash.begin(); cit != _trash.end(); cit++){
+ Cell* cell = *cit;
+ delete cell;
+ }
+ _trash.clear();
+ }
// get the number of certain dimensional cells
int getSize(int dim){ return _cells[dim].size(); }
@@ -768,7 +984,7 @@ class CellComplex
bool inSameDomain(Cell* c1, Cell* c2) const { return
( (!c1->inSubdomain() && !c2->inSubdomain()) || (c1->inSubdomain() && c2->inSubdomain()) ); }
-
+
// useful functions for (co)reduction of cell complex
int reduceComplex(int omit = 0);
int coreduceComplex(int omit = 0);
@@ -799,6 +1015,11 @@ class CellComplex
// also check whether all (co)boundary cells of a cell belong to this cell complex
bool checkCoherence();
+ int eulerCharacteristic(){
+ return getSize(0) - getSize(1) + getSize(2) - getSize(3);
+ }
+ void printEuler(){ printf("Euler characteristic: %d. \n", eulerCharacteristic()); }
+
// change roles of boundaries and coboundaries of the cells in this cell complex
// equivalent to transposing boundary operator matrices
void makeDualComplex(){
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index 650ac8a1dc..c3ff48d103 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -77,7 +77,11 @@ void Homology::findGenerators(std::string fileName){
Msg::Info("Reducing Cell Complex...");
double t1 = Cpu();
+ //_cellComplex->printEuler();
int omitted = _cellComplex->reduceComplex(_omit);
+ //_cellComplex->printEuler();
+
+ _cellComplex->emptyTrash();
if(getCombine()){
_cellComplex->combine(3);
@@ -87,6 +91,9 @@ void Homology::findGenerators(std::string fileName){
_cellComplex->combine(1);
}
_cellComplex->checkCoherence();
+ //_cellComplex->printEuler();
+
+ _cellComplex->emptyTrash();
double t2 = Cpu();
Msg::Info("Cell Complex reduction complete (%g s).", t2 - t1);
@@ -130,7 +137,7 @@ void Homology::findGenerators(std::string fileName){
Msg::Info("H1 = %d", HRank[1]);
Msg::Info("H2 = %d", HRank[2]);
Msg::Info("H3 = %d", HRank[3]);
- if(omitted != 0) Msg::Info("Computation of generators in %d highest dimensions was omitted.", _omit);
+ if(omitted != 0) Msg::Info("The computation of generators in %d highest dimensions was omitted.", _omit);
Msg::Info("Wrote results to %s.", fileName.c_str());
@@ -144,13 +151,15 @@ void Homology::findGenerators(std::string fileName){
return;
}
-void Homology::findThickCuts(std::string fileName){
+void Homology::findDualGenerators(std::string fileName){
//for(int i = 0; i < 4; i++) { printf("Dim %d: \n", i); _cellComplex->printComplex(i); }
Msg::Info("Reducing Cell Complex...");
double t1 = Cpu();
int omitted = _cellComplex->coreduceComplex(_omit);
+ _cellComplex->emptyTrash();
+
/*
_cellComplex->makeDualComplex();
int omitted = _cellComplex->reduceComplex(_omit);
@@ -168,6 +177,8 @@ void Homology::findThickCuts(std::string fileName){
_cellComplex->cocombine(2);
}
+ _cellComplex->emptyTrash();
+
_cellComplex->checkCoherence();
double t2 = Cpu();
Msg::Info("Cell Complex reduction complete (%g s).", t2 - t1);
@@ -200,12 +211,12 @@ void Homology::findThickCuts(std::string fileName){
convert(i, generator);
convert(dim-j, dimension);
- std::string name = dimension + "D Thick cut " + generator;
+ std::string name = dimension + "D Dual generator " + generator;
Chain* chain = new Chain(_cellComplex->getCells(j), chains->getCoeffVector(j,i), _cellComplex, name, chains->getTorsion(j,i));
chain->writeChainMSH(fileName);
if(chain->getSize() != 0){
HRank[dim-j] = HRank[dim-j] + 1;
- if(chain->getTorsion() != 1) Msg::Warning("%dD Thick cut %d has torsion coefficient %d!", j, i, chain->getTorsion());
+ if(chain->getTorsion() != 1) Msg::Warning("%dD Dual generator %d has torsion coefficient %d!", j, i, chain->getTorsion());
}
delete chain;
@@ -217,7 +228,7 @@ void Homology::findThickCuts(std::string fileName){
Msg::Info("H1 = %d", HRank[1]);
Msg::Info("H2 = %d", HRank[2]);
Msg::Info("H3 = %d", HRank[3]);
- if(omitted != 0) Msg::Info("Computation of %d highest dimension thick cuts was omitted.", _omit);
+ if(omitted != 0) Msg::Info("The computation of %d highest dimension dual generators was omitted.", _omit);
Msg::Info("Wrote results to %s.", fileName.c_str());
@@ -230,4 +241,21 @@ void Homology::findThickCuts(std::string fileName){
return;
}
+
+void Homology::computeBettiNumbers(){
+
+ Msg::Info("Running coreduction...");
+ double t1 = Cpu();
+ _cellComplex->computeBettiNumbers();
+ double t2 = Cpu();
+ Msg::Info("Betti number computation complete (%g s).", t2- t1);
+
+ Msg::Info("b0 = %d \n", _cellComplex->getBettiNumber(0));
+ Msg::Info("b1 = %d \n", _cellComplex->getBettiNumber(1));
+ Msg::Info("b2 = %d \n", _cellComplex->getBettiNumber(2));
+ Msg::Info("b3 = %d \n", _cellComplex->getBettiNumber(3));
+
+ return;
+}
+
#endif
diff --git a/Geo/Homology.h b/Geo/Homology.h
index c5dad2dbe2..8be0e76cb5 100644
--- a/Geo/Homology.h
+++ b/Geo/Homology.h
@@ -40,8 +40,9 @@ class Homology
bool setCombine(bool combine) { _combine = combine; return _combine; }
void findGenerators(std::string fileName);
- void findThickCuts(std::string fileName);
-
+ void findDualGenerators(std::string fileName);
+ void computeBettiNumbers();
+
void swapSubdomain() { _cellComplex->swapSubdomain(); }
int getOmit() {return _omit; }
diff --git a/Plugin/HomologyComputation.cpp b/Plugin/HomologyComputation.cpp
index 0f0fa194c1..337b2a79c8 100644
--- a/Plugin/HomologyComputation.cpp
+++ b/Plugin/HomologyComputation.cpp
@@ -21,8 +21,9 @@ StringXNumber HomologyComputationOptions_Number[] = {
{GMSH_FULLRC, "PhysicalGroupForSubdomain1", NULL, 0.},
{GMSH_FULLRC, "PhysicalGroupForSubdomain2", NULL, 0.},
{GMSH_FULLRC, "ComputeGenerators", NULL, 1.},
- {GMSH_FULLRC, "ComputeThickCuts", NULL, 0.},
- {GMSH_FULLRC, "OmitDimensions", NULL, 1.},
+ {GMSH_FULLRC, "ComputeDualGenerators", NULL, 0.},
+ {GMSH_FULLRC, "ComputeBettiNumbers", NULL, 0.},
+ //{GMSH_FULLRC, "OmitDimensions", NULL, 1.},
//{GMSH_FULLRC, "CombineCells", NULL, 1.},
};
@@ -84,7 +85,8 @@ PView *GMSH_HomologyComputationPlugin::execute(PView *v)
int gens = (int)HomologyComputationOptions_Number[4].def;
int cuts = (int)HomologyComputationOptions_Number[5].def;
- int omit = (int)HomologyComputationOptions_Number[6].def;
+ int betti = (int)HomologyComputationOptions_Number[6].def;
+ //int omit = (int)HomologyComputationOptions_Number[6].def;
//int combine = (int)HomologyComputationOptions_Number[7].def;
@@ -92,19 +94,22 @@ PView *GMSH_HomologyComputationPlugin::execute(PView *v)
Homology* homology = new Homology(m, domain, subdomain);
//if(combine == 0) homology->setCombine(false);
- homology->setOmit(omit);
+ homology->setOmit(1);
//if(swap == 1) homology->swapSubdomain();
- if(gens == 1 && cuts != 1) {
+ if(gens == 1 && cuts != 1 && betti != 1) {
homology->findGenerators(fileName);
GmshMergeFile(fileName);
}
- else if(cuts == 1 && gens != 1) {
- homology->findThickCuts(fileName);
+ else if(cuts == 1 && gens != 1 && betti != 1) {
+ homology->findDualGenerators(fileName);
GmshMergeFile(fileName);
}
- else Msg::Error("Choose either generators or thick cuts to compute.");
+ else if(cuts != 1 && gens != 1 && betti == 1) {
+ homology->computeBettiNumbers();
+ }
+ else Msg::Error("Choose either generators, dual generators or Betti numbers to compute.");
delete homology;
--
GitLab