From 5bed86a1d4dfd3c1fbb3e12ebe086e39cfda20be Mon Sep 17 00:00:00 2001
From: Matti Pellika <matti.pellikka@tut.fi>
Date: Sun, 3 Jan 2010 16:28:59 +0000
Subject: [PATCH] Rewrote chain local deformation.
---
Geo/CellComplex.cpp | 5 +
Geo/CellComplex.h | 5 +
Geo/ChainComplex.cpp | 352 ++++++++-----------------------------------
Geo/ChainComplex.h | 11 +-
4 files changed, 79 insertions(+), 294 deletions(-)
diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index ccab617075..a0fe39fb5f 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -16,6 +16,7 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su
_subdomain = subdomain;
_dim = 0;
+ _simplicial = true;
// find boundary elements
find_boundary(domain, subdomain);
@@ -167,16 +168,20 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){
}
else if(type == MSH_QUA_4 || type == MSH_QUA_8 || type == MSH_QUA_9){
cell = new CQuadrangle(vertices, tag, subdomain, boundary);
+ _simplicial = false;
}
/*
else if(type == MSH_HEX_8 || type == MSH_HEX_27 || type == MSH_HEX_20){
cell = new CHexahedron(vertices, tag, subdomain, boundary);
+ _simplicial = false;
}
else if(type == MSH_PRI_6 || type == MSH_PRI_18 || type == MSH_PRI_15){
cell = new CPrism(vertices, tag, 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);
+ _simplicial = false;
}*/
else printf("Error: mesh element %d not implemented yet! \n", type);
tag++;
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index 41a9bb0e85..da55e8a68a 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -903,6 +903,9 @@ class CellComplex
int _dim;
+ // Is the cell complex simplicial
+ bool _simplicial;
+
public:
// iterator for the cells of same dimension
typedef std::set<Cell*, Less_Cell>::iterator citer;
@@ -998,6 +1001,8 @@ class CellComplex
int getDim() {return _dim; }
+ bool simplicial() { return _simplicial; }
+
std::set<Cell*, Less_Cell> getCells(int dim){ return _cells[dim]; }
std::set<Cell*, Less_Cell> getOrgCells(int dim){ return _ocells[dim]; }
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index 0b8e2ac284..90c2c9f7fc 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -483,325 +483,104 @@ Chain::Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs, CellComp
}
-Cell* Chain::findCommonCbdCell(Cell* c1, Cell* c2, Cell* c3){
- /*
- c1->printCell();
- c1->printOrgCbd();
- c2->printOrgCbd();
- c2->printCell();
- printf("------ \n");
- */
-
- std::map< Cell*, int, Less_Cell > c1Cbd = c1->getOrgCbd();
- for(std::map< Cell*, int, Less_Cell >::iterator it = c1Cbd.begin(); it != c1Cbd.end(); it++){
- Cell* c1CbdCell = (*it).first;
- if(c3 == NULL){
- if(c2->hasCoboundary(c1CbdCell, true)) return c1CbdCell;
- }
- else{
- if(c2->hasCoboundary(c1CbdCell, true) && c3->hasCoboundary(c1CbdCell, true)) return c1CbdCell;
- }
- }
+bool Chain::deform(std::map<Cell*, int, Less_Cell> &cellsInChain, std::map<Cell*, int, Less_Cell> &cellsNotInChain){
- return NULL;
-
-}
-
-std::pair<Cell*, int> Chain::findRemainingBoundary(Cell* b, Cell* c1, Cell* c2, Cell* c3){
- std::map<Cell*, int, Less_Cell> bBd = b->getOrgBd();
- for(std::map<Cell*, int, Less_Cell >::iterator it = bBd.begin(); it != bBd.end(); it++){
- Cell* bBdCell = (*it).first;
- if(c3 == NULL) {
- if( !(*bBdCell == *c1) && !(*bBdCell == *c2)) return *it;
- }
- else{
- if( !(*bBdCell == *c1) && !(*bBdCell == *c2) && !(*bBdCell == *c3) ) return *it;
+ std::vector<int> cc;
+ std::vector<int> bc;
+
+ for(citer cit = cellsInChain.begin(); cit != cellsInChain.end(); cit++){
+ Cell* c = (*cit).first;
+ c->setImmune(false);
+ if(!c->inSubdomain()) {
+ cc.push_back(getCoeff(c));
+ bc.push_back((*cit).second);
}
+ removeCell(c);
}
- return std::make_pair(b, 0);
-}
-
-int Chain::findOrientation(Cell* b, Cell* c){
- std::map< Cell*, int, Less_Cell > bBd = b->getOrgBd();
- std::map< Cell*, int, Less_Cell >::iterator it = bBd.find(c);
- if(it != bBd.end()) return (*it).second;
- return 0;
-}
-
-std::map<Cell*, int, Less_Cell> Chain::getBdCellsInChain(Cell* cell){
-
- std::map<Cell*, int, Less_Cell> cells;
- std::map<Cell*, int, Less_Cell> boundary = cell->getOrgBd();
- for(citer cit = boundary.begin(); cit != boundary.end(); cit++){
- Cell* BdCell = (*cit).first;
- int BdCellO = (*cit).second;
- if(hasCell(BdCell) || BdCell->inSubdomain() ) cells.insert( std::make_pair(BdCell, BdCellO) );
+ // FIXME: orientations don't get right on 2-chains
+ int flip = 1;
+ if(cc[0] != bc[0]) flip = flip*-1;
+
+ int n = 1;
+ for(citer cit = cellsNotInChain.begin(); cit != cellsNotInChain.end(); cit++){
+ Cell* c = (*cit).first;
+ if(n == 2) c->setImmune(true);
+ else c->setImmune(false);
+ int coeff = -1*(*cit).second*flip;
+ addCell(c, coeff);
+ n++;
}
- return cells;
+
+ return true;
}
-bool Chain::removeBoundary( std::pair<Cell*, int> cell ){
+bool Chain::deformChain(std::pair<Cell*, int> cell, bool straighten, bool bend){
Cell* c1 = cell.first;
- int c1c = cell.second;
- if(c1c == 0) return false;
std::map<Cell*, int, Less_Cell> c1Cbd = c1->getOrgCbd();
for(citer cit = c1Cbd.begin(); cit != c1Cbd.end(); cit++){
- Cell* c1CbdCell = (*cit).first;
- std::map<Cell*, int, Less_Cell> cells = getBdCellsInChain(c1CbdCell);
+ std::map<Cell*, int, Less_Cell> cellsInChain;
+ std::map<Cell*, int, Less_Cell> cellsNotInChain;
- if( (getDim() == 1 && cells.size() == 3) || (getDim() == 2 && cells.size() == 4)){
- for(citer cit2 = cells.begin(); cit2 != cells.end(); cit2++){
- Cell* cell = (*cit2).first;
- removeCell(cell);
- }
- return true;
- }
- }
- return false;
-}
-
-bool Chain::straightenChain( std::pair<Cell*, int> cell ){
-
- Cell* c1 = cell.first;
- int c1c = cell.second;
- int c1o = 0;
- if(c1c == 0 || c1->getImmune()) return false;
-
- if(cell.second == 0 || cell.first->getImmune()) return false;
- /*
- Cell* c1 = NULL;
- int c1c = 0;
- int c1o = 0;
- */
- Cell* c2 = NULL;
- int c2c = 0;
- int c2o = 0;
-
- Cell* c3 = NULL;
- int c3c = 0;
- int c3o = 0;
-
- Cell* b = NULL;
-
- std::map<Cell*, int, Less_Cell> c1Cbd = cell.first->getOrgCbd();
- for(citer cit = c1Cbd.begin(); cit != c1Cbd.end(); cit++){
Cell* c1CbdCell = (*cit).first;
- c1o = (*cit).second;
-
- /*
- std::map<Cell*, int, Less_Cell> cells = getBdCellsInChain(c1CbdCell);
- //if((getDim() == 1 && cells.size() != 2) || (getDim() == 2 && cells.size() != 3) ) continue;
- if( getDim() == 1 && cells.size() == 2){
-
- citer cit2 = cells.end();
- c1 = (*cit2).first;
- //c1->printCell();
-
- c1c = getCoeff(c1);
- c1o = (*cit2).second;
-
- cit2 = cells.begin();
- c2 = (*cit2).first;
- c2c = getCoeff(c2);
- c2o = (*cit2).second;
- b = c1CbdCell;
- printf("kakak4 \n");
- }
- else if( getDim() == 2 && cells.size() == 3){
-
- }
- */
-
-
std::map<Cell*, int, Less_Cell> c1CbdBd = c1CbdCell->getOrgBd();
- int count = 0;
for(citer cit2 = c1CbdBd.begin(); cit2 != c1CbdBd.end(); cit2++){
Cell* c1CbdBdCell = (*cit2).first;
- int c1CbdBdCellO = (*cit2).second;
- int coeff = getCoeff(c1CbdBdCell);
- if( (coeff != 0 || c1CbdBdCell->inSubdomain() ) && !(*c1CbdBdCell == *c1) && !c1CbdBdCell->getImmune()){
- if(c1->getDim() == 1){
- count++;
- c2 = c1CbdBdCell; c2c = coeff; c2o = c1CbdBdCellO;
- b = c1CbdCell; break;
- }
- else if(c1->getDim() == 2){
- count++;
- if(count == 1) { c2 = c1CbdBdCell; c2c = coeff; c2o = c1CbdBdCellO; }
- else if(count == 2) { c3 = c1CbdBdCell; c3c = coeff; c3o = c1CbdBdCellO; b = c1CbdCell; break;}
- }
- }
+ int coeff = (*cit2).second;
+ if( (hasCell(c1CbdBdCell) && getCoeff(c1CbdBdCell) != 0) || c1CbdBdCell->inSubdomain()) cellsInChain.insert(std::make_pair(c1CbdBdCell, coeff));
+ else cellsNotInChain.insert(std::make_pair(c1CbdBdCell, coeff));
}
- if (b != NULL) break;
+ bool next = false;
- }
-
- if(c1->getDim() == 1 &&
- b != NULL && c2 != NULL && !(*c2 == *c1)){
-
- int temp1 = c1c - c1o;
-
- std::pair<Cell*, int> c4p = std::make_pair(b, 0);
- c4p = findRemainingBoundary(b, c1, c2);
- Cell* c4 = c4p.first;
- int c4o = c4p.second;
+ for(citer cit2 = cellsInChain.begin(); cit2 != cellsInChain.end(); cit2++){
+ Cell* c = (*cit2).first;
+ int coeff = getCoeff(c);
+ if(c->getImmune()) next = true;
+ if(bend && c->inSubdomain()) next = true;
+ if(coeff > 1 || coeff < -1) next = true;
+ }
- if(c4o != 0 && !c2->getImmune() && !c4->getImmune()
- && ( hasCell(c1) || c1->inSubdomain() ) && (hasCell(c2) || c2->inSubdomain() ) && !hasCell(c4) ){
-
- int c4c = -c4o;
- if(temp1 != 0) c4c= c4c*-1;
-
- this->removeCell(c1);
- this->removeCell(c2);
- c1->setImmune(false);
- c2->setImmune(false);
- c4->setImmune(false);
- if(!c4->inSubdomain()) this->addCell(c4, c4c);
- return true;
+ for(citer cit2 = cellsNotInChain.begin(); cit2 != cellsNotInChain.end(); cit2++){
+ Cell* c = (*cit2).first;
+ if(c->inSubdomain()) next = true;
}
- }
-
- else if(c1->getDim() == 2 &&
- b != NULL && c2 != NULL && c3 != NULL && !(*c2 == *c1) && !(*c1 == *c3) && !(*c2 == *c3)){
- int temp1 = c1c - c1o;
+ if(next) continue;
- std::pair<Cell*, int> c4p = std::make_pair(b, 0);
- c4p = findRemainingBoundary(b, c1, c2, c3);
- Cell* c4 = c4p.first;
- int c4o = c4p.second;
+ //printf("dim: %d, in chain: %d, not in chain: %d \n", getDim(), cellsInChain.size(), cellsNotInChain.size());
- if(c4o != 0 && !c2->getImmune() && !c3->getImmune() && !c4->getImmune()
- && (hasCell(c1) || c1->inSubdomain()) && (hasCell(c2) || c2->inSubdomain())
- && (hasCell(c3) || c3->inSubdomain()) && !hasCell(c4)) {
-
- int c4c = -c4o;
- if(temp1 != 0) c4c= c4c*-1;
-
- this->removeCell(c1);
- this->removeCell(c2);
- this->removeCell(c3);
- c1->setImmune(false);
- c2->setImmune(false);
- c3->setImmune(false);
- c4->setImmune(false);
- if(!c4->inSubdomain()) this->addCell(c4, c4c);
+ if( (getDim() == 1 && cellsInChain.size() == 2 && cellsNotInChain.size() == 1 && straighten) ||
+ (getDim() == 2 && cellsInChain.size() == 3 && cellsNotInChain.size() == 1 && straighten)){
+ //printf("straighten \n");
+ deform(cellsInChain, cellsNotInChain);
return true;
}
-
- }
- return false;
-}
-
-bool Chain::bendChain( std::pair<Cell*, int> cell ){
-
- Cell* c1 = cell.first;
- int c1c = cell.second;
- if(c1c == 0 || c1->getImmune()) return false;
- int c1o = 0;
-
- Cell* c2 = NULL;
- int c2c = 0;
- int c2o = 0;
-
- Cell* c3 = NULL;
- int c3c = 0;
- int c3o = 0;
-
- Cell* b = NULL;
-
- std::map<Cell*, int, Less_Cell> c1Cbd = c1->getOrgCbd();
- for(citer cit = c1Cbd.begin(); cit != c1Cbd.end(); cit++){
- Cell* c1CbdCell = (*cit).first;
- c1o = (*cit).second;
- std::map<Cell*, int, Less_Cell> c1CbdBd = c1CbdCell->getOrgBd();
-
- int count = 0;
- for(citer cit2 = c1CbdBd.begin(); cit2 != c1CbdBd.end(); cit2++){
- Cell* c1CbdBdCell = (*cit2).first;
- int c1CbdBdCellO = (*cit2).second;
- if(!hasCell(c1CbdBdCell) && !c1CbdBdCell->getImmune() ){
- count++;
- if(count == 1) { c2 = c1CbdBdCell; c2o = c1CbdBdCellO; }
- else if(count == 2) { c3 = c1CbdBdCell; c3o = c1CbdBdCellO; b = c1CbdCell; break;}
- }
- }
-
- if (b != NULL) break;
- else c2 = NULL;
- }
-
- if(c1->getDim() == 2 &&
- b != NULL && c2 != NULL && c3 != NULL && !(*c2 == *c1) && !(*c1 == *c3) && !(*c2 == *c3) &&
- (c1c == 1 || c1c == -1) && !c2->getImmune() && !c3->getImmune() ){
-
- std::pair<Cell*, int> c4p = std::make_pair(b, 0);
- c4p = findRemainingBoundary(b, c1, c2, c3);
- int c4c = getCoeff(c4p.first);
- Cell* c4 = c4p.first;
-
- int temp1 = c1c - c1o;
-
- if(c4p.second != 0 && c4c != 0 && !c2->inSubdomain() && !c3->inSubdomain()
- && hasCell(c1) && hasCell(c4) && !hasCell(c2) && !hasCell(c3)) {
-
- c2c = -c2o;
- c3c = -c3o;
- if(temp1 != 0) c2c= c2c*-1;
- if(temp1 != 0) c3c= c3c*-1;
-
- this->removeCell(c1);
- this->removeCell(c4);
- this->addCell(c2, c2c);
- this->addCell(c3, c3c);
- c1->setImmune(false);
- c2->setImmune(true);
- c3->setImmune(false);
- c4->setImmune(false);
+ else if ( (getDim() == 1 && cellsInChain.size() == 1 && cellsNotInChain.size() == 2 && bend) ||
+ (getDim() == 2 && cellsInChain.size() == 2 && cellsNotInChain.size() == 2 && bend)){
+ //printf("bend \n");
+ deform(cellsInChain, cellsNotInChain);
return true;
}
-
- }
-
- else if(c1->getDim() == 1 &&
- b != NULL && c2 != NULL && c3 != NULL && !(*c2 == *c1) && !(*c1 == *c3) && !(*c2 == *c3) &&
- (c1c == 1 || c1c == -1) && !c2->getImmune() && !c3->getImmune() ){
-
- int temp1 = c1c - c1o;
-
- if(!c2->inSubdomain() && !c3->inSubdomain() && hasCell(c1) && !hasCell(c2) && !hasCell(c3)) {
-
- //printf("c2: %d, c3; %d \n", getCoeff(c2), getCoeff(c3));
-
- c2c = -c2o;
- c3c = -c3o;
- if(temp1 != 0) c2c= c2c*-1;
- if(temp1 != 0) c3c= c3c*-1;
-
- this->removeCell(c1);
- this->addCell(c2, c2c);
- this->addCell(c3, c3c);
- c1->setImmune(false);
- c2->setImmune(true);
- c3->setImmune(false);
+ else if ((getDim() == 1 && cellsInChain.size() == 3 && cellsNotInChain.size() == 0) ||
+ (getDim() == 2 && cellsInChain.size() == 4 && cellsNotInChain.size() == 0)){
+ //printf("remove boundary \n");
+ deform(cellsInChain, cellsNotInChain);
return true;
}
-
}
-
return false;
}
void Chain::smoothenChain(){
+ if(!_cellComplex->simplicial()) return;
+
int start = getSize();
double t1 = Cpu();
@@ -809,8 +588,8 @@ void Chain::smoothenChain(){
for(int i = 0; i < 20; i++){
int size = getSize();
for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
- if(!straightenChain(*cit) && getDim() == 2) bendChain(*cit);
- removeBoundary(*cit);
+ if(!deformChain(*cit, true, false) && getDim() == 2) deformChain(*cit, false, true);
+ deformChain(*cit, false, false);
}
deImmuneCells();
eraseNullCells();
@@ -819,9 +598,8 @@ void Chain::smoothenChain(){
if (useless > 5) break;
}
- for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
- straightenChain(*cit);
- }
+ deImmuneCells();
+ for(citer cit = _cells.begin(); cit != _cells.end(); cit++) deformChain(*cit, true, false);
eraseNullCells();
double t2 = Cpu();
@@ -877,15 +655,15 @@ void Chain::createPView(){
Cell* cell = (*cit).first;
int coeff = (*cit).second;
std::vector<MVertex*> v = cell->getVertexVector();
- if(cell->getDim() > 0 && coeff == -1){ // flip orientation
+ if(cell->getDim() > 0 && coeff < 0){ // flip orientation
MVertex* temp = v[0];
v[0] = v[1];
v[1] = temp;
}
MElement *e = factory.create(cell->getTypeForMSH(), v, cell->getNum(), cell->getPartition());
- elements.push_back(e);
+ for(int i = 0; i < abs(coeff); i++) elements.push_back(e);
- std::vector<double> coeffs (1,1);
+ std::vector<double> coeffs (1,abs(coeff));
data[e->getNum()] = coeffs;
}
diff --git a/Geo/ChainComplex.h b/Geo/ChainComplex.h
index 3481a4f19e..1b34844e57 100644
--- a/Geo/ChainComplex.h
+++ b/Geo/ChainComplex.h
@@ -150,13 +150,10 @@ class Chain{
int _dim;
- std::pair<Cell*, int> findRemainingBoundary(Cell* b, Cell* c1, Cell* c2, Cell* c3=NULL);
- Cell* findCommonCbdCell(Cell* c1, Cell* c2, Cell* c3=NULL);
- int findOrientation(Cell* b, Cell* c);
- bool straightenChain( std::pair<Cell*, int> cell);
- bool bendChain( std::pair<Cell*, int> cell);
- bool removeBoundary( std::pair<Cell*, int> cell);
- std::map<Cell*, int, Less_Cell> getBdCellsInChain(Cell* cell);
+
+ bool deform(std::map<Cell*, int, Less_Cell> &cellsInChain, std::map<Cell*, int, Less_Cell> &cellsNotInChain);
+ bool deformChain(std::pair<Cell*, int> cell, bool straighten, bool bend);
+
public:
Chain(){}
--
GitLab