diff --git a/Geo/Cell.cpp b/Geo/Cell.cpp
index 962a53b3ad9c2239fea66986b1416cfd863c1d34..fdfb1a6fd44c1d9e0c0f1c42b16647a5e4105ef8 100755
--- a/Geo/Cell.cpp
+++ b/Geo/Cell.cpp
@@ -243,7 +243,6 @@ CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell() {
}
_index = c1->getIndex();
- //_tag = c1->getTag();
_dim = c1->getDim();
_num = c1->getNum();
_inSubdomain = c1->inSubdomain();
diff --git a/Geo/Cell.h b/Geo/Cell.h
index 3fb2439e8300f8ffb4c4b76ab0f65a5634216d6a..dd5b41bb3b2f5595df8f44e71062c32741c9f34e 100644
--- a/Geo/Cell.h
+++ b/Geo/Cell.h
@@ -169,7 +169,7 @@ class Cell
virtual void printOrgCbd();
virtual bool inSubdomain() const { return _inSubdomain; }
- //virtual void setInSubdomain(bool subdomain) { _inSubdomain = subdomain; }
+ virtual void setInSubdomain(bool subdomain) { _inSubdomain = subdomain; }
virtual bool onDomainBoundary() const { return _onDomainBoundary; }
virtual void setOnDomainBoundary(bool domainboundary) { _onDomainBoundary = domainboundary; }
diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index cb2b779a5a8bb1914192f955c3044561acfe870f..cd88192b041548c497c53806307439547d7d9018 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -18,9 +18,20 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su
_dim = 0;
_simplicial = true;
- // find boundary entities
- find_boundary(domain, subdomain);
+ _multidim = false;
+ int dim = 0;
+ for(int i = 0; i < domain.size(); i++){
+ GEntity* entity = domain.at(i);
+ if(i == 0) dim = entity->dim();
+ if(dim != entity->dim()){
+ _multidim = true;
+ Msg::Warning("Domain is not a manifold.");
+ }
+ }
+ // find boundary entities
+ if(!_multidim) find_boundary(domain, subdomain);
+
// insert cells into cell complex
// subdomain need to be inserted first!
insert_cells(true, true);
@@ -34,7 +45,6 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su
if(getSize(i) > _dim) _dim = i;
}
-
}
void CellComplex::find_boundary(std::vector<GEntity*>& domain, std::vector<GEntity*>& subdomain){
@@ -97,8 +107,6 @@ void CellComplex::find_boundary(std::vector<GEntity*>& domain, std::vector<GEnti
}
void CellComplex::insert_cells(bool subdomain, bool boundary){
-
- // works only for simplicial meshes at the moment!
std::vector<GEntity*> domain;
@@ -134,17 +142,17 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){
else if(type == MSH_QUA_4 || type == MSH_QUA_8 || type == MSH_QUA_9){
cell = new Cell(domain.at(j)->getMeshElement(i), subdomain, boundary);
_simplicial = false;
- }/* FIXME: reduction doesn't work for these (but should).
+ }
else if(type == MSH_HEX_8 || type == MSH_HEX_27 || type == MSH_HEX_20){
- cell = new CHexahedron(vertices, tag, subdomain, boundary);
+ cell = new Cell(domain.at(j)->getMeshElement(i), subdomain, boundary);
_simplicial = false;
- }
+ }/* FIXME: no getFaceInfo methods for these MElements
else if(type == MSH_PRI_6 || type == MSH_PRI_18 || type == MSH_PRI_15){
- cell = new CPrism(vertices, tag, subdomain, boundary);
+ cell = new Cell(domain.at(j)->getMeshElement(i), subdomain, boundary);
_simplicial = false;
}
else if(type == MSH_PYR_5 || type == MSH_PYR_14 || type == MSH_PYR_13){
- cell = new CPyramid(vertices, tag, subdomain, boundary);
+ cell = new Cell(domain.at(j)->getMeshElement(i), subdomain, boundary);
_simplicial = false;
}*/
else Msg::Error("Error: mesh element %d not implemented yet! \n", type);
@@ -163,23 +171,31 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){
for(int i = 0; i < cell->getNumFacets(); i++){
cell->getFacetVertices(i, vertices);
int type = cell->getTypeForMSH();
-
int newtype = 0;
//FIXME: add missing boundary cell type relations
+ //FIXME: high order meshes don't work
if(dim == 3){
if(type == MSH_TET_4) newtype = MSH_TRI_3;
- else if(type == MSH_TET_10) newtype = MSH_TRI_6;
+ /*else if(type == MSH_TET_10) newtype = MSH_TRI_6;
+ else if(type == MSH_TET_20) newtype = MSH_TRI_9;*/
+ else if(type == MSH_HEX_8) newtype = MSH_QUA_4;
+ /*else if(type == MSH_HEX_20) newtype = MSH_QUA_8;
+ else if(type == MSH_HEX_27) newtype = MSH_QUA_9;*/
}
else if(dim == 2){
if(type == MSH_TRI_3 || type == MSH_QUA_4) newtype = MSH_LIN_2;
- else if(type == MSH_TRI_6 || type == MSH_QUA_8) newtype = MSH_LIN_3;
+ /*else if(type == MSH_TRI_6 || type == MSH_QUA_8) newtype = MSH_LIN_3;
+ else if(type == MSH_TRI_9) newtype = MSH_LIN_4;
+ else if(type == MSH_QUA_9) newtype = MSH_LIN_3;*/
}
else if(dim == 1){
- if(type == MSH_LIN_2 || type == MSH_LIN_3 || type == MSH_LIN_4 ||
- type == MSH_LIN_5 || type == MSH_LIN_6) newtype = MSH_PNT;
+ if(type == MSH_LIN_2) newtype = MSH_PNT;
+ /*else if(type == MSH_LIN_3 || type == MSH_LIN_4 ||
+ type == MSH_LIN_5 || type == MSH_LIN_6) newtype = MSH_PNT;*/
}
if(newtype == 0) Msg::Error("Error: mesh element %d not implemented yet! \n", type);
+ //printf("dim: %d type: %d, newtype: %d, vertices: %d \n", dim, type, newtype, vertices.size());
MElement* element = factory.create(newtype, vertices, 0, cell->getPartition());
Cell* newCell = new Cell(element, subdomain, boundary);
newCell->setImmune(cell->getImmune());
@@ -206,7 +222,6 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){
}
}
}
-
}
CellComplex::~CellComplex(){
@@ -893,6 +908,18 @@ bool CellComplex::hasCell(Cell* cell, bool org){
}
}
+bool CellComplex::swapSubdomain(){
+ if(_multidim) return false;
+ 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 true;
+}
+
void CellComplex::makeDualComplex(){
std::set<Cell*, Less_Cell> temp = _cells[0];
_cells[0] = _cells[3];
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index 49c9ffc07270df86eb1c188ed93d9d21e6157c3d..017925888365673ca31c0b2c4782cfd4c7201318 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -59,6 +59,9 @@ class CellComplex
// is the cell complex simplicial
bool _simplicial;
+ // does the domain contain entities of different dimensions
+ bool _multidim;
+
// enqueue cells in queue if they are not there already
void enqueueCells(std::list<Cell*>& cells,
std::queue<Cell*>& Q, std::set<Cell*, Less_Cell>& Qset);
@@ -129,10 +132,6 @@ class CellComplex
// print the vertices of cells of certain dimension
void printComplex(int dim);
- // write this cell complex in 2.0 MSH ASCII format
- // for debugging only
- //int writeComplexMSH(const std::string &name);
-
// Cell combining for reduction and coreduction
int combine(int dim);
int cocombine(int dim);
@@ -145,6 +144,9 @@ class CellComplex
// also check whether all (co)boundary cells of a cell belong to this cell complex
bool checkCoherence();
+ // make cells on the domain boundary that are not in the subdomain subdomain cells and vice versa
+ bool swapSubdomain();
+
int eulerCharacteristic(){ return getSize(0) - getSize(1) + getSize(2) - getSize(3);}
void printEuler(){ printf("Euler characteristic: %d. \n", eulerCharacteristic()); }
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index 231bee797bbc7685b6601040e13f9ec92fc195cf..cbfd41bb1897c5d54cb0f17563f4473222007dcd 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -602,7 +602,6 @@ int Chain::writeChainMSH(const std::string &name){
void Chain::createPView(){
std::vector<MElement*> elements;
- MElementFactory factory;
std::map<int, std::vector<double> > data;
for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index fe8850d0dbfc146d6727a6737a6fe8a59d8f2411..1a5344f3ea95c33466c064feb59943c8873bf325 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -111,7 +111,6 @@ void Homology::findGenerators(std::string fileName){
_cellComplex->getSize(3), _cellComplex->getSize(2), _cellComplex->getSize(1), _cellComplex->getSize(0));
Msg::StatusBar(2, false, "%d V, %d F, %d E, %d N.",
_cellComplex->getSize(3), _cellComplex->getSize(2), _cellComplex->getSize(1), _cellComplex->getSize(0));
- //for(int i = 0; i < 4; i++) { printf("Dim %d: \n", i); _cellComplex->printComplex(i); }
Msg::Info("Computing homology spaces...");
Msg::StatusBar(1, false, "Computing...");
@@ -190,7 +189,6 @@ void Homology::findDualGenerators(std::string fileName){
Msg::StatusBar(1, false, "Reducing...");
double t1 = Cpu();
int omitted = _cellComplex->coreduceComplex();
- //_cellComplex->emptyTrash();
/*
_cellComplex->makeDualComplex();
@@ -203,14 +201,10 @@ void Homology::findDualGenerators(std::string fileName){
_cellComplex->makeDualComplex();
*/
-
_cellComplex->cocombine(0);
_cellComplex->cocombine(1);
_cellComplex->cocombine(2);
-
- //_cellComplex->emptyTrash();
-
_cellComplex->checkCoherence();
double t2 = Cpu();
Msg::Info("Cell Complex reduction complete (%g s).", t2 - t1);
@@ -220,8 +214,6 @@ void Homology::findDualGenerators(std::string fileName){
_cellComplex->getSize(3), _cellComplex->getSize(2), _cellComplex->getSize(1), _cellComplex->getSize(0));
- //_cellComplex->writeComplexMSH(fileName);
-
Msg::Info("Computing homology spaces...");
Msg::StatusBar(1, false, "Computing...");
t1 = Cpu();
@@ -360,16 +352,8 @@ void Homology::createPViews(){
bool Homology::writeGeneratorsMSH(std::string fileName, bool binary){
if(!_model->writeMSH(fileName, 2.0, binary)) return false;
- /*
- for(int i = 0; i < 4; i++){
- for(int j = 0; j < _generators[i].size(); j++){
- Chain* chain = _generators[i].at(j);
- if(!chain->writeChainMSH(fileName)) return false;
- }
- }*/
Msg::Info("Wrote homology computation results to %s.", fileName.c_str());
- Msg::Debug("Wrote homology computation results to %s. \n", fileName.c_str());
-
+ Msg::Debug("Wrote homology computation results to %s. \n", fileName.c_str());
return true;
}
diff --git a/Geo/Homology.h b/Geo/Homology.h
index fbc47be58c61358189172e4a422f9c0e423f5a59..0443ba40fdcf4f486572dc36440350953100b8b0 100644
--- a/Geo/Homology.h
+++ b/Geo/Homology.h
@@ -47,9 +47,8 @@ class Homology
void findGenerators(std::string fileName);
void findDualGenerators(std::string fileName);
void computeBettiNumbers();
-
-
- //void swapSubdomain() { _cellComplex->swapSubdomain(); }
+
+ bool swapSubdomain() { return _cellComplex->swapSubdomain(); }
// Restore the cell complex to its original state before cell reductions
void restoreHomology();
diff --git a/Parser/Gmsh.y b/Parser/Gmsh.y
index 813628da5e1c52a00239fb044f13eabe0c2e2d52..dd1cf5962c28f125a76e5e18ee03e20795a85955 100644
--- a/Parser/Gmsh.y
+++ b/Parser/Gmsh.y
@@ -3333,7 +3333,7 @@ Homology :
#if defined(HAVE_KBIPACK)
Homology* homology = new Homology(GModel::current(), domain, subdomain);
- homology->findGenerators(fileName);
+ homology->findGenerators(fileName);
delete homology;
#else
yymsg(0, "Gmsh needs to be configured with option Kbipack to use homology computation.");