diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 46eb4b36c28e16b615a6940133aeaeb5dc6f627d..2eb31a86ac99368be821596280cc6da95fd2b598 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 ad55b013c9efe8050472a722058712f4a88dbff1..3f6bdd09830136854a19ee3c742f534840b018eb 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 650ac8a1dcbb8bebf5955dd932f4049e00f9c38d..c3ff48d10313e2972d636b8e01a64c53bee2eee4 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 c5dad2dbe24d8c203f4b8529e2981c07be5ba32b..8be0e76cb5ea9ff5803a371c526ad00e51b85ed7 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 0f0fa194c1c144c2b7885ffb3bde9aa538ec0239..337b2a79c83a3782e2ca36c4712bcf66988d0329 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;