diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 2f2fd3afde076c0acdfa5c3f746c20847fe7bede..226a7b1273f2db7a8372128b58107f0ef186c558 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -10,8 +10,9 @@
CellComplex::CellComplex(GModel* model,
std::vector<MElement*>& domainElements,
- std::vector<MElement*>& subdomainElements) :
- _model(model), _dim(0), _simplicial(true), _saveorig(true)
+ std::vector<MElement*>& subdomainElements,
+ bool saveOriginalComplex) :
+ _model(model), _dim(0), _simplicial(true), _saveorig(saveOriginalComplex)
{
_insertCells(subdomainElements, 1);
@@ -20,7 +21,7 @@ CellComplex::CellComplex(GModel* model,
int num = 0;
for(int dim = 0; dim < 4; dim++){
if(getSize(dim) != 0) _dim = dim;
- _ocells[dim] = _cells[dim];
+ if(_saveorig) _ocells[dim] = _cells[dim];
for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
Cell* cell = *cit;
cell->setNum(++num);
@@ -71,9 +72,17 @@ bool CellComplex::_insertCells(std::vector<MElement*>& elements,
CellComplex::~CellComplex()
{
for(int i = 0; i < 4; i++){
- for(citer cit = _ocells[i].begin(); cit != _ocells[i].end(); cit++){
- Cell* cell = *cit;
- delete cell;
+ if(_saveorig) {
+ for(citer cit = _ocells[i].begin(); cit != _ocells[i].end(); cit++){
+ Cell* cell = *cit;
+ delete cell;
+ }
+ }
+ else {
+ for(citer cit = _cells[i].begin(); cit != _cells[i].end(); cit++){
+ Cell* cell = *cit;
+ delete cell;
+ }
}
}
for(unsigned int i = 0; i < _newcells.size(); i++) delete _newcells.at(i);
@@ -94,7 +103,7 @@ void CellComplex::insertCell(Cell* cell)
void CellComplex::removeCell(Cell* cell, bool other)
{
- if(!hasCell(cell)) return;
+ //if(!hasCell(cell)) return;
std::map<Cell*, short int, Less_Cell > coboundary;
cell->getCoboundary(coboundary);
std::map<Cell*, short int, Less_Cell > boundary;
@@ -120,7 +129,8 @@ void CellComplex::removeCell(Cell* cell, bool other)
bdCell->removeCoboundaryCell(cell, false);
}*/
- _cells[cell->getDim()].erase(cell);
+ int erased = _cells[cell->getDim()].erase(cell);
+ if(!erased) Msg::Debug("Tried to remove a cell from the cell complex \n");
if(!_saveorig && !cell->isCombined()) _removedcells.push_back(cell);
}
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index 7968a153e6783192d294fbcad7ba4daf3a17bfe6..c32124083edff2fcad8574e0e9b0344d70eaf17e 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -26,8 +26,8 @@ class CellComplex
private:
GModel* _model;
-
- // sorted containers of unique cells in this cell complex
+
+ // sorted containers of unique cells in this cell complex
// one for each dimension
std::set<Cell*, Less_Cell> _cells[4];
@@ -42,73 +42,74 @@ class CellComplex
int _dim;
bool _saveorig;
- // for constructor
+ // for constructor
bool _insertCells(std::vector<MElement*>& elements, int domain);
-
+
// enqueue cells in queue if they are not there already
- void enqueueCells(std::map<Cell*, short int, Less_Cell>& cells,
+ void enqueueCells(std::map<Cell*, short int, Less_Cell>& cells,
std::queue<Cell*>& Q, std::set<Cell*, Less_Cell>& Qset);
-
+
// insert/remove a cell from this cell complex
void removeCell(Cell* cell, bool other=true);
void insertCell(Cell* cell);
-
+
// queued coreduction
- int coreduction(Cell* startCell, bool omit,
+ int coreduction(Cell* startCell, bool omit,
std::vector<Cell*>& omittedCells);
public:
CellComplex(GModel* model,
- std::vector<MElement*>& domainElements,
- std::vector<MElement*>& subdomainElements);
+ std::vector<MElement*>& domainElements,
+ std::vector<MElement*>& subdomainElements,
+ bool saveOriginalComplex=true);
~CellComplex();
-
+
GModel* getModel() const { return _model; }
- int getDim() { return _dim; }
+ int getDim() { return _dim; }
bool simplicial() { return _simplicial; }
// get the number of certain dimensional cells
- int getSize(int dim, bool orig=false){
+ int getSize(int dim, bool orig=false){
if(!orig) return _cells[dim].size();
else return _ocells[dim].size(); }
-
+
// get dim-dimensional cells
// domain = 0: cells in domain relative to subdomain
// domain = 1: cells in domain
// domain = 2: cells in subdomain
void getCells(std::set<Cell*, Less_Cell>& cells, int dim, int domain=0);
//std::set<Cell*, Less_Cell> getOrigCells(int dim){ return _ocells[dim]; }
-
+
// iterator for the cells of same dimension
typedef std::set<Cell*, Less_Cell>::iterator citer;
-
+
// iterators to the first and last cells of certain dimension
citer firstCell(int dim, bool orig=false) {
return orig ? _ocells[dim].begin() : _cells[dim].begin(); }
citer lastCell(int dim, bool orig=false) {
return orig ? _ocells[dim].end() : _cells[dim].end(); }
-
+
// true if cell complex has given cell
- bool hasCell(Cell* cell, bool orig=false);
-
+ bool hasCell(Cell* cell, bool orig=false);
+
// check whether two cells both belong to subdomain or if neither one does
- bool inSameDomain(Cell* c1, Cell* c2) const {
+ bool inSameDomain(Cell* c1, Cell* c2) const {
return (c1->getDomain() == c2->getDomain()); }
-
+
// remove cells in subdomain from this cell complex
void removeSubdomain();
// (co)reduction of this cell complex
// removes (co)reduction pairs of cell of dimension dim and dim-1
int reduction(int dim, bool omit, std::vector<Cell*>& omittedCells);
- int coreduction(int dim, bool omit, std::vector<Cell*>& omittedCells);
-
+ int coreduction(int dim, bool omit, std::vector<Cell*>& omittedCells);
+
// Cell combining for reduction and coreduction
int combine(int dim);
int cocombine(int dim);
-
- // check whether all boundary cells of a cell has this cell
+
+ // check whether all boundary cells of a cell has this cell
// as coboundary cell and vice versa
// also check whether all (co)boundary cells of a cell
// belong to this cell complex
@@ -118,12 +119,12 @@ class CellComplex
// (with highest dimensional cell omitting?)
int reduceComplex(bool docombine=true, bool omit=true);
int coreduceComplex(bool docombine=true, bool omit=true);
-
- int eulerCharacteristic(){
+
+ int eulerCharacteristic(){
return getSize(0) - getSize(1) + getSize(2) - getSize(3);}
- void printEuler(){
+ void printEuler(){
printf("Euler characteristic: %d. \n", eulerCharacteristic()); }
-
+
// restore the cell complex to its original state before (co)reduction
bool restoreComplex();
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index f516d59953a2cf1593c1be5ef9b8a7499086bc20..69d2047fc6656f88ccc32c4dfac20b89a14a31a9 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -2762,8 +2762,9 @@ void GModel::computeHomology()
std::pair<std::multimap<dpair, std::string>::iterator,
std::multimap<dpair, std::string>::iterator> itp =
_homologyRequests.equal_range(*it);
+ bool prepareToRestore = (itp.first != itp.second);
Homology* homology = new Homology(this, itp.first->first.first,
- itp.first->first.second, false);
+ itp.first->first.second, false, prepareToRestore);
CellComplex *cellcomplex = homology->createCellComplex();
if(cellcomplex->getSize(0)){
for(std::multimap<dpair, std::string>::iterator itt = itp.first;
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index f82c55861a46dbf503540cadd541ecd943f5c99d..ff87c781abaa21bb85a7d03a385a877f52e1e289 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -15,10 +15,12 @@
#if defined(HAVE_KBIPACK)
Homology::Homology(GModel* model, std::vector<int> physicalDomain,
- std::vector<int> physicalSubdomain, bool save0N,
+ std::vector<int> physicalSubdomain,
+ bool save0N, bool saveOrig,
bool combine, bool omit, bool smoothen) :
_model(model), _domain(physicalDomain), _subdomain(physicalSubdomain),
- _save0N(save0N), _combine(combine), _omit(omit), _smoothen(smoothen)
+ _save0N(save0N), _saveOrig(saveOrig),
+ _combine(combine), _omit(omit), _smoothen(smoothen)
{
_fileName = "";
@@ -90,7 +92,8 @@ CellComplex* Homology::createCellComplex(std::vector<GEntity*>& domainEntities,
CellComplex* cellComplex = new CellComplex(_model,
domainElements,
- subdomainElements);
+ subdomainElements,
+ _saveOrig);
if(cellComplex->getSize(0) == 0){
Msg::Error("Cell Complex is empty: check the domain and the mesh");
diff --git a/Geo/Homology.h b/Geo/Homology.h
index fc0962fe7f1d83afc277ac24fc7ddf18a478a44f..054f8309b078e87bde5e05f0af8b2b0873f48a2a 100644
--- a/Geo/Homology.h
+++ b/Geo/Homology.h
@@ -52,6 +52,7 @@ class Homology
// save chains of 0 and highest dimension
bool _save0N;
+ bool _saveOrig;
int _maxdomain;
int _maxnum;
@@ -64,7 +65,8 @@ class Homology
public:
Homology(GModel* model, std::vector<int> physicalDomain,
- std::vector<int> physicalSubdomain, bool save0N=false,
+ std::vector<int> physicalSubdomain,
+ bool save0N=false, bool saveOrig=true,
bool combine=true, bool omit=true, bool smoothen=true);
~Homology();