diff --git a/Geo/Cell.cpp b/Geo/Cell.cpp
index 3b7895cbdd6e69e70de950ae538f42119208fb43..d8c773c4f8b7eae18272cf7b1c070907c1082187 100644
--- a/Geo/Cell.cpp
+++ b/Geo/Cell.cpp
@@ -374,7 +374,18 @@ void Cell::printCell()
printf("combined: %d. \n" , isCombined() );
};
-void Cell::restoreCell(){
+void Cell::saveCellBoundary()
+{
+ for(biter it = firstCoboundary(); it != lastCoboundary(); it++){
+ it->second.init();
+ }
+ for(biter it = firstBoundary(); it != lastBoundary(); it++){
+ it->second.init();
+ }
+}
+
+void Cell::restoreCellBoundary()
+{
std::vector<Cell*> toRemove;
for(biter it = firstCoboundary(true); it != lastCoboundary(); it++){
it->second.reset();
@@ -473,6 +484,74 @@ bool Cell::hasCoboundary(Cell* cell, bool orig)
}
}
+Cell::biter Cell::firstBoundary(bool orig)
+{
+ biter it = _bd.begin();
+ if(!orig) while(it->second.get() == 0 && it != _bd.end()) it++;
+ else while(it->second.geto() == 0 && it != _bd.end()) it++;
+ return it;
+}
+
+Cell::biter Cell::lastBoundary()
+{
+ return _bd.end();
+}
+
+Cell::biter Cell::firstCoboundary(bool orig)
+{
+ biter it = _cbd.begin();
+ if(!orig) while(it->second.get() == 0 && it != _cbd.end()) it++;
+ else while(it->second.geto() == 0 && it != _cbd.end()) it++;
+ return it;
+}
+
+Cell::biter Cell::lastCoboundary()
+{
+ return _cbd.end();
+}
+
+int Cell::getBoundarySize(bool orig)
+{
+ int size = 0;
+ for(biter bit = _bd.begin(); bit != _bd.end(); bit++){
+ if(!orig && bit->second.get() != 0) size++;
+ else if(orig && bit->second.geto() != 0) size++;
+ }
+ return size;
+}
+
+int Cell::getCoboundarySize(bool orig)
+{
+ int size = 0;
+ for(biter bit = _cbd.begin(); bit != _cbd.end(); bit++){
+ if(!orig && bit->second.get() != 0) size++;
+ else if(orig && bit->second.geto() != 0) size++;
+ }
+ return size;
+}
+
+void Cell::getBoundary(std::map<Cell*, short int, Less_Cell>& boundary,
+ bool orig)
+{
+ boundary.clear();
+ for(biter it = firstBoundary(); it != lastBoundary(); it++){
+ Cell* cell = it->first;
+ if(!orig && it->second.get() != 0) boundary[cell] = it->second.get();
+ if(orig && it->second.geto() != 0) boundary[cell] = it->second.geto();
+ }
+}
+
+void Cell::getCoboundary(std::map<Cell*, short int, Less_Cell>& coboundary,
+ bool orig)
+{
+ coboundary.clear();
+ for(biter it = firstCoboundary(); it != lastCoboundary(); it++){
+ Cell* cell = it->first;
+ if(!orig && it->second.get() != 0) coboundary[cell] = it->second.get();
+ if(orig && it->second.geto() != 0) coboundary[cell] = it->second.geto();
+ }
+}
+
void Cell::printBoundary()
{
for(biter it = firstBoundary(); it != lastBoundary(); it++){
diff --git a/Geo/Cell.h b/Geo/Cell.h
index 62d4fd6f8b0da76ec977eb9f5f8b084c00aa56c3..5c999d81c829665303af075ec10c848d9e581281 100644
--- a/Geo/Cell.h
+++ b/Geo/Cell.h
@@ -22,17 +22,16 @@ public:
// Class to save cell boundary orientation information
class BdInfo {
private:
- signed char _ori;
- signed char _origOri;
+ signed char _ori[2];
public:
- BdInfo(int ori) { _ori = ori; _origOri = 0; }
+ BdInfo(int ori) { _ori[0] = ori; _ori[1] = 0; }
- int get() const { return _ori; }
- void reset() { _ori = _origOri; }
- void init() { _origOri = _ori; }
- void set(int ori) { _ori = ori; }
- int geto() const { return _origOri; }
+ int get() const { return _ori[0]; }
+ void reset() { _ori[0] = _ori[1]; }
+ void init() { _ori[1] = _ori[0]; }
+ void set(int ori) { _ori[0] = ori; }
+ int geto() const { return _ori[1]; }
};
@@ -100,81 +99,37 @@ class Cell {
void increaseGlobalNum() { _globalNum++; }
- // restores the cell information to its original state before reduction
- void restoreCell();
+ // save/restore the original boundary information of the cell
+ void saveCellBoundary();
+ void restoreCellBoundary();
// true if this cell has given vertex
virtual bool hasVertex(int vertex) const;
// (co)boundary cell iterator
typedef std::map<Cell*, BdInfo, Less_Cell>::iterator biter;
+
// iterators to (first/last (co)boundary cells of this cell
// (orig: to original (co)boundary cells of this cell)
- biter firstBoundary(bool orig=false){
- biter it = _bd.begin();
- if(!orig) while(it->second.get() == 0 && it != _bd.end()) it++;
- else while(it->second.geto() == 0 && it != _bd.end()) it++;
- return it;
- }
- biter lastBoundary(){ return _bd.end(); }
- biter firstCoboundary(bool orig=false){
- biter it = _cbd.begin();
- if(!orig) while(it->second.get() == 0 && it != _cbd.end()) it++;
- else while(it->second.geto() == 0 && it != _cbd.end()) it++;
- return it;
- }
- biter lastCoboundary(){ return _cbd.end(); }
-
- int getBoundarySize(bool orig=false) {
- int size = 0;
- for(biter bit = _bd.begin(); bit != _bd.end(); bit++){
- if(!orig && bit->second.get() != 0) size++;
- else if(orig && bit->second.geto() != 0) size++;
- }
- return size; }
- int getCoboundarySize(bool orig=false) {
- int size = 0;
- for(biter bit = _cbd.begin(); bit != _cbd.end(); bit++){
- if(!orig && bit->second.get() != 0) size++;
- else if(orig && bit->second.geto() != 0) size++;
- }
- return size; }
+ biter firstBoundary(bool orig=false);
+ biter lastBoundary();
+ biter firstCoboundary(bool orig=false);
+ biter lastCoboundary();
+
+ int getBoundarySize(bool orig=false);
+ int getCoboundarySize(bool orig=false);
// get the (orig: original) cell boundary
void getBoundary(std::map<Cell*, short int, Less_Cell>& boundary,
- bool orig=false){
- boundary.clear();
- for(biter it = firstBoundary(); it != lastBoundary(); it++){
- Cell* cell = it->first;
- if(!orig && it->second.get() != 0) boundary[cell] = it->second.get();
- if(orig && it->second.geto() != 0) boundary[cell] = it->second.geto();
- }
- }
+ bool orig=false);
void getCoboundary(std::map<Cell*, short int, Less_Cell>& coboundary,
- bool orig=false){
- coboundary.clear();
- for(biter it = firstCoboundary(); it != lastCoboundary(); it++){
- Cell* cell = it->first;
- if(!orig && it->second.get() != 0) coboundary[cell] = it->second.get();
- if(orig && it->second.geto() != 0) coboundary[cell] = it->second.geto();
- }
- }
-
+ bool orig=false);
// add (co)boundary cell
// (other: reciprocally also add this cell from the other cell's (co)boundary)
void addBoundaryCell(int orientation, Cell* cell, bool other);
void addCoboundaryCell(int orientation, Cell* cell, bool other);
- void saveOriginalBd() {
- for(biter it = firstCoboundary(); it != lastCoboundary(); it++){
- it->second.init();
- }
- for(biter it = firstBoundary(); it != lastBoundary(); it++){
- it->second.init();
- }
- }
-
// remove (co)boundary cell
// (other: reciprocally also revove this cell from the other cell's (co)boundary)
void removeBoundaryCell(Cell* cell, bool other);
diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 29b31d0bc4e4f8ae3271af1c5a105442d1a52ea6..f08f804863698311127a32b1e386b97ba26ebe7c 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -44,7 +44,7 @@ CellComplex::CellComplex(GModel* model,
Cell* cell = *cit;
cell->setNum(++num);
cell->increaseGlobalNum();
- cell->saveOriginalBd();
+ cell->saveCellBoundary();
}
}
@@ -316,7 +316,8 @@ int CellComplex::coreduction(Cell* startCell, int omit,
Q.pop();
Qset.erase(s);
if(s->getBoundarySize() == 1 &&
- inSameDomain(s, s->firstBoundary()->first)){
+ inSameDomain(s, s->firstBoundary()->first) &&
+ abs(s->firstBoundary()->second.get()) < 2){
s->getBoundary(bd_s);
removeCell(s);
bd_s.begin()->first->getCoboundary(cbd_c);
@@ -341,6 +342,9 @@ int CellComplex::reduction(int dim, int omit,
{
if(dim < 1 || dim > 3) return 0;
+ int numCells[4];
+ for(int i = 0; i < 4; i++) numCells[i] = getSize(i);
+
int count = 0;
bool reduced = true;
@@ -353,7 +357,8 @@ int CellComplex::reduction(int dim, int omit,
if(cell->getCoboundarySize() == 1 &&
inSameDomain(cell, cell->firstCoboundary()->first) &&
!cell->getImmune() &&
- !cell->firstCoboundary()->first->getImmune()){
+ !cell->firstCoboundary()->first->getImmune() &&
+ abs(cell->firstCoboundary()->second.get()) < 2){
cit++;
if(dim == omit){
omittedCells.push_back(cell->firstCoboundary()->first);
@@ -369,6 +374,9 @@ int CellComplex::reduction(int dim, int omit,
}
}
_reduced = true;
+ Msg::Debug("Cell complex %d-reduction removed %dv, %df, %de, %dn", dim,
+ numCells[3]-getSize(3), numCells[2]-getSize(2),
+ numCells[1]-getSize(1), numCells[0]-getSize(0));
return count;
}
@@ -377,6 +385,9 @@ int CellComplex::coreduction(int dim, int omit,
{
if(dim < 1 || dim > 3) return 0;
+ int numCells[4];
+ for(int i = 0; i < 4; i++) numCells[i] = getSize(i);
+
int count = 0;
bool reduced = true;
@@ -387,7 +398,8 @@ int CellComplex::coreduction(int dim, int omit,
while(cit != lastCell(dim)){
Cell* cell = *cit;
if(cell->getBoundarySize() == 1 &&
- inSameDomain(cell, cell->firstBoundary()->first)) {
+ inSameDomain(cell, cell->firstBoundary()->first) &&
+ abs(cell->firstBoundary()->second.get()) < 2) {
++cit;
if(dim-1 == omit){
omittedCells.push_back(cell->firstBoundary()->first);
@@ -403,6 +415,9 @@ int CellComplex::coreduction(int dim, int omit,
}
}
_reduced = true;
+ Msg::Debug("Cell complex %d-coreduction removed %dv, %df, %de, %dn", dim,
+ numCells[3]-getSize(3), numCells[2]-getSize(2),
+ numCells[1]-getSize(1), numCells[0]-getSize(0));
return count;
}
@@ -428,11 +443,15 @@ int CellComplex::getDomain(Cell* cell, std::string& str)
Cell* CellComplex::_omitCell(Cell* cell, bool dual)
{
+ Msg::Debug("Omitting %d-cell from the cell complex", cell->getDim());
removeCell(cell, false);
std::vector<Cell*> omittedCells;
omittedCells.push_back(cell);
int count = 0;
+ int numCells[4];
+ for(int i = 0; i < 4; i++) numCells[i] = getSize(i);
+
if(!dual) {
for(int j = 3; j > 0; j--)
count += reduction(j, cell->getDim(), omittedCells);
@@ -449,8 +468,10 @@ Cell* CellComplex::_omitCell(Cell* cell, bool dual)
std::string domainstr = "";
int domain = getDomain(cell, domainstr);
- Msg::Debug("Omitted %d-cell in %s that caused %d reductions",
- cell->getDim(), domainstr.c_str(), count);
+ Msg::Debug("Cell complex %d-omit removed %dv, %df, %de, %dn",
+ cell->getDim(),
+ numCells[3]-getSize(3), numCells[2]-getSize(2),
+ numCells[1]-getSize(1), numCells[0]-getSize(0));
Msg::Debug(" - number of %d-cells left in %s: %d",
cell->getDim(), domainstr.c_str(),
getNumCells(cell->getDim(), domain));
@@ -460,10 +481,6 @@ Cell* CellComplex::_omitCell(Cell* cell, bool dual)
int CellComplex::reduceComplex(int combine, bool omit, bool homseq)
{
- Msg::Debug("Cell Complex reduction:");
- Msg::Debug(" %d volumes, %d faces, %d edges, and %d vertices",
- getSize(3), getSize(2), getSize(1), getSize(0));
-
if(!getSize(0)) return 0;
double t1 = Cpu();
@@ -510,9 +527,6 @@ int CellComplex::reduceComplex(int combine, bool omit, bool homseq)
if(combine > 2) for(int i = 3; i > 0; i--) reduction(i, -1, empty);
else if(combine > 1) reduction(0, -1, empty);
- Msg::Debug(" %d volumes, %d faces, %d edges, and %d vertices",
- getSize(3), getSize(2), getSize(1), getSize(0));
-
_reduced = true;
return count;
}
@@ -543,10 +557,6 @@ void CellComplex::removeCells(int dim)
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",
- getSize(3), getSize(2), getSize(1), getSize(0));
-
if(!getSize(0)) return 0;
double t1 = Cpu();
@@ -618,8 +628,6 @@ int CellComplex::coreduceComplex(int combine, bool omit, int heuristic)
else if(combine > 1) coreduction(3, -1, empty);
coherent();
- Msg::Debug(" %d volumes, %d faces, %d edges, and %d vertices",
- getSize(3), getSize(2), getSize(1), getSize(0));
_reduced = true;
return count;
@@ -682,11 +690,11 @@ std::vector<int> CellComplex::bettiCoreduceComplex()
int CellComplex::combine(int dim)
{
- //Msg::Debug("Cell complex before combining:");
- //Msg::Debug(" %d volumes, %d faces, %d edges, and %d vertices",
- // getSize(3), getSize(2), getSize(1), getSize(0));
if(dim < 1 || dim > 3) return 0;
+ int numCells[4];
+ for(int i = 0; i < 4; i++) numCells[i] = getSize(i);
+
double t1 = Cpu();
std::queue<Cell*> Q;
@@ -755,21 +763,20 @@ int CellComplex::combine(int dim)
}
}
- //Msg::Debug("Cell complex after combining:");
- //Msg::Debug(" %d volumes, %d faces, %d edges, and %d vertices",
- // getSize(3), getSize(2), getSize(1), getSize(0));
+ Msg::Debug("Cell complex %d-combine removed %dv, %df, %de, %dn", dim,
+ numCells[3]-getSize(3), numCells[2]-getSize(2),
+ numCells[1]-getSize(1), numCells[0]-getSize(0));
_reduced = true;
return count;
}
int CellComplex::cocombine(int dim)
{
- //Msg::Debug("Cell complex before cocombining:");
- //Msg::Debug(" %d volumes, %d faces, %d edges, and %d vertices",
- //getSize(3), getSize(2), getSize(1), getSize(0));
-
if(dim < 0 || dim > 2) return 0;
+ int numCells[4];
+ for(int i = 0; i < 4; i++) numCells[i] = getSize(i);
+
double t1 = Cpu();
std::queue<Cell*> Q;
@@ -839,9 +846,9 @@ int CellComplex::cocombine(int dim)
}
}
- //Msg::Debug("Cell complex after cocombining:");
- //Msg::Debug(" %d volumes, %d faces, %d edges, and %d vertices",
- // getSize(3), getSize(2), getSize(1), getSize(0));
+ Msg::Debug("Cell complex %d-cocombine removed %dv, %df, %de, %dn", dim,
+ numCells[3]-getSize(3), numCells[2]-getSize(2),
+ numCells[1]-getSize(1), numCells[0]-getSize(0));
_reduced = true;
return count;
}
@@ -973,7 +980,7 @@ bool CellComplex::restoreComplex()
_cells[i] = _ocells[i];
for(citer cit = firstCell(i); cit != lastCell(i); cit++){
Cell* cell = *cit;
- cell->restoreCell();
+ cell->restoreCellBoundary();
if(relative()) {
if(cell->inSubdomain()) _numSubdomainCells[i] += 1;
else _numRelativeCells[i] += 1;