diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index a0fe39fb5f37c85d9353fa7c1abec5847da6af01..3338006ed7ac58cfa35cf3255fa228502a2efd83 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -169,8 +169,7 @@ 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;
- }
- /*
+ }/* FIXME: reduction doesn't work for these
else if(type == MSH_HEX_8 || type == MSH_HEX_27 || type == MSH_HEX_20){
cell = new CHexahedron(vertices, tag, subdomain, boundary);
_simplicial = false;
@@ -215,7 +214,8 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){
delete newCell;
Cell* oldCell = *(insertInfo.first);
if(!subdomain && !boundary){
- int ori = cell->kappa(oldCell);
+ //int ori = cell->kappa(oldCell);
+ int ori = cell->getFacetOri(oldCell);
oldCell->addCoboundaryCell( ori, cell );
oldCell->addOrgCbdCell( ori, cell );
cell->addBoundaryCell( ori, oldCell);
@@ -223,7 +223,8 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){
}
}
else if(!subdomain && !boundary) {
- int ori = cell->kappa(newCell);
+ //int ori = cell->kappa(newCell);
+ int ori = cell->getFacetOri(vertices);
cell->addBoundaryCell( ori, newCell );
cell->addOrgBdCell( ori, newCell );
newCell->addCoboundaryCell( ori, cell);
@@ -256,7 +257,7 @@ int Simplex::kappa(Cell* tau) const{
return value;
}
*/
-
+/*
int Cell::kappa(Cell* tau) const{
for(int i=0; i < tau->getNumVertices(); i++){
if( !(this->hasVertex(tau->getVertex(i)->getNum())) ) return 0;
@@ -306,8 +307,7 @@ int OneSimplex::kappa(Cell* tau) const{
if(tau->getVertex(0) == this->getVertex(0)) return -1;
else return 1;
-}
-
+}*/
void CellComplex::removeCell(Cell* cell, bool other){
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index da55e8a68afdacf0c31cbcb3f66f966a04ccb635..c5a901e6c7c8382ddde06eb8748a43d0ce74f54c 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -99,7 +99,7 @@ class Cell
virtual void setPartition(int num){}
virtual void setImmune(bool immune) { _immune = immune; };
virtual bool getImmune() const { return _immune; };
-
+
// get the number of vertices this cell has
virtual int getNumVertices() const = 0;
virtual MVertex* getVertex(int vertex) const = 0; //{return _vertices.at(vertex);}
@@ -295,7 +295,10 @@ class Cell
virtual int getNumFacets() const { return 0; }
virtual void getFacetVertices(const int num, std::vector<MVertex*> &v) const {};
-
+ // get boundary cell orientation
+ virtual int getFacetOri(Cell* cell) { std::vector<MVertex*> v = cell->getVertexVector(); return getFacetOri(v); }
+ virtual int getFacetOri(std::vector<MVertex*> &v) { return 0; }
+
virtual bool isCombined() { return _combined; }
virtual std::list< std::pair<int, Cell*> > getCells() { std::list< std::pair<int, Cell*> >cells; cells.push_back( std::make_pair(1, this)); return cells; }
virtual int getNumCells() {return 1;}
@@ -326,7 +329,7 @@ class Cell
}
// checks whether lower dimensional simplex tau is on the boundary of this cell
- virtual int kappa(Cell* tau) const;
+ //virtual int kappa(Cell* tau) const;
virtual void writeMSH(FILE *fp, double version=1.0, bool binary=false, int num=0, int elementary=1, int physical=1, int parentNum=0) {}
};
@@ -416,8 +419,14 @@ class OneSimplex : public Simplex, public MLine
v.resize(1);
v[0] = _v[num];
}
-
- int kappa(Cell* tau) const;
+ int getFacetOri(std::vector<MVertex*> &v){
+ if(v.size() != 1) return 0;
+ else if(v[0] == _v[0]) return -1;
+ else if(v[0] == _v[1]) return 1;
+ else return 0;
+ }
+
+ //int kappa(Cell* tau) const;
void printCell() const { printf("Cell dimension: %d, Vertices: %d %d, in subdomain: %d \n", getDim(), _v[0]->getNum(), _v[1]->getNum(), _inSubdomain); }
@@ -466,7 +475,15 @@ class TwoSimplex : public Simplex, public MTriangle
void getFacetVertices(const int num, std::vector<MVertex*> &v) const {
MTriangle::getEdgeVertices(num, v);
}
-
+ int getFacetOri(std::vector<MVertex*> &v){
+ if(v.size() != 2) return 0;
+ MEdge facet = MEdge(v[0], v[1]);
+ int ithFacet = 0;
+ int sign = 0;
+ MTriangle::getEdgeInfo(facet, ithFacet, sign);
+ return sign;
+ }
+
void printCell() const { printf("Cell dimension: %d, Vertices: %d %d %d, in subdomain: %d\n", getDim(), _v[0]->getNum(), _v[1]->getNum(), _v[2]->getNum(), _inSubdomain); }
};
@@ -514,6 +531,14 @@ class CQuadrangle : public Cell, public MQuadrangle
void getFacetVertices(const int num, std::vector<MVertex*> &v) const {
MQuadrangle::getEdgeVertices(num, v);
}
+ int getFacetOri(std::vector<MVertex*> &v){
+ if(v.size() != 2) return 0;
+ MEdge facet = MEdge(v[0], v[1]);
+ int ithFacet = 0;
+ int sign = 0;
+ MQuadrangle::getEdgeInfo(facet, ithFacet, sign);
+ return sign;
+ }
void printCell() const { printf("Cell dimension: %d, Vertices: %d %d %d %d, in subdomain: %d\n",
getDim(), _v[0]->getNum(), _v[1]->getNum(), _v[2]->getNum(),
@@ -562,7 +587,16 @@ class ThreeSimplex : public Simplex, public MTetrahedron
void getFacetVertices(const int num, std::vector<MVertex*> &v) const {
MTetrahedron::getFaceVertices(num, v);
}
-
+ int getFacetOri(std::vector<MVertex*> &v){
+ if(v.size() != 3) return 0;
+ MFace facet = MFace(v);
+ int ithFacet = 0;
+ int sign = 0;
+ int rot = 0;
+ MTetrahedron::getFaceInfo(facet, ithFacet, sign, rot);
+ return sign;
+ }
+
virtual void printCell() const { printf("Cell dimension: %d, Vertices: %d %d %d %d, in subdomain: %d \n", getDim(), _v[0]->getNum(), _v[1]->getNum(), _v[2]->getNum(), _v[3]->getNum(), _inSubdomain); }
};
@@ -612,6 +646,15 @@ class CHexahedron : public Cell, public MHexahedron
void getFacetVertices(const int num, std::vector<MVertex*> &v) const {
MHexahedron::getFaceVertices(num, v);
}
+ int getFacetOri(std::vector<MVertex*> &v){
+ if(v.size() != 4) return 0;
+ MFace facet = MFace(v);
+ int ithFacet = 0;
+ int sign = 0;
+ int rot = 0;
+ MHexahedron::getFaceInfo(facet, ithFacet, sign, rot);
+ return sign;
+ }
virtual void printCell() const { printf("Cell dimension: %d, Vertices: %d %d %d %d %d %d %d %d, in subdomain: %d \n", getDim(), _v[0]->getNum(), _v[1]->getNum(), _v[2]->getNum(), _v[3]->getNum(), _v[4]->getNum(), _v[5]->getNum(), _v[6]->getNum(), _v[7]->getNum(), _inSubdomain); }
};
@@ -661,6 +704,15 @@ class CPrism : public Cell, public MPrism
void getFacetVertices(const int num, std::vector<MVertex*> &v) const {
MPrism::getFaceVertices(num, v);
}
+ int getFacetOri(std::vector<MVertex*> &v){
+ if(v.size() != 4 && v.size() != 3) return 0;
+ MFace facet = MFace(v);
+ int ithFacet = 0;
+ int sign = 0;
+ int rot = 0;
+ MPrism::getFaceInfo(facet, ithFacet, sign, rot);
+ return sign;
+ }
virtual void printCell() const { printf("Cell dimension: %d, Vertices: %d %d %d %d %d %d, in subdomain: %d \n", getDim(), _v[0]->getNum(), _v[1]->getNum(), _v[2]->getNum(), _v[3]->getNum(), _v[4]->getNum(), _v[5]->getNum(), _inSubdomain); }
};
@@ -709,6 +761,15 @@ class CPyramid : public Cell, public MPyramid
void getFacetVertices(const int num, std::vector<MVertex*> &v) const {
MPyramid::getFaceVertices(num, v);
}
+ int getFacetOri(std::vector<MVertex*> &v){
+ if(v.size() != 4 && v.size() != 3) return 0;
+ MFace facet = MFace(v);
+ int ithFacet = 0;
+ int sign = 0;
+ int rot = 0;
+ MPyramid::getFaceInfo(facet, ithFacet, sign, rot);
+ return sign;
+ }
virtual void printCell() const { printf("Cell dimension: %d, Vertices: %d %d %d %d %d, in subdomain: %d \n", getDim(), _v[0]->getNum(), _v[1]->getNum(), _v[2]->getNum(), _v[3]->getNum(), _v[4]->getNum(), _inSubdomain); }
};
@@ -841,7 +902,7 @@ class CombinedCell : public Cell{
int getSortedVertex(int vertex) const { return _vs.at(vertex); }
std::vector<MVertex*> getVertexVector() const { return _v; }
- int kappa(Cell* tau) const { return 0; }
+ //int kappa(Cell* tau) const { return 0; }
// true if this cell has given vertex
bool hasVertex(int vertex) const {
@@ -1030,7 +1091,7 @@ class CellComplex
// kappa for two cells of this cell complex
// implementation will vary depending on cell type
- inline int kappa(Cell* sigma, Cell* tau) const { return sigma->kappa(tau); }
+ //inline int kappa(Cell* sigma, Cell* tau) const { return sigma->kappa(tau); }
// check whether two cells both belong to subdomain or if neither one does
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index 53d0dce5e9e8f5552505669f1b02442270dc9097..42ee1ee0d18d8afc2e958b70e43a18ce03fe119f 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -471,7 +471,6 @@ Chain::Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs, CellComp
int coeff = (*it).first;
_cells.insert( std::make_pair(subCell, coeffs.at(i)*coeff));
}
- //_cells.push_back( std::make_pair(cell, coeffs.at(i)) );
}
i++;
}
@@ -497,27 +496,23 @@ bool Chain::deform(std::map<Cell*, int, Less_Cell> &cellsInChain, std::map<Cell*
cc.push_back(getCoeff(c));
bc.push_back((*cit).second);
}
+ removeCell(c);
}
-
- // FIXME: orientations don't get right on 2-chains with bending or subdomain cells (disabled for now).
-
- bool flip = false;
- if(cc[0] != bc[0]) flip = true;
+ if(cc.empty()) return false;
+ int inout = cc[0]*bc[0];
for(int i = 0; i < cc.size(); i++){
//printf("cc: %d, bc: %d \n", cc.at(i), bc.at(i));
- if(flip && cc[i] == bc[i]) return false;
- if(!flip && cc[i] != bc[i]) return false;
+ if(cc[i]*bc[i] != inout) printf("Error: Chain smoothening orientation mismatch! \n");
}
- for(citer cit = cellsInChain.begin(); cit != cellsInChain.end(); cit++) removeCell((*cit).first);
+ //for(citer cit = cellsInChain.begin(); cit != cellsInChain.end(); cit++) removeCell((*cit).first);
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;
- if(flip) coeff = coeff*-1;
+ int coeff = -1*inout*(*cit).second;
addCell(c, coeff);
//c->printCell();
n++;
@@ -527,17 +522,17 @@ bool Chain::deform(std::map<Cell*, int, Less_Cell> &cellsInChain, std::map<Cell*
return true;
}
-bool Chain::deformChain(std::pair<Cell*, int> cell, bool straighten, bool bend){
+bool Chain::deformChain(std::pair<Cell*, int> cell, bool bend){
Cell* c1 = cell.first;
-
std::map<Cell*, int, Less_Cell> c1Cbd = c1->getOrgCbd();
for(citer cit = c1Cbd.begin(); cit != c1Cbd.end(); cit++){
std::map<Cell*, int, Less_Cell> cellsInChain;
std::map<Cell*, int, Less_Cell> cellsNotInChain;
-
Cell* c1CbdCell = (*cit).first;
+ //c1CbdCell->printCell();
+ //c1CbdCell->printOrgBd();
std::map<Cell*, int, Less_Cell> c1CbdBd = c1CbdCell->getOrgBd();
for(citer cit2 = c1CbdBd.begin(); cit2 != c1CbdBd.end(); cit2++){
Cell* c1CbdBdCell = (*cit2).first;
@@ -552,7 +547,7 @@ bool Chain::deformChain(std::pair<Cell*, int> cell, bool straighten, bool bend){
Cell* c = (*cit2).first;
int coeff = getCoeff(c);
if(c->getImmune()) next = true;
- if(/* FIXME: bend && */c->inSubdomain()) next = true;
+ if(c->inSubdomain()) bend = false;
if(coeff > 1 || coeff < -1) next = true;
}
@@ -563,18 +558,19 @@ bool Chain::deformChain(std::pair<Cell*, int> cell, bool straighten, bool bend){
if(next) continue;
+ if(getDim() != 2) bend = false;
+
//printf("dim: %d, in chain: %d, not in chain: %d \n", getDim(), cellsInChain.size(), cellsNotInChain.size());
- if( (getDim() == 1 && cellsInChain.size() == 2 && cellsNotInChain.size() == 1 && straighten) ||
- (getDim() == 2 && cellsInChain.size() == 3 && cellsNotInChain.size() == 1 && straighten)){
+ if( (getDim() == 1 && cellsInChain.size() == 2 && cellsNotInChain.size() == 1) ||
+ (getDim() == 2 && cellsInChain.size() == 3 && cellsNotInChain.size() == 1)){
//printf("straighten \n");
return deform(cellsInChain, cellsNotInChain);
}
else if ( (getDim() == 1 && cellsInChain.size() == 1 && cellsNotInChain.size() == 2 && bend) ||
(getDim() == 2 && cellsInChain.size() == 2 && cellsNotInChain.size() == 2 && bend)){
//printf("bend \n");
- //FIXME: return deform(cellsInChain, cellsNotInChain);
- return false;
+ return deform(cellsInChain, cellsNotInChain);
}
else if ((getDim() == 1 && cellsInChain.size() == 3 && cellsNotInChain.size() == 0) ||
(getDim() == 2 && cellsInChain.size() == 4 && cellsNotInChain.size() == 0)){
@@ -587,6 +583,7 @@ bool Chain::deformChain(std::pair<Cell*, int> cell, bool straighten, bool bend){
}
void Chain::smoothenChain(){
+
if(!_cellComplex->simplicial()) return;
int start = getSize();
@@ -596,8 +593,8 @@ void Chain::smoothenChain(){
for(int i = 0; i < 20; i++){
int size = getSize();
for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
- if(!deformChain(*cit, true, false) && getDim() == 2) deformChain(*cit, false, true);
- deformChain(*cit, false, false);
+ deformChain(*cit, true);
+ deformChain(*cit, false);
}
deImmuneCells();
eraseNullCells();
@@ -607,7 +604,7 @@ void Chain::smoothenChain(){
}
deImmuneCells();
- for(citer cit = _cells.begin(); cit != _cells.end(); cit++) deformChain(*cit, true, false);
+ for(citer cit = _cells.begin(); cit != _cells.end(); cit++) deformChain(*cit, false);
eraseNullCells();
double t2 = Cpu();
diff --git a/Geo/ChainComplex.h b/Geo/ChainComplex.h
index 1b34844e577d62f2f60feaf900f918c5ba38a17a..3a474b3a2ca0a4463419f37a143030031b819980 100644
--- a/Geo/ChainComplex.h
+++ b/Geo/ChainComplex.h
@@ -152,7 +152,7 @@ class Chain{
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);
+ bool deformChain(std::pair<Cell*, int> cell, bool bend);
public:
diff --git a/Geo/MHexahedron.cpp b/Geo/MHexahedron.cpp
index 8c53b9438c6ceac03466a2461b54c74f52b9577c..88adfc2c08038ed596dd4b2fc8757e841a3556a7 100644
--- a/Geo/MHexahedron.cpp
+++ b/Geo/MHexahedron.cpp
@@ -26,3 +26,38 @@ void MHexahedron::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
*npts = getNGQHPts(pOrder);
*pts = getGQHPts(pOrder);
}
+
+void MHexahedron::getFaceInfo (const MFace & face, int &ithFace, int &sign, int &rot)const{
+ for (ithFace=0;ithFace<6;ithFace++){
+ MVertex *v0 = _v[faces_hexa(ithFace,0)];
+ MVertex *v1 = _v[faces_hexa(ithFace,1)];
+ MVertex *v2 = _v[faces_hexa(ithFace,2)];
+ MVertex *v3 = _v[faces_hexa(ithFace,3)];
+
+ if (v0 == face.getVertex(0) && v1 == face.getVertex(1) && v2 == face.getVertex(2) && v3 == face.getVertex(3)){
+ sign = 1; rot = 0; return;
+ }
+ if (v0 == face.getVertex(1) && v1 == face.getVertex(2) && v3 == face.getVertex(3) && v2 == face.getVertex(0)){
+ sign = 1; rot = 1; return;
+ }
+ if (v0 == face.getVertex(2) && v3 == face.getVertex(3) && v1 == face.getVertex(0) && v2 == face.getVertex(1)){
+ sign = 1; rot = 2; return;
+ }
+ if (v0 == face.getVertex(3) && v3 == face.getVertex(0) && v1 == face.getVertex(1) && v2 == face.getVertex(2)){
+ sign = 1; rot = 3; return;
+ }
+ if (v0 == face.getVertex(0) && v1 == face.getVertex(3) && v2 == face.getVertex(2) && v3 == face.getVertex(1)){
+ sign = -1; rot = 0; return;
+ }
+ if (v0 == face.getVertex(3) && v1 == face.getVertex(2) && v2 == face.getVertex(1) && v3 == face.getVertex(0)){
+ sign = -1; rot = 1; return;
+ }
+ if (v0 == face.getVertex(2) && v1 == face.getVertex(1) && v2 == face.getVertex(0) && v3 == face.getVertex(3)){
+ sign = -1; rot = 2; return;
+ }
+ if (v0 == face.getVertex(1) && v1 == face.getVertex(0) && v2 == face.getVertex(3) && v3 == face.getVertex(2)){
+ sign = -1; rot = 3; return;
+ }
+ }
+ throw;
+}
diff --git a/Geo/MHexahedron.h b/Geo/MHexahedron.h
index 350a23e23bcea6aebbbca4afc10174c8552f5d07..6be76e7c8e01b9931aa552b05d1fdda79b563d27 100644
--- a/Geo/MHexahedron.h
+++ b/Geo/MHexahedron.h
@@ -92,6 +92,7 @@ class MHexahedron : public MElement {
_v[faces_hexa(num, 2)],
_v[faces_hexa(num, 3)]);
}
+ virtual void getFaceInfo (const MFace & face, int &ithFace, int &sign, int &rot)const;
virtual int getNumFacesRep(){ return 12; }
virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
{