diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 74abf563ac22c3cca4abe63f683349a9a1dcce1e..87b9526a47cbaec00c0434ada0d352a253cc783a 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -18,6 +18,8 @@ CellComplex::CellComplex(GModel* model,
_model(model), _dim(0), _simplicial(true), _saveorig(saveOriginalComplex),
_relative(false)
{
+ _smallestCell.second = -1.;
+ _biggestCell.second = -1.;
_deleteCount = 0;
_createCount = 0;
_insertCells(subdomainElements, 1);
@@ -62,8 +64,17 @@ CellComplex::CellComplex(GModel* model,
bool CellComplex::_insertCells(std::vector<MElement*>& elements,
int domain)
{
+ std::pair<Cell*, double> smallestElement[4];
+ std::pair<Cell*, double> biggestElement[4];
+ for(int i = 0; i < 4; i++) {
+ smallestElement[i].second = -1.;
+ biggestElement[i].second = -1.;
+ }
+ _dim = 0;
+
for(unsigned int i=0; i < elements.size(); i++){
MElement* element = elements.at(i);
+ int dim = element->getDim();
int type = element->getType();
if(type == TYPE_PYR || type == TYPE_PRI ||
type == TYPE_POLYG || type == TYPE_POLYH) {
@@ -77,6 +88,8 @@ bool CellComplex::_insertCells(std::vector<MElement*>& elements,
delete maybeCell.first;
continue;
}
+
+ if(_dim < dim) _dim = dim;
Cell* cell = maybeCell.first;
std::pair<citer, bool> insert =
_cells[cell->getDim()].insert(cell);
@@ -86,7 +99,17 @@ bool CellComplex::_insertCells(std::vector<MElement*>& elements,
if(domain) cell->setDomain(domain);
}
else _createCount++;
+
+ if(domain == 0) {
+ double size = fabs(element->getVolume());
+ if(smallestElement[dim].second < 0. || smallestElement[dim].second > size)
+ smallestElement[dim] = std::make_pair(cell, size);
+ if(biggestElement[dim].second < 0. || biggestElement[dim].second < size)
+ biggestElement[dim] = std::make_pair(cell, size);
+ }
}
+ _smallestCell = smallestElement[_dim];
+ _biggestCell = biggestElement[_dim];
for (int dim = 3; dim > 0; dim--){
for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
@@ -109,6 +132,10 @@ bool CellComplex::_insertCells(std::vector<MElement*>& elements,
if(domain == 0) {
int ori = cell->findBdCellOrientation(newCell, i);
cell->addBoundaryCell( ori, newCell, true);
+ if(_smallestCell.first == cell)
+ _smallestCell = std::make_pair(newCell, _smallestCell.second);
+ if(_biggestCell.first == cell)
+ _biggestCell = std::make_pair(newCell, _biggestCell.second);
}
}
}
@@ -508,7 +535,7 @@ void CellComplex::removeCells(int dim)
_reduced = true;
}
-int CellComplex::coreduceComplex(int combine, bool omit)
+int CellComplex::coreduceComplex(int combine, bool omit, int heuristic)
{
Msg::Debug("Cell Complex coreduction:");
Msg::Debug(" %d volumes, %d faces, %d edges, and %d vertices",
@@ -537,9 +564,23 @@ int CellComplex::coreduceComplex(int combine, bool omit)
std::vector<Cell*> newCells;
while (getSize(0) != 0){
+
citer cit = firstCell(0);
Cell* cell = *cit;
+ if(heuristic == -1 && _smallestCell.second != 0. &&
+ hasCell(_smallestCell.first)) {
+ Msg::Debug("Omitted a cell in the smallest mesh element with volume %g",
+ _smallestCell.second);
+ cell = _smallestCell.first;
+ }
+ else if(heuristic == 1 && _biggestCell.second != 0. &&
+ hasCell(_biggestCell.first)) {
+ Msg::Debug("Omitted a cell in the biggest mesh element with volume %g",
+ _biggestCell.second);
+ cell = _biggestCell.first;
+ }
+
newCells.push_back(_omitCell(cell, true));
}
for(unsigned int i = 0; i < newCells.size(); i++){
@@ -599,7 +640,7 @@ int CellComplex::combine(int dim)
inSameDomain(s, c1) && inSameDomain(s, c2) &&
c1->getImmune() == c2->getImmune() ) {
- removeCell(s, true, true);
+ removeCell(s, true, false);
c1->getBoundary(bd_c);
enqueueCells(bd_c, Q, Qset);
@@ -615,10 +656,6 @@ int CellComplex::combine(int dim)
cit = firstCell(dim);
count++;
- if(s->isCombined()) {
- delete s;
- _deleteCount++;
- }
if(c1->isCombined()) {
delete c1;
_deleteCount++;
@@ -656,6 +693,7 @@ int CellComplex::cocombine(int dim)
for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
Cell* cell = *cit;
+
cell->getCoboundary(cbd_c);
enqueueCells(cbd_c, Q, Qset);
@@ -674,7 +712,8 @@ int CellComplex::cocombine(int dim)
if(!(*c1 == *c2) && abs(or1) == abs(or2)
&& inSameDomain(s, c1) && inSameDomain(s, c2)){
- removeCell(s, true, true);
+
+ removeCell(s, true, false);
c1->getCoboundary(cbd_c);
enqueueCells(cbd_c, Q, Qset);
@@ -691,10 +730,6 @@ int CellComplex::cocombine(int dim)
cit = firstCell(dim);
count++;
- if(s->isCombined()) {
- delete s;
- _deleteCount++;
- }
if(c1->isCombined()) {
delete c1;
_deleteCount++;
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index d8b97eb8b3f788655f97c10faa4e92f08bfd3a83..e75c8c3d65ce4b0ab83605e0e4b2b29d715a1057 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -25,6 +25,10 @@ class CellComplex
{
private:
+
+ std::pair<Cell*, double> _smallestCell;
+ std::pair<Cell*, double> _biggestCell;
+
GModel* _model;
// sorted containers of unique cells in this cell complex
@@ -158,8 +162,11 @@ class CellComplex
// (combine = 1 -> with combining)
// (omit = true -> with highest dimensional cell omitting?)
// (homseq = true -> homology sequence splitting possible after reduction)
+ // (heuristic = 0 -> no heuristic, let mesh indexing determine
+ // 1 -> omit 0-cell in biggest element
+ // -1 -> omit 0-cell in smallest element)
int reduceComplex(int combine=1, bool omit=true, bool homseq=false);
- int coreduceComplex(int combine=1, bool omit=true);
+ int coreduceComplex(int combine=1, bool omit=true, int heuristic=0);
bool isReduced() const { return _reduced; }
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index 19c943ac721d9f0420b7bcf24b399d3895b44c3f..6ec9f1e218e4e7aeff96db8360dfc28ae43d9938 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -19,10 +19,11 @@ Homology::Homology(GModel* model,
std::vector<int> physicalSubdomain,
std::vector<int> physicalImdomain,
bool saveOrig,
- int combine, bool omit, bool smoothen) :
+ int combine, bool omit, bool smoothen, int heuristic) :
_model(model), _domain(physicalDomain), _subdomain(physicalSubdomain),
_imdomain(physicalImdomain), _saveOrig(saveOrig),
- _combine(combine), _omit(omit), _smoothen(smoothen), _cellComplex(NULL)
+ _combine(combine), _omit(omit), _smoothen(smoothen),
+ _heuristic(heuristic), _cellComplex(NULL)
{
_fileName = "";
@@ -49,6 +50,9 @@ Homology::Homology(GModel* model,
_cohomologyComputed[i] = false;
_betti[i] = -1;
}
+
+ if(abs(_heuristic) > 1) _heuristic = 0;
+
}
void Homology::_getEntities(const std::vector<int>& physicalGroups,
@@ -257,7 +261,7 @@ void Homology::findCohomologyBasis(std::vector<int> dim)
double t1 = Cpu();
- int omitted = _cellComplex->coreduceComplex(_combine, _omit);
+ int omitted = _cellComplex->coreduceComplex(_combine, _omit, _heuristic);
std::sort(dim.begin(), dim.end());
if(!dim.empty() && _combine > 1) {
diff --git a/Geo/Homology.h b/Geo/Homology.h
index 386e6a7e7f122d448a4d6588b1b462079de87a62..99241b0292ff22dddbe0f367bb160a18d4b8f86f 100644
--- a/Geo/Homology.h
+++ b/Geo/Homology.h
@@ -49,6 +49,8 @@ class Homology
bool _omit;
// use chain smoothning
bool _smoothen;
+ // corecution heuristic
+ int _heuristic;
// file name to store the results
std::string _fileName;
@@ -93,7 +95,8 @@ class Homology
std::vector<int> physicalSubdomain,
std::vector<int> physicalIm,
bool saveOrig=true,
- int combine=2, bool omit=true, bool smoothen=true);
+ int combine=2, bool omit=true, bool smoothen=true,
+ int heuristic=1);
~Homology();
GModel* getModel() const { return _model; }
diff --git a/Plugin/HomologyComputation.cpp b/Plugin/HomologyComputation.cpp
index e3794c78d4521336fdb3f6daf6368d0b306d2e27..24dbb8e7db03fe3f87d1c334fdd5f45336065661 100644
--- a/Plugin/HomologyComputation.cpp
+++ b/Plugin/HomologyComputation.cpp
@@ -25,7 +25,8 @@ StringXNumber HomologyComputationOptions_Number[] = {
{GMSH_FULLRC, "CreatePostProcessingViews", NULL, 1.},
{GMSH_FULLRC, "ReductionOmit", NULL, 1.},
{GMSH_FULLRC, "ReductionCombine", NULL, 2.},
- {GMSH_FULLRC, "PostProcessSimplify", NULL, 1.}
+ {GMSH_FULLRC, "PostProcessSimplify", NULL, 1.},
+ {GMSH_FULLRC, "ReductionHeuristic", NULL, 1.}
};
StringXString HomologyComputationOptions_String[] = {
@@ -111,6 +112,7 @@ PView *GMSH_HomologyComputationPlugin::execute(PView *v)
bool omit = (bool)HomologyComputationOptions_Number[5].def;
int combine = (int)HomologyComputationOptions_Number[6].def;
bool smoothen = (bool)HomologyComputationOptions_Number[7].def;
+ int heuristic = (int)HomologyComputationOptions_Number[8].def;
std::vector<int> domain;
std::vector<int> subdomain;
@@ -124,7 +126,7 @@ PView *GMSH_HomologyComputationPlugin::execute(PView *v)
GModel* m = GModel::current();
Homology* homology = new Homology(m, domain, subdomain, imdomain,
- true, combine, omit, smoothen);
+ true, combine, omit, smoothen, heuristic);
homology->setFileName(fileName);
if(hom != 0) homology->findHomologyBasis(dimsave);