diff --git a/Geo/Cell.cpp b/Geo/Cell.cpp
index ea9ccbc8aeef8abf1eebeef578ab86bb866bf8dd..ae9edf9f988dcd1723abbf1e428efc176c8e1ecd 100755
--- a/Geo/Cell.cpp
+++ b/Geo/Cell.cpp
@@ -6,6 +6,12 @@
// Contributed by Matti Pellikka <matti.pellikka@tut.fi>.
#include "Cell.h"
+#include "MTriangle.h"
+#include "MQuadrangle.h"
+#include "MTetrahedron.h"
+#include "MPyramid.h"
+#include "MHexahedron.h"
+#include "MPrism.h"
bool Less_Cell::operator()(const Cell* c1, const Cell* c2) const
{
@@ -23,35 +29,26 @@ bool Less_Cell::operator()(const Cell* c1, const Cell* c2) const
return false;
}
-int Cell::_globalNum = 0;
-Cell::Cell(int dim, std::vector<MVertex*>& v)
+bool Less_VertexNumIndex::operator()(const std::pair<MVertex*, int> v1,
+ const std::pair<MVertex*, int> v2) const
{
- _dim = dim;
- _domain = 0;
- _combined = false;
- _immune = false;
- _v = v;
- std::sort(_v.begin(), _v.end(), MVertexLessThanNum());
-
- _num = 0;
- _index = 0;
+ return (v1.first->getNum() < v2.first->getNum());
}
+int Cell::_globalNum = 0;
+
Cell::Cell(MElement* element, int domain)
{
_dim = element->getDim();
_domain = domain;
_combined = false;
_immune = false;
+ _num = 0;
- for(int i = 0; i < element->getNumPrimaryVertices(); i++){
+ for(int i = 0; i < element->getNumPrimaryVertices(); i++)
_v.push_back(element->getVertex(i));
- }
- std::sort(_v.begin(), _v.end(), MVertexLessThanNum());
-
- _num = 0;
- _index = 0;
+ _sortVertexIndices();
}
Cell::Cell(Cell* parent, int i)
@@ -60,144 +57,240 @@ Cell::Cell(Cell* parent, int i)
_domain = parent->getDomain();
_combined = false;
_immune = false;
+ _num = 0;
parent->findBdElement(i, _v);
- std::sort(_v.begin(), _v.end(), MVertexLessThanNum());
+ _sortVertexIndices();
+}
- _num = 0;
- _index = 0;
+void Cell::_sortVertexIndices()
+{
+ std::vector< std::pair<MVertex*, int> > si;
+
+ for(unsigned int i = 0; i < _v.size(); i++)
+ si.push_back( std::make_pair(_v[i], i) );
+ std::sort(si.begin(), si.end(), Less_VertexNumIndex());
+
+ for(unsigned int i = 0; i < si.size(); i++)
+ _si.push_back(si[i].second);
+}
+
+inline int Cell::getSortedVertex(int vertex) const
+{
+ return _v[(int)_si[vertex]]->getNum();
}
void Cell::findBdElement(int i, std::vector<MVertex*>& vertices) const
{
vertices.clear();
- if(_dim == 1){ // i = 0 -> ori = -1, i = 1 -> ori = 1
+ switch (_dim) {
+ case 1:
vertices.push_back(_v[i]);
- }
- else if(_dim == 2){
- if(i == 0){ // ori = 1
- vertices.push_back(_v[0]);
- vertices.push_back(_v[1]);
- }
- else if(i == 1){ // ori = 1
- vertices.push_back(_v[1]);
- vertices.push_back(_v[2]);
- }
- else if(i == 2){ // ori = 1
- vertices.push_back(_v[2]);
- vertices.push_back(_v[0]);
- }
- }
- else if(_dim == 3){
- if(i == 0){
- vertices.push_back(_v[0]);
- vertices.push_back(_v[2]);
- vertices.push_back(_v[1]);
- }
- else if(i == 1){
- vertices.push_back(_v[0]);
- vertices.push_back(_v[1]);
- vertices.push_back(_v[3]);
- }
- else if(i == 2){
- vertices.push_back(_v[0]);
- vertices.push_back(_v[3]);
- vertices.push_back(_v[2]);
+ return;
+ case 2:
+ switch (getNumVertices()) {
+ case 3:
+ for(int j = 0; j < 2; j++)
+ vertices.push_back(_v[MTriangle::edges_tri(i, j)]);
+ return;
+ case 4:
+ for(int j = 0; j < 2; j++)
+ vertices.push_back(_v[MQuadrangle::edges_quad(i, j)]);
+ return;
+ default: return;
}
- else if(i == 3){
- vertices.push_back(_v[3]);
- vertices.push_back(_v[1]);
- vertices.push_back(_v[2]);
+ case 3:
+ switch (getNumVertices()) {
+ case 4:
+ for(int j = 0; j < 3; j++)
+ vertices.push_back(_v[MTetrahedron::faces_tetra(i, j)]);
+ return;
+ case 5:
+ if(i < 3) {
+ for(int j = 0; j < 3; j++)
+ vertices.push_back(_v[MPyramid::faces_pyramid(i, j)]);
+ }
+ else {
+ vertices.push_back(_v[0]);
+ vertices.push_back(_v[3]);
+ vertices.push_back(_v[2]);
+ vertices.push_back(_v[1]);
+ }
+ return;
+ case 6:
+ if(i < 2)
+ for(int j = 0; j < 3; j++)
+ vertices.push_back(_v[MPrism::faces_prism(i, j)]);
+ else
+ for(int j = 0; j < 4; j++)
+ vertices.push_back(_v[MPrism::faces_prism(i, j)]);
+ return;
+ case 8:
+ for(int j = 0; j < 4; j++)
+ vertices.push_back(_v[MHexahedron::faces_hexa(i, j)]);
+ return;
+ default: return;
}
+ default: return;
}
}
int Cell::getNumBdElements() const
{
- if(_dim == 0) return 0;
- else if(_dim == 1) return 2;
- else if(_dim == 2) return 3;
- else if(_dim == 3) return 4;
- else return 0;
+ switch (_dim) {
+ case 0: return 0;
+ case 1: return 2;
+ case 2:
+ switch (getNumVertices()) {
+ case 3: return 3;
+ case 4: return 4;
+ default: return 0;
+ }
+ case 3:
+ switch (getNumVertices()) {
+ case 4: return 4;
+ case 5: return 5;
+ case 6: return 5;
+ case 8: return 6;
+ default: return 0;
+ }
+ default: return 0;
+ }
}
-int Cell::findBdCellOrientation(Cell* cell) const
+
+int Cell::findBdCellOrientation(Cell* cell, int i) const
{
std::vector<MVertex*> v;
cell->getMeshVertices(v);
- if(_dim == 0) return 0;
- else if(_dim == 1){
+ switch (_dim) {
+ case 0: return 0;
+ case 1:
if(v[0]->getNum() == _v[0]->getNum()) return -1;
else if(v[0]->getNum() == _v[1]->getNum()) return 1;
+ return 0;
+ case 2:
+ switch (getNumVertices()) {
+ case 3:
+ if(_v[MTriangle::edges_tri(i, 0)]->getNum() == v[0]->getNum() &&
+ _v[MTriangle::edges_tri(i, 1)]->getNum() == v[1]->getNum())
+ return 1;
+ if(_v[MTriangle::edges_tri(i, 1)]->getNum() == v[0]->getNum() &&
+ _v[MTriangle::edges_tri(i, 0)]->getNum() == v[1]->getNum())
+ return -1;
+ case 4:
+ if(_v[MQuadrangle::edges_quad(i, 0)]->getNum() == v[0]->getNum() &&
+ _v[MQuadrangle::edges_quad(i, 1)]->getNum() == v[1]->getNum())
+ return 1;
+ if(_v[MQuadrangle::edges_quad(i, 0)]->getNum() == v[1]->getNum() &&
+ _v[MQuadrangle::edges_quad(i, 1)]->getNum() == v[0]->getNum())
+ return -1;
+ default: return 0;
+ }
+ case 3:
+ switch (getNumVertices()) {
+ case 4:
+ if(_v[MTetrahedron::faces_tetra(i, 0)]->getNum() == v[0]->getNum() &&
+ _v[MTetrahedron::faces_tetra(i, 1)]->getNum() == v[1]->getNum() &&
+ _v[MTetrahedron::faces_tetra(i, 2)]->getNum() == v[2]->getNum())
+ return 1;
+ if(_v[MTetrahedron::faces_tetra(i, 0)]->getNum() == v[1]->getNum() &&
+ _v[MTetrahedron::faces_tetra(i, 1)]->getNum() == v[2]->getNum() &&
+ _v[MTetrahedron::faces_tetra(i, 2)]->getNum() == v[0]->getNum())
+ return 1;
+ if(_v[MTetrahedron::faces_tetra(i, 0)]->getNum() == v[2]->getNum() &&
+ _v[MTetrahedron::faces_tetra(i, 1)]->getNum() == v[0]->getNum() &&
+ _v[MTetrahedron::faces_tetra(i, 2)]->getNum() == v[1]->getNum())
+ return 1;
+ if(_v[MTetrahedron::faces_tetra(i, 0)]->getNum() == v[0]->getNum() &&
+ _v[MTetrahedron::faces_tetra(i, 1)]->getNum() == v[2]->getNum() &&
+ _v[MTetrahedron::faces_tetra(i, 2)]->getNum() == v[1]->getNum())
+ return -1;
+ if(_v[MTetrahedron::faces_tetra(i, 0)]->getNum() == v[1]->getNum() &&
+ _v[MTetrahedron::faces_tetra(i, 1)]->getNum() == v[0]->getNum() &&
+ _v[MTetrahedron::faces_tetra(i, 2)]->getNum() == v[2]->getNum())
+ return -1;
+ if(_v[MTetrahedron::faces_tetra(i, 0)]->getNum() == v[2]->getNum() &&
+ _v[MTetrahedron::faces_tetra(i, 1)]->getNum() == v[1]->getNum() &&
+ _v[MTetrahedron::faces_tetra(i, 2)]->getNum() == v[0]->getNum())
+ return -1;
+ case 5: return 0;
+ case 6: return 0;
+ case 8:
+ if(_v[MHexahedron::faces_hexa(i, 0)]->getNum() == v[0]->getNum() &&
+ _v[MHexahedron::faces_hexa(i, 1)]->getNum() == v[1]->getNum() &&
+ _v[MHexahedron::faces_hexa(i, 2)]->getNum() == v[2]->getNum() &&
+ _v[MHexahedron::faces_hexa(i, 3)]->getNum() == v[3]->getNum())
+ return 1;
+ if(_v[MHexahedron::faces_hexa(i, 0)]->getNum() == v[1]->getNum() &&
+ _v[MHexahedron::faces_hexa(i, 1)]->getNum() == v[2]->getNum() &&
+ _v[MHexahedron::faces_hexa(i, 2)]->getNum() == v[3]->getNum() &&
+ _v[MHexahedron::faces_hexa(i, 3)]->getNum() == v[0]->getNum())
+ return 1;
+ if(_v[MHexahedron::faces_hexa(i, 0)]->getNum() == v[2]->getNum() &&
+ _v[MHexahedron::faces_hexa(i, 1)]->getNum() == v[3]->getNum() &&
+ _v[MHexahedron::faces_hexa(i, 2)]->getNum() == v[0]->getNum() &&
+ _v[MHexahedron::faces_hexa(i, 3)]->getNum() == v[1]->getNum())
+ return 1;
+ if(_v[MHexahedron::faces_hexa(i, 0)]->getNum() == v[3]->getNum() &&
+ _v[MHexahedron::faces_hexa(i, 1)]->getNum() == v[0]->getNum() &&
+ _v[MHexahedron::faces_hexa(i, 2)]->getNum() == v[1]->getNum() &&
+ _v[MHexahedron::faces_hexa(i, 3)]->getNum() == v[2]->getNum())
+ return 1;
+ if(_v[MHexahedron::faces_hexa(i, 0)]->getNum() == v[0]->getNum() &&
+ _v[MHexahedron::faces_hexa(i, 1)]->getNum() == v[3]->getNum() &&
+ _v[MHexahedron::faces_hexa(i, 2)]->getNum() == v[2]->getNum() &&
+ _v[MHexahedron::faces_hexa(i, 3)]->getNum() == v[1]->getNum())
+ return -1;
+ if(_v[MHexahedron::faces_hexa(i, 0)]->getNum() == v[3]->getNum() &&
+ _v[MHexahedron::faces_hexa(i, 1)]->getNum() == v[2]->getNum() &&
+ _v[MHexahedron::faces_hexa(i, 2)]->getNum() == v[1]->getNum() &&
+ _v[MHexahedron::faces_hexa(i, 3)]->getNum() == v[0]->getNum())
+ return -1;
+ if(_v[MHexahedron::faces_hexa(i, 0)]->getNum() == v[2]->getNum() &&
+ _v[MHexahedron::faces_hexa(i, 1)]->getNum() == v[1]->getNum() &&
+ _v[MHexahedron::faces_hexa(i, 2)]->getNum() == v[0]->getNum() &&
+ _v[MHexahedron::faces_hexa(i, 3)]->getNum() == v[3]->getNum())
+ return -1;
+ if(_v[MHexahedron::faces_hexa(i, 0)]->getNum() == v[1]->getNum() &&
+ _v[MHexahedron::faces_hexa(i, 1)]->getNum() == v[0]->getNum() &&
+ _v[MHexahedron::faces_hexa(i, 2)]->getNum() == v[3]->getNum() &&
+ _v[MHexahedron::faces_hexa(i, 3)]->getNum() == v[2]->getNum())
+ return -1;
+ default: return 0;
+ }
+ default: return 0;
}
- if(_dim == 2){
- if (v[0]->getNum() == _v[0]->getNum() && v[1]->getNum() == _v[1]->getNum() ) return 1;
- else if(v[0]->getNum() == _v[1]->getNum() && v[1]->getNum() == _v[0]->getNum() ) return -1;
-
- else if(v[0]->getNum() == _v[1]->getNum() && v[1]->getNum() == _v[2]->getNum() ) return 1;
- else if(v[0]->getNum() == _v[2]->getNum() && v[1]->getNum() == _v[1]->getNum() ) return -1;
-
- else if(v[0]->getNum() == _v[2]->getNum() && v[1]->getNum() == _v[0]->getNum() ) return 1;
- else if(v[0]->getNum() == _v[0]->getNum() && v[1]->getNum() == _v[2]->getNum() ) return -1;
- }
- else if(_dim == 3) {
- if (v[0]->getNum() == _v[0]->getNum() && v[1]->getNum() == _v[2]->getNum() && v[2]->getNum() == _v[1]->getNum() ) return 1;
- else if(v[0]->getNum() == _v[0]->getNum() && v[1]->getNum() == _v[1]->getNum() && v[2]->getNum() == _v[2]->getNum()) return -1;
-
- else if(v[0]->getNum() == _v[2]->getNum() && v[1]->getNum() == _v[1]->getNum() && v[2]->getNum() == _v[0]->getNum()) return 1;
- else if(v[0]->getNum() == _v[2]->getNum() && v[1]->getNum() == _v[0]->getNum() && v[2]->getNum() == _v[1]->getNum()) return -1;
-
- else if(v[0]->getNum() == _v[1]->getNum() && v[1]->getNum() == _v[0]->getNum() && v[2]->getNum() == _v[2]->getNum()) return 1;
- else if(v[0]->getNum() == _v[1]->getNum() && v[1]->getNum() == _v[2]->getNum() && v[2]->getNum() == _v[0]->getNum()) return -1;
-
-
- else if(v[0]->getNum() == _v[0]->getNum() && v[1]->getNum() == _v[1]->getNum() && v[2]->getNum() == _v[3]->getNum() ) return 1;
- else if(v[0]->getNum() == _v[0]->getNum() && v[1]->getNum() == _v[3]->getNum() && v[2]->getNum() == _v[1]->getNum()) return -1;
-
- else if(v[0]->getNum() == _v[1]->getNum() && v[1]->getNum() == _v[3]->getNum() && v[2]->getNum() == _v[0]->getNum()) return 1;
- else if(v[0]->getNum() == _v[1]->getNum() && v[1]->getNum() == _v[0]->getNum() && v[2]->getNum() == _v[3]->getNum()) return -1;
-
- else if(v[0]->getNum() == _v[3]->getNum() && v[1]->getNum() == _v[0]->getNum() && v[2]->getNum() == _v[1]->getNum()) return 1;
- else if(v[0]->getNum() == _v[3]->getNum() && v[1]->getNum() == _v[1]->getNum() && v[2]->getNum() == _v[0]->getNum()) return -1;
-
-
- else if(v[0]->getNum() == _v[0]->getNum() && v[1]->getNum() == _v[3]->getNum() && v[2]->getNum() == _v[2]->getNum() ) return 1;
- else if(v[0]->getNum() == _v[0]->getNum() && v[1]->getNum() == _v[2]->getNum() && v[2]->getNum() == _v[3]->getNum()) return -1;
-
- else if(v[0]->getNum() == _v[3]->getNum() && v[1]->getNum() == _v[2]->getNum() && v[2]->getNum() == _v[3]->getNum()) return 1;
- else if(v[0]->getNum() == _v[3]->getNum() && v[1]->getNum() == _v[3]->getNum() && v[2]->getNum() == _v[2]->getNum()) return -1;
-
- else if(v[0]->getNum() == _v[2]->getNum() && v[1]->getNum() == _v[0]->getNum() && v[2]->getNum() == _v[3]->getNum()) return 1;
- else if(v[0]->getNum() == _v[2]->getNum() && v[1]->getNum() == _v[3]->getNum() && v[2]->getNum() == _v[0]->getNum()) return -1;
-
-
- else if(v[0]->getNum() == _v[3]->getNum() && v[1]->getNum() == _v[1]->getNum() && v[2]->getNum() == _v[2]->getNum() ) return 1;
- else if(v[0]->getNum() == _v[3]->getNum() && v[1]->getNum() == _v[2]->getNum() && v[2]->getNum() == _v[1]->getNum()) return -1;
-
- else if(v[0]->getNum() == _v[1]->getNum() && v[1]->getNum() == _v[2]->getNum() && v[2]->getNum() == _v[3]->getNum()) return 1;
- else if(v[0]->getNum() == _v[1]->getNum() && v[1]->getNum() == _v[3]->getNum() && v[2]->getNum() == _v[2]->getNum()) return -1;
-
- else if(v[0]->getNum() == _v[2]->getNum() && v[1]->getNum() == _v[3]->getNum() && v[2]->getNum() == _v[1]->getNum()) return 1;
- else if(v[0]->getNum() == _v[2]->getNum() && v[1]->getNum() == _v[1]->getNum() && v[2]->getNum() == _v[3]->getNum()) return -1;
- }
- return 0;
}
int Cell::getTypeMSH() const
{
- if(_dim == 0) return 15;
- else if(_dim == 1) return 1;
- else if(_dim == 2) return 2;
- else if(_dim == 3) return 4;
- else return 0;
+ switch (_dim) {
+ case 0: return MSH_PNT;
+ case 1: return MSH_LIN_2;
+ case 2:
+ switch (getNumVertices()) {
+ case 3: return MSH_TRI_3;
+ case 4: return MSH_QUA_4;
+ default: return 0;
+ }
+ case 3:
+ switch (getNumVertices()) {
+ case 4: return MSH_TET_4;
+ case 5: return MSH_PYR_5;
+ case 6: return MSH_PRI_6;
+ case 8: return MSH_HEX_8;
+ default: return 0;
+ }
+ default: return 0;
+ }
}
bool Cell::hasVertex(int vertex) const
{
- //_vs
std::vector<int> v;
for(unsigned int i = 0; i < _v.size(); i++) {
- v.push_back(_v[i]->getNum());
+ v.push_back(_v[(int)_si[i]]->getNum());
}
std::vector<int>::const_iterator it = std::find(v.begin(), v.end(),
vertex);
@@ -219,7 +312,7 @@ void Cell::printCell()
printf("%d-cell %d: \n" , getDim(), getNum());
printf("Vertices: ");
for(int i = 0; i < this->getNumVertices(); i++){
- printf("%d ", this->getVertex(i));
+ printf("%d ", this->getMeshVertex(i)->getNum());
}
printf(", in subdomain: %d, ", inSubdomain());
printf("combined: %d. \n" , isCombined() );
@@ -239,7 +332,6 @@ void Cell::restoreCell(){
}
for(unsigned int i = 0; i < toRemove.size(); i++) _bd.erase(toRemove[i]);
_combined = false;
- _index = 0;
_immune = false;
}
@@ -351,7 +443,7 @@ void Cell::printCoboundary(bool orig)
}
}*/
-CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell()
+CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co)
{
// use "smaller" cell as c2
if(c1->getNumCells() < c2->getNumCells()){
@@ -361,7 +453,6 @@ CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell()
}
_num = ++_globalNum;
- _index = c1->getIndex();
_domain = c1->getDomain();
_combined = true;
_immune = false;
@@ -417,12 +508,12 @@ CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell()
}
-CombinedCell::CombinedCell(std::vector<Cell*>& cells) : Cell()
+CombinedCell::CombinedCell(std::vector<Cell*>& cells)
{
_num = ++_globalNum;
- _index = cells.at(0)->getIndex();
_domain = cells.at(0)->getDomain();
_combined = true;
+ _immune = false;
// cells
for(unsigned int i = 0; i < cells.size(); i++){
diff --git a/Geo/Cell.h b/Geo/Cell.h
index 7fb65e077b3ef2d99d6eb9f73b6877f352fe556d..da83de9f6ccd25de3868be2250b705c48ae03fec 100644
--- a/Geo/Cell.h
+++ b/Geo/Cell.h
@@ -19,11 +19,17 @@ public:
bool operator()(const Cell* c1, const Cell* c2) const;
};
+class Less_VertexNumIndex {
+public:
+ bool operator()(const std::pair<MVertex*, int> v1,
+ const std::pair<MVertex*, int> v2) const;
+};
+
// Class to save cell boundary orientation information
class BdInfo {
private:
- short int _ori;
- short int _origOri;
+ signed char _ori;
+ signed char _origOri;
public:
BdInfo(int ori) { _ori = ori; _origOri = 0; }
@@ -41,15 +47,13 @@ class BdInfo {
class Cell {
protected:
+
static int _globalNum;
int _num;
- // mutable index for each cell (used to create boundary operator matrices)
- int _index;
-
- int _domain;
+ char _domain;
- // whether this cell a combinded cell of elemetary cells
+ // whether this cell a combinded cell of elementary cells
bool _combined;
// for some algorithms to omit this cell
bool _immune;
@@ -60,15 +64,16 @@ class Cell {
private:
- int _dim;
- // sorted vertices of this cell (used for ordering of the cells)
+ char _dim;
std::vector<MVertex*> _v;
+ // sorted vertices of this cell (used for ordering of the cells)
+ std::vector<char> _si;
- public:
+ inline void _sortVertexIndices();
- Cell(int dim, std::vector<MVertex*>& v);
- Cell() : _num(0), _dim(0), _index(0), _domain(0), _combined(false), _immune(false) {}
+ public:
+ Cell() {}
Cell(MElement* element, int domain);
Cell(Cell* parent, int i);
@@ -77,26 +82,21 @@ class Cell {
void setNum(int num) { _num = num; };
int getTypeMSH() const;
virtual int getDim() const { return _dim; }
- bool inSubdomain() const { if(_domain != 0) return true; else return false; }
+ bool inSubdomain() const { return _domain ? true : false; }
void getMeshVertices(std::vector<MVertex*>& v) const { v = _v; }
- int getIndex() const { return _index; };
- void setIndex(int index) { _index = index; };
void setImmune(bool immune) { _immune = immune; };
bool getImmune() const { return _immune; };
int getNumSortedVertices() const { return _v.size(); }
- int getSortedVertex(int vertex) const { return _v.at(vertex)->getNum(); }
+ inline int getSortedVertex(int vertex) const;
int getNumVertices() const { return _v.size(); }
- int getVertex(int vertex) const { return _v.at(vertex)->getNum(); }
MVertex* getMeshVertex(int vertex) const { return _v.at(vertex); }
- void clearSortedVertices() { _v.clear(); }
void findBdElement(int i, std::vector<MVertex*>& vertices) const;
int getNumBdElements() const;
- int findBdCellOrientation(Cell* cell) const;
+ int findBdCellOrientation(Cell* cell, int i) const;
- void revertGlobalNum() { _globalNum--; }
void increaseGlobalNum() { _globalNum++; }
// restores the cell information to its original state before reduction
diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 109e4a56f99eda5406401ffece4a597a2111222e..dfc6ff9a5defadff219a358da367609c39fa426e 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -37,11 +37,13 @@ bool CellComplex::_insertCells(std::vector<MElement*>& elements,
for(unsigned int i=0; i < elements.size(); i++){
MElement* element = elements.at(i);
int type = element->getType();
- if(type != TYPE_PNT && type != TYPE_LIN &&
- type != TYPE_TRI && type != TYPE_TET) {
+ if(type == TYPE_PYR || type == TYPE_PRI ||
+ type == TYPE_POLYG || type == TYPE_POLYH) {
Msg::Error("Mesh element type %d not implemented in homology solver", type);
return false;
}
+ if(type == TYPE_QUA || type == TYPE_HEX)
+ _simplicial = false;
Cell* cell = new Cell(element, domain);
bool insert = _cells[cell->getDim()].insert(cell).second;
if(!insert) delete cell;
@@ -50,8 +52,7 @@ bool CellComplex::_insertCells(std::vector<MElement*>& elements,
for (int dim = 3; dim > 0; dim--){
for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
Cell* cell = *cit;
- int numBdElements = cell->getNumBdElements();
- for(int i = 0; i < numBdElements; i++){
+ for(int i = 0; i < cell->getNumBdElements(); i++){
Cell* newCell = new Cell(cell, i);
std::pair<citer, bool> insert =
_cells[newCell->getDim()].insert(newCell);
@@ -60,7 +61,7 @@ bool CellComplex::_insertCells(std::vector<MElement*>& elements,
newCell = *(insert.first);
}
if(domain == 0) {
- int ori = cell->findBdCellOrientation(newCell);
+ int ori = cell->findBdCellOrientation(newCell, i);
cell->addBoundaryCell( ori, newCell, true);
}
}
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index 0da4364848d02e040f5803ecea42ee1389aa33ca..93f4e3237f7ace378b7cafe837846ca80d3abc09 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -34,15 +34,14 @@ ChainComplex::ChainComplex(CellComplex* cellComplex, int domain)
for(CellComplex::citer cit = cellComplex->firstCell(dim);
cit != cellComplex->lastCell(dim); cit++){
Cell* cell = *cit;
- cell->setIndex(0);
cols--;
if((domain == 0 && !cell->inSubdomain()) || domain == 1
|| (domain == 2 && cell->inSubdomain()) ){
cols++;
- cell->setIndex(index);
- index++;
- _cellIndices[dim][cell] = cell->getIndex();
+_cellIndices[dim][cell] = index;
+ index++;
}
+ else _cellIndices[dim][cell] = 0;
}
if(dim > 0) rows = lastCols;
lastCols = cols;
@@ -73,23 +72,24 @@ ChainComplex::ChainComplex(CellComplex* cellComplex, int domain)
if((domain == 0 && !bdCell->inSubdomain()) || domain == 1
|| (domain == 2 && cell->inSubdomain()) ){
int old_elem = 0;
-
- if(bdCell->getIndex() > (int)gmp_matrix_rows( _HMatrix[dim])
- || bdCell->getIndex() < 1
- || cell->getIndex() > (int)gmp_matrix_cols( _HMatrix[dim])
- || cell->getIndex() < 1){
+ int bdCellIndex = getCellIndex(bdCell);
+ int cellIndex = getCellIndex(cell);
+ if(bdCellIndex > (int)gmp_matrix_rows( _HMatrix[dim])
+ || bdCellIndex < 1
+ || cellIndex > (int)gmp_matrix_cols( _HMatrix[dim])
+ || cellIndex < 1){
Msg::Debug("Index out of bound! HMatrix: %d", dim);
}
else{
- gmp_matrix_get_elem(elem, bdCell->getIndex(),
- cell->getIndex(), _HMatrix[dim]);
+ gmp_matrix_get_elem(elem, bdCellIndex,
+ cellIndex, _HMatrix[dim]);
old_elem = mpz_get_si(elem);
mpz_set_si(elem, old_elem + it->second.get());
if( abs((old_elem + it->second.get())) > 1){
//printf("Incidence index: %d, in HMatrix: %d. \n", (old_elem + (*it).second), dim);
}
- gmp_matrix_set_elem(elem, bdCell->getIndex(),
- cell->getIndex(), _HMatrix[dim]);
+ gmp_matrix_set_elem(elem, bdCellIndex,
+ cellIndex, _HMatrix[dim]);
}
}
}
@@ -746,7 +746,7 @@ void HomologySequence::findIcMaps()
cit != _complex->lastCell(i); cit++){
Cell* cell = cit->first;
int row = cit->second;
- int col = _subcomplex->cellIndex(cell);
+ int col = _subcomplex->getCellIndex(cell);
//printf("row %d, col %d. \n", row, col);
if(col != 0) gmp_matrix_set_elem(one, row, col, _Ic_sub[i]);
}
@@ -761,7 +761,7 @@ void HomologySequence::findIcMaps()
cit != _complex->lastCell(i); cit++){
Cell* cell = cit->first;
int row = cit->second;
- int col = _relcomplex->cellIndex(cell);
+ int col = _relcomplex->getCellIndex(cell);
//printf("row %d, col %d. \n", row, col);
if(col != 0) gmp_matrix_set_elem(one, row, col, _Ic_rel[i]);
}
diff --git a/Geo/ChainComplex.h b/Geo/ChainComplex.h
index 03f8b15f9293a821707f55041fa0464f499a180d..45a7b3e9b0a5e1f25d0d457eb51c3e057402fa94 100644
--- a/Geo/ChainComplex.h
+++ b/Geo/ChainComplex.h
@@ -18,7 +18,7 @@
#include <queue>
#include "CellComplex.h"
-#if defined(HAVE_GMP)
+#if defined(HAVE_GMP)
#include "gmp.h"
#include "gmp_normal_form.h"
#else
@@ -30,7 +30,7 @@ class CellComplex;
// A class representing a chain complex of a cell complex.
// This should only be constructed for a reduced cell complex because of
-// dense matrix representations and great computational complexity in
+// dense matrix representations and great computational complexity in
// its methods.
class ChainComplex
{
@@ -38,12 +38,12 @@ class ChainComplex
// boundary operator matrices for this chain complex
// h_k: C_k -> C_(k-1)
gmp_matrix* _HMatrix[5];
-
+
// Basis matrices for the kernels and codomains of the boundary operator
// matrices
gmp_matrix* _kerH[5];
gmp_matrix* _codH[5];
-
+
// matrix of the mapping B_k -> Z_k
gmp_matrix* _JMatrix[5];
// matrix of the mapping H_k -> Z_k
@@ -54,68 +54,68 @@ class ChainComplex
// torsion coefficients of homology generators
// corresponding the columns of _Hbasis
std::vector<long int> _torsion[5];
-
+
int _dim;
CellComplex* _cellComplex;
- // index to cell map
+ // index to cell map
// matrix indices correspond to these cells in _cellComplex
std::map<Cell*, int, Less_Cell> _cellIndices[4];
-
+
// set the matrices
- void setHMatrix(int dim, gmp_matrix* matrix) {
+ void setHMatrix(int dim, gmp_matrix* matrix) {
if(dim > -1 && dim < 5) _HMatrix[dim] = matrix;}
- void setKerHMatrix(int dim, gmp_matrix* matrix) {
+ void setKerHMatrix(int dim, gmp_matrix* matrix) {
if(dim > -1 && dim < 5) _kerH[dim] = matrix;}
- void setCodHMatrix(int dim, gmp_matrix* matrix) {
+ void setCodHMatrix(int dim, gmp_matrix* matrix) {
if(dim > -1 && dim < 5) _codH[dim] = matrix;}
- void setJMatrix(int dim, gmp_matrix* matrix) {
+ void setJMatrix(int dim, gmp_matrix* matrix) {
if(dim > -1 && dim < 5) _JMatrix[dim] = matrix;}
- void setQMatrix(int dim, gmp_matrix* matrix) {
+ void setQMatrix(int dim, gmp_matrix* matrix) {
if(dim > -1 && dim < 5) _QMatrix[dim] = matrix;}
- void setHbasis(int dim, gmp_matrix* matrix) {
+ void setHbasis(int dim, gmp_matrix* matrix) {
if(dim > -1 && dim < 5) _Hbasis[dim] = matrix;}
// get the boundary operator matrix dim->dim-1
- gmp_matrix* getHMatrix(int dim) const {
+ gmp_matrix* getHMatrix(int dim) const {
if(dim > -1 && dim < 5) return _HMatrix[dim]; else return NULL;}
- gmp_matrix* getKerHMatrix(int dim) const {
+ gmp_matrix* getKerHMatrix(int dim) const {
if(dim > -1 && dim < 5) return _kerH[dim]; else return NULL;}
- gmp_matrix* getCodHMatrix(int dim) const {
+ gmp_matrix* getCodHMatrix(int dim) const {
if(dim > -1 && dim < 5) return _codH[dim]; else return NULL;}
- gmp_matrix* getJMatrix(int dim) const {
+ gmp_matrix* getJMatrix(int dim) const {
if(dim > -1 && dim < 5) return _JMatrix[dim]; else return NULL;}
- gmp_matrix* getQMatrix(int dim) const {
+ gmp_matrix* getQMatrix(int dim) const {
if(dim > -1 && dim < 5) return _QMatrix[dim]; else return NULL;}
- gmp_matrix* getHbasis(int dim) const {
+ gmp_matrix* getHbasis(int dim) const {
if(dim > -1 && dim < 5) return _Hbasis[dim]; else return NULL;}
// local deformation tools for chains
bool deformChain(std::map<Cell*, int, Less_Cell>& cells,
std::pair<Cell*, int> cell, bool bend);
bool deform(std::map<Cell*, int, Less_Cell>& cells,
- std::map<Cell*, int, Less_Cell>& cellsInChain,
+ std::map<Cell*, int, Less_Cell>& cellsInChain,
std::map<Cell*, int, Less_Cell>& cellsNotInChain);
void smoothenChain(std::map<Cell*, int, Less_Cell>& cells);
void eraseNullCells(std::map<Cell*, int, Less_Cell>& cells);
void deImmuneCells(std::map<Cell*, int, Less_Cell>& cells);
-
- public:
+
+ public:
// domain = 0 : relative chain space
// domain = 1 : absolute chain space of all cells in cellComplex
// domain = 2 : absolute chain space of cells in subdomain
- ChainComplex(CellComplex* cellComplex, int domain=0);
+ ChainComplex(CellComplex* cellComplex, int domain=0);
~ChainComplex();
-
+
int getDim() const { return _dim; }
-
+
// 0 : C basis (chains)
// 1 : Z basis (cycles)
// 2 : B basis (boundaries)
// 3 : H basis (homology)
// get the bases for various spaces
- gmp_matrix* getBasis(int dim, int basis);
+ gmp_matrix* getBasis(int dim, int basis);
gmp_matrix* getBoundaryOp(int dim) {
if(dim > -1 && dim < 5) return _HMatrix[dim]; else return NULL;}
@@ -123,30 +123,30 @@ class ChainComplex
// of dimension dim (ie. ker(h_dim) and cod(h_dim) )
void KerCod(int dim);
// compute matrix representation J for inclusion relation from dim-cells
- // who are boundary of dim+1-cells to cycles of dim-cells
+ // who are boundary of dim+1-cells to cycles of dim-cells
// (ie. j: cod(h_(dim+1)) -> ker(h_dim) )
void Inclusion(int lowDim, int highDim);
// compute quotient problem for given inclusion relation j to find
- // representatives of homology group generators and possible
+ // representatives of homology group generators and possible
// torsion coeffcients
void Quotient(int dim);
-
- // transpose the boundary operator matrices, these are boundary operator
+
+ // transpose the boundary operator matrices, these are boundary operator
// matrices for the dual mesh
- void transposeHMatrices() {
+ void transposeHMatrices() {
for(int i = 0; i < 5; i++) gmp_matrix_transp(_HMatrix[i]); }
- void transposeHMatrix(int dim) {
+ void transposeHMatrix(int dim) {
if(dim > -1 && dim < 5) gmp_matrix_transp(_HMatrix[dim]); }
-
- // Compute bases for the homology groups of this chain complex
+
+ // Compute bases for the homology groups of this chain complex
void computeHomology(bool dual=false);
-
-
+
+
typedef std::map<Cell*, int, Less_Cell>::iterator citer;
- citer firstCell(int dim){ return _cellIndices[dim].begin(); }
- citer lastCell(int dim){ return _cellIndices[dim].end(); }
+ citer firstCell(int dim){ return _cellIndices[dim].begin(); }
+ citer lastCell(int dim){ return _cellIndices[dim].end(); }
// get the cell index
- int cellIndex(Cell* cell){
+ int getCellIndex(Cell* cell){
citer cit = _cellIndices[cell->getDim()].find(cell);
if(cit != lastCell(cell->getDim())) return cit->second;
else return 0;
@@ -158,13 +158,13 @@ class ChainComplex
// (deform: with local deformations to make chain smoother and to have
// smaller support, deformed chain is homologous to the old one,
// only works for chains of the primary chain complex)
- void getBasisChain(std::map<Cell*, int, Less_Cell>& chain,
+ void getBasisChain(std::map<Cell*, int, Less_Cell>& chain,
int num, int dim, int basis, bool deform=false);
// get rank of a basis
int getBasisSize(int dim, int basis);
- // homology torsion coefficient for dim-dimensional chain num
+ // homology torsion coefficient for dim-dimensional chain num
int getTorsion(int dim, int num);
-
+
// apply a transformation T to a basis (T should be unimodular)
void transformBasis(gmp_matrix* T, int dim, int basis){
if(basis == 3 && _Hbasis[dim] != NULL) {
@@ -174,10 +174,10 @@ class ChainComplex
//void printBasisChain(std::map<Cell*, int, Less_Cell>& chain);
// debugging aid
- int printMatrix(gmp_matrix* matrix){
- printf("%d rows and %d columns\n",
- (int)gmp_matrix_rows(matrix), (int)gmp_matrix_cols(matrix));
- return gmp_matrix_printf(matrix); }
+ int printMatrix(gmp_matrix* matrix){
+ printf("%d rows and %d columns\n",
+ (int)gmp_matrix_rows(matrix), (int)gmp_matrix_cols(matrix));
+ return gmp_matrix_printf(matrix); }
void matrixTest();
};
@@ -191,7 +191,7 @@ class HomologySequence
ChainComplex* _subcomplex;
ChainComplex* _complex;
ChainComplex* _relcomplex;
-
+
gmp_matrix* _Ic_sub[4];
gmp_matrix* _Ic_rel[4];
@@ -199,7 +199,7 @@ class HomologySequence
gmp_matrix* _Jh[4];
gmp_matrix* _invIh[4];
gmp_matrix* _invJh[4];
-
+
gmp_matrix* _Dh[4];
gmp_matrix* _invDh[4];
@@ -213,29 +213,29 @@ class HomologySequence
void findInvDhMap(int i);
public:
-
- HomologySequence(ChainComplex* subcomplex, ChainComplex* complex,
+
+ HomologySequence(ChainComplex* subcomplex, ChainComplex* complex,
ChainComplex* relcomplex);
~HomologySequence();
// create an inclusion map from domBasis to codBasis
// (deletes domBasis, leaves codBasis unaffected)
- gmp_matrix* createIncMap(gmp_matrix* domBasis,
+ gmp_matrix* createIncMap(gmp_matrix* domBasis,
gmp_matrix* codBasis);
gmp_matrix* removeZeroCols(gmp_matrix* matrix);
- void blockHBasis(gmp_matrix* block1T, gmp_matrix* block2T,
+ void blockHBasis(gmp_matrix* block1T, gmp_matrix* block2T,
ChainComplex* complex, int dim);
-
+
int printMatrix(gmp_matrix* matrix){
if(matrix == NULL){ printf("NULL matrix. \n"); return 0; }
printf("%d rows and %d columns\n",
(int)gmp_matrix_rows(matrix), (int)gmp_matrix_cols(matrix));
- return gmp_matrix_printf(matrix); }
-
+ return gmp_matrix_printf(matrix); }
+
};
#endif
-#endif
-
+#endif
+
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index d635a22e3b6ba06f97a313dd52db547d8e047113..9596a1715d9b90d05194100b59259ce72c1473f5 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -313,6 +313,9 @@ void Homology::_createPhysicalGroup(std::map<Cell*, int, Less_Cell>& chain, std:
MElementFactory factory;
int dim = 0;
+ std::map<Cell*, int, Less_Cell> chain2;
+ std::string name2 = name + "i";
+
typedef std::map<Cell*, int, Less_Cell>::iterator citer;
for(citer cit = chain.begin(); cit != chain.end(); cit++){
Cell* cell = cit->first;
@@ -335,8 +338,15 @@ void Homology::_createPhysicalGroup(std::map<Cell*, int, Less_Cell>& chain, std:
elements.push_back(ecopy);
}
+ //std::vector<double> coeffs (1,1);
std::vector<double> coeffs (1,abs(coeff));
data[e->getNum()] = coeffs;
+
+ /*if(abs(coeff) > 1) {
+ if(coeff < 0) chain2[cell] = coeff+1;
+ else chain2[cell] = coeff-1;
+ }*/
+
}
GModel* m = this->getModel();
@@ -373,6 +383,7 @@ void Homology::_createPhysicalGroup(std::map<Cell*, int, Less_Cell>& chain, std:
view->setOptions(opt);
#endif
}
+ //if(!chain2.empty()) _createPhysicalGroup(chain2, name2);
}
bool Homology::writeBasisMSH(bool binary)