Skip to content
Snippets Groups Projects
Commit e88541e7 authored by Matti Pellika's avatar Matti Pellika
Browse files

Homology solver to work with quad and hex meshes, and some cleanup

parent ba4acffd
Branches
Tags
No related merge requests found
Show whitespace changes
Inline Side-by-side
Geo/Cell.cpp
+ 218
− 127
View file @ e88541e7
...@@ -6,6 +6,12 @@ ...@@ -6,6 +6,12 @@
// Contributed by Matti Pellikka <matti.pellikka@tut.fi>. // Contributed by Matti Pellikka <matti.pellikka@tut.fi>.
#include "Cell.h" #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 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 ...@@ -23,35 +29,26 @@ bool Less_Cell::operator()(const Cell* c1, const Cell* c2) const
return false; return false;
} }
int Cell::_globalNum = 0; bool Less_VertexNumIndex::operator()(const std::pair<MVertex*, int> v1,
Cell::Cell(int dim, std::vector<MVertex*>& v) const std::pair<MVertex*, int> v2) const
{ {
_dim = dim; return (v1.first->getNum() < v2.first->getNum());
_domain = 0;
_combined = false;
_immune = false;
_v = v;
std::sort(_v.begin(), _v.end(), MVertexLessThanNum());
_num = 0;
_index = 0;
} }
int Cell::_globalNum = 0;
Cell::Cell(MElement* element, int domain) Cell::Cell(MElement* element, int domain)
{ {
_dim = element->getDim(); _dim = element->getDim();
_domain = domain; _domain = domain;
_combined = false; _combined = false;
_immune = 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)); _v.push_back(element->getVertex(i));
}
std::sort(_v.begin(), _v.end(), MVertexLessThanNum());
_num = 0;
_index = 0;
_sortVertexIndices();
} }
Cell::Cell(Cell* parent, int i) Cell::Cell(Cell* parent, int i)
...@@ -60,144 +57,240 @@ Cell::Cell(Cell* parent, int i) ...@@ -60,144 +57,240 @@ Cell::Cell(Cell* parent, int i)
_domain = parent->getDomain(); _domain = parent->getDomain();
_combined = false; _combined = false;
_immune = false; _immune = false;
_num = 0;
parent->findBdElement(i, _v); parent->findBdElement(i, _v);
std::sort(_v.begin(), _v.end(), MVertexLessThanNum()); _sortVertexIndices();
}
_num = 0; void Cell::_sortVertexIndices()
_index = 0; {
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 void Cell::findBdElement(int i, std::vector<MVertex*>& vertices) const
{ {
vertices.clear(); vertices.clear();
if(_dim == 1){ // i = 0 -> ori = -1, i = 1 -> ori = 1 switch (_dim) {
case 1:
vertices.push_back(_v[i]); vertices.push_back(_v[i]);
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(_dim == 2){ case 3:
if(i == 0){ // ori = 1 switch (getNumVertices()) {
vertices.push_back(_v[0]); case 4:
vertices.push_back(_v[1]); for(int j = 0; j < 3; j++)
} vertices.push_back(_v[MTetrahedron::faces_tetra(i, j)]);
else if(i == 1){ // ori = 1 return;
vertices.push_back(_v[1]); case 5:
vertices.push_back(_v[2]); if(i < 3) {
} for(int j = 0; j < 3; j++)
else if(i == 2){ // ori = 1 vertices.push_back(_v[MPyramid::faces_pyramid(i, j)]);
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){ else {
vertices.push_back(_v[0]); vertices.push_back(_v[0]);
vertices.push_back(_v[3]); vertices.push_back(_v[3]);
vertices.push_back(_v[2]); vertices.push_back(_v[2]);
}
else if(i == 3){
vertices.push_back(_v[3]);
vertices.push_back(_v[1]); vertices.push_back(_v[1]);
vertices.push_back(_v[2]);
} }
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 int Cell::getNumBdElements() const
{ {
if(_dim == 0) return 0; switch (_dim) {
else if(_dim == 1) return 2; case 0: return 0;
else if(_dim == 2) return 3; case 1: return 2;
else if(_dim == 3) return 4; case 2:
else return 0; switch (getNumVertices()) {
case 3: return 3;
case 4: return 4;
default: return 0;
} }
int Cell::findBdCellOrientation(Cell* cell) const 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, int i) const
{ {
std::vector<MVertex*> v; std::vector<MVertex*> v;
cell->getMeshVertices(v); cell->getMeshVertices(v);
if(_dim == 0) return 0; switch (_dim) {
else if(_dim == 1){ case 0: return 0;
case 1:
if(v[0]->getNum() == _v[0]->getNum()) return -1; if(v[0]->getNum() == _v[0]->getNum()) return -1;
else if(v[0]->getNum() == _v[1]->getNum()) return 1; else if(v[0]->getNum() == _v[1]->getNum()) return 1;
}
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; 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;
}
} }
int Cell::getTypeMSH() const int Cell::getTypeMSH() const
{ {
if(_dim == 0) return 15; switch (_dim) {
else if(_dim == 1) return 1; case 0: return MSH_PNT;
else if(_dim == 2) return 2; case 1: return MSH_LIN_2;
else if(_dim == 3) return 4; case 2:
else return 0; 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 bool Cell::hasVertex(int vertex) const
{ {
//_vs
std::vector<int> v; std::vector<int> v;
for(unsigned int i = 0; i < _v.size(); i++) { 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(), std::vector<int>::const_iterator it = std::find(v.begin(), v.end(),
vertex); vertex);
...@@ -219,7 +312,7 @@ void Cell::printCell() ...@@ -219,7 +312,7 @@ void Cell::printCell()
printf("%d-cell %d: \n" , getDim(), getNum()); printf("%d-cell %d: \n" , getDim(), getNum());
printf("Vertices: "); printf("Vertices: ");
for(int i = 0; i < this->getNumVertices(); i++){ 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(", in subdomain: %d, ", inSubdomain());
printf("combined: %d. \n" , isCombined() ); printf("combined: %d. \n" , isCombined() );
...@@ -239,7 +332,6 @@ void Cell::restoreCell(){ ...@@ -239,7 +332,6 @@ void Cell::restoreCell(){
} }
for(unsigned int i = 0; i < toRemove.size(); i++) _bd.erase(toRemove[i]); for(unsigned int i = 0; i < toRemove.size(); i++) _bd.erase(toRemove[i]);
_combined = false; _combined = false;
_index = 0;
_immune = false; _immune = false;
} }
...@@ -351,7 +443,7 @@ void Cell::printCoboundary(bool orig) ...@@ -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 // use "smaller" cell as c2
if(c1->getNumCells() < c2->getNumCells()){ if(c1->getNumCells() < c2->getNumCells()){
...@@ -361,7 +453,6 @@ CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell() ...@@ -361,7 +453,6 @@ CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell()
} }
_num = ++_globalNum; _num = ++_globalNum;
_index = c1->getIndex();
_domain = c1->getDomain(); _domain = c1->getDomain();
_combined = true; _combined = true;
_immune = false; _immune = false;
...@@ -417,12 +508,12 @@ CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell() ...@@ -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; _num = ++_globalNum;
_index = cells.at(0)->getIndex();
_domain = cells.at(0)->getDomain(); _domain = cells.at(0)->getDomain();
_combined = true; _combined = true;
_immune = false;
// cells // cells
for(unsigned int i = 0; i < cells.size(); i++){ for(unsigned int i = 0; i < cells.size(); i++){
......
Geo/Cell.h
+ 20
− 20
View file @ e88541e7
...@@ -19,11 +19,17 @@ public: ...@@ -19,11 +19,17 @@ public:
bool operator()(const Cell* c1, const Cell* c2) const; 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 to save cell boundary orientation information
class BdInfo { class BdInfo {
private: private:
short int _ori; signed char _ori;
short int _origOri; signed char _origOri;
public: public:
BdInfo(int ori) { _ori = ori; _origOri = 0; } BdInfo(int ori) { _ori = ori; _origOri = 0; }
...@@ -41,15 +47,13 @@ class BdInfo { ...@@ -41,15 +47,13 @@ class BdInfo {
class Cell { class Cell {
protected: protected:
static int _globalNum; static int _globalNum;
int _num; int _num;
// mutable index for each cell (used to create boundary operator matrices) char _domain;
int _index;
int _domain;
// whether this cell a combinded cell of elemetary cells // whether this cell a combinded cell of elementary cells
bool _combined; bool _combined;
// for some algorithms to omit this cell // for some algorithms to omit this cell
bool _immune; bool _immune;
...@@ -60,15 +64,16 @@ class Cell { ...@@ -60,15 +64,16 @@ class Cell {
private: private:
int _dim; char _dim;
// sorted vertices of this cell (used for ordering of the cells)
std::vector<MVertex*> _v; 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); public:
Cell() : _num(0), _dim(0), _index(0), _domain(0), _combined(false), _immune(false) {}
Cell() {}
Cell(MElement* element, int domain); Cell(MElement* element, int domain);
Cell(Cell* parent, int i); Cell(Cell* parent, int i);
...@@ -77,26 +82,21 @@ class Cell { ...@@ -77,26 +82,21 @@ class Cell {
void setNum(int num) { _num = num; }; void setNum(int num) { _num = num; };
int getTypeMSH() const; int getTypeMSH() const;
virtual int getDim() const { return _dim; } 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; } 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; }; void setImmune(bool immune) { _immune = immune; };
bool getImmune() const { return _immune; }; bool getImmune() const { return _immune; };
int getNumSortedVertices() const { return _v.size(); } 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 getNumVertices() const { return _v.size(); }
int getVertex(int vertex) const { return _v.at(vertex)->getNum(); }
MVertex* getMeshVertex(int vertex) const { return _v.at(vertex); } MVertex* getMeshVertex(int vertex) const { return _v.at(vertex); }
void clearSortedVertices() { _v.clear(); }
void findBdElement(int i, std::vector<MVertex*>& vertices) const; void findBdElement(int i, std::vector<MVertex*>& vertices) const;
int getNumBdElements() const; int getNumBdElements() const;
int findBdCellOrientation(Cell* cell) const; int findBdCellOrientation(Cell* cell, int i) const;
void revertGlobalNum() { _globalNum--; }
void increaseGlobalNum() { _globalNum++; } void increaseGlobalNum() { _globalNum++; }
// restores the cell information to its original state before reduction // restores the cell information to its original state before reduction
......
Geo/CellComplex.cpp
+ 6
− 5
View file @ e88541e7
...@@ -37,11 +37,13 @@ bool CellComplex::_insertCells(std::vector<MElement*>& elements, ...@@ -37,11 +37,13 @@ bool CellComplex::_insertCells(std::vector<MElement*>& elements,
for(unsigned int i=0; i < elements.size(); i++){ for(unsigned int i=0; i < elements.size(); i++){
MElement* element = elements.at(i); MElement* element = elements.at(i);
int type = element->getType(); int type = element->getType();
if(type != TYPE_PNT && type != TYPE_LIN && if(type == TYPE_PYR || type == TYPE_PRI ||
type != TYPE_TRI && type != TYPE_TET) { type == TYPE_POLYG || type == TYPE_POLYH) {
Msg::Error("Mesh element type %d not implemented in homology solver", type); Msg::Error("Mesh element type %d not implemented in homology solver", type);
return false; return false;
} }
if(type == TYPE_QUA || type == TYPE_HEX)
_simplicial = false;
Cell* cell = new Cell(element, domain); Cell* cell = new Cell(element, domain);
bool insert = _cells[cell->getDim()].insert(cell).second; bool insert = _cells[cell->getDim()].insert(cell).second;
if(!insert) delete cell; if(!insert) delete cell;
...@@ -50,8 +52,7 @@ bool CellComplex::_insertCells(std::vector<MElement*>& elements, ...@@ -50,8 +52,7 @@ bool CellComplex::_insertCells(std::vector<MElement*>& elements,
for (int dim = 3; dim > 0; dim--){ for (int dim = 3; dim > 0; dim--){
for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){ for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
Cell* cell = *cit; Cell* cell = *cit;
int numBdElements = cell->getNumBdElements(); for(int i = 0; i < cell->getNumBdElements(); i++){
for(int i = 0; i < numBdElements; i++){
Cell* newCell = new Cell(cell, i); Cell* newCell = new Cell(cell, i);
std::pair<citer, bool> insert = std::pair<citer, bool> insert =
_cells[newCell->getDim()].insert(newCell); _cells[newCell->getDim()].insert(newCell);
...@@ -60,7 +61,7 @@ bool CellComplex::_insertCells(std::vector<MElement*>& elements, ...@@ -60,7 +61,7 @@ bool CellComplex::_insertCells(std::vector<MElement*>& elements,
newCell = *(insert.first); newCell = *(insert.first);
} }
if(domain == 0) { if(domain == 0) {
int ori = cell->findBdCellOrientation(newCell); int ori = cell->findBdCellOrientation(newCell, i);
cell->addBoundaryCell( ori, newCell, true); cell->addBoundaryCell( ori, newCell, true);
} }
} }
......
Geo/ChainComplex.cpp
+ 15
− 15
View file @ e88541e7
...@@ -34,15 +34,14 @@ ChainComplex::ChainComplex(CellComplex* cellComplex, int domain) ...@@ -34,15 +34,14 @@ ChainComplex::ChainComplex(CellComplex* cellComplex, int domain)
for(CellComplex::citer cit = cellComplex->firstCell(dim); for(CellComplex::citer cit = cellComplex->firstCell(dim);
cit != cellComplex->lastCell(dim); cit++){ cit != cellComplex->lastCell(dim); cit++){
Cell* cell = *cit; Cell* cell = *cit;
cell->setIndex(0);
cols--; cols--;
if((domain == 0 && !cell->inSubdomain()) || domain == 1 if((domain == 0 && !cell->inSubdomain()) || domain == 1
|| (domain == 2 && cell->inSubdomain()) ){ || (domain == 2 && cell->inSubdomain()) ){
cols++; cols++;
cell->setIndex(index); _cellIndices[dim][cell] = index;
index++; index++;
_cellIndices[dim][cell] = cell->getIndex();
} }
else _cellIndices[dim][cell] = 0;
} }
if(dim > 0) rows = lastCols; if(dim > 0) rows = lastCols;
lastCols = cols; lastCols = cols;
...@@ -73,23 +72,24 @@ ChainComplex::ChainComplex(CellComplex* cellComplex, int domain) ...@@ -73,23 +72,24 @@ ChainComplex::ChainComplex(CellComplex* cellComplex, int domain)
if((domain == 0 && !bdCell->inSubdomain()) || domain == 1 if((domain == 0 && !bdCell->inSubdomain()) || domain == 1
|| (domain == 2 && cell->inSubdomain()) ){ || (domain == 2 && cell->inSubdomain()) ){
int old_elem = 0; int old_elem = 0;
int bdCellIndex = getCellIndex(bdCell);
if(bdCell->getIndex() > (int)gmp_matrix_rows( _HMatrix[dim]) int cellIndex = getCellIndex(cell);
|| bdCell->getIndex() < 1 if(bdCellIndex > (int)gmp_matrix_rows( _HMatrix[dim])
|| cell->getIndex() > (int)gmp_matrix_cols( _HMatrix[dim]) || bdCellIndex < 1
|| cell->getIndex() < 1){ || cellIndex > (int)gmp_matrix_cols( _HMatrix[dim])
|| cellIndex < 1){
Msg::Debug("Index out of bound! HMatrix: %d", dim); Msg::Debug("Index out of bound! HMatrix: %d", dim);
} }
else{ else{
gmp_matrix_get_elem(elem, bdCell->getIndex(), gmp_matrix_get_elem(elem, bdCellIndex,
cell->getIndex(), _HMatrix[dim]); cellIndex, _HMatrix[dim]);
old_elem = mpz_get_si(elem); old_elem = mpz_get_si(elem);
mpz_set_si(elem, old_elem + it->second.get()); mpz_set_si(elem, old_elem + it->second.get());
if( abs((old_elem + it->second.get())) > 1){ if( abs((old_elem + it->second.get())) > 1){
//printf("Incidence index: %d, in HMatrix: %d. \n", (old_elem + (*it).second), dim); //printf("Incidence index: %d, in HMatrix: %d. \n", (old_elem + (*it).second), dim);
} }
gmp_matrix_set_elem(elem, bdCell->getIndex(), gmp_matrix_set_elem(elem, bdCellIndex,
cell->getIndex(), _HMatrix[dim]); cellIndex, _HMatrix[dim]);
} }
} }
} }
...@@ -746,7 +746,7 @@ void HomologySequence::findIcMaps() ...@@ -746,7 +746,7 @@ void HomologySequence::findIcMaps()
cit != _complex->lastCell(i); cit++){ cit != _complex->lastCell(i); cit++){
Cell* cell = cit->first; Cell* cell = cit->first;
int row = cit->second; int row = cit->second;
int col = _subcomplex->cellIndex(cell); int col = _subcomplex->getCellIndex(cell);
//printf("row %d, col %d. \n", row, col); //printf("row %d, col %d. \n", row, col);
if(col != 0) gmp_matrix_set_elem(one, row, col, _Ic_sub[i]); if(col != 0) gmp_matrix_set_elem(one, row, col, _Ic_sub[i]);
} }
...@@ -761,7 +761,7 @@ void HomologySequence::findIcMaps() ...@@ -761,7 +761,7 @@ void HomologySequence::findIcMaps()
cit != _complex->lastCell(i); cit++){ cit != _complex->lastCell(i); cit++){
Cell* cell = cit->first; Cell* cell = cit->first;
int row = cit->second; int row = cit->second;
int col = _relcomplex->cellIndex(cell); int col = _relcomplex->getCellIndex(cell);
//printf("row %d, col %d. \n", row, col); //printf("row %d, col %d. \n", row, col);
if(col != 0) gmp_matrix_set_elem(one, row, col, _Ic_rel[i]); if(col != 0) gmp_matrix_set_elem(one, row, col, _Ic_rel[i]);
} }
......
Geo/ChainComplex.h
+ 57
− 57
View file @ e88541e7
...@@ -146,7 +146,7 @@ class ChainComplex ...@@ -146,7 +146,7 @@ class ChainComplex
citer firstCell(int dim){ return _cellIndices[dim].begin(); } citer firstCell(int dim){ return _cellIndices[dim].begin(); }
citer lastCell(int dim){ return _cellIndices[dim].end(); } citer lastCell(int dim){ return _cellIndices[dim].end(); }
// get the cell index // get the cell index
int cellIndex(Cell* cell){ int getCellIndex(Cell* cell){
citer cit = _cellIndices[cell->getDim()].find(cell); citer cit = _cellIndices[cell->getDim()].find(cell);
if(cit != lastCell(cell->getDim())) return cit->second; if(cit != lastCell(cell->getDim())) return cit->second;
else return 0; else return 0;
......
Geo/Homology.cpp
+ 11
− 0
View file @ e88541e7
...@@ -313,6 +313,9 @@ void Homology::_createPhysicalGroup(std::map<Cell*, int, Less_Cell>& chain, std: ...@@ -313,6 +313,9 @@ void Homology::_createPhysicalGroup(std::map<Cell*, int, Less_Cell>& chain, std:
MElementFactory factory; MElementFactory factory;
int dim = 0; int dim = 0;
std::map<Cell*, int, Less_Cell> chain2;
std::string name2 = name + "i";
typedef std::map<Cell*, int, Less_Cell>::iterator citer; typedef std::map<Cell*, int, Less_Cell>::iterator citer;
for(citer cit = chain.begin(); cit != chain.end(); cit++){ for(citer cit = chain.begin(); cit != chain.end(); cit++){
Cell* cell = cit->first; Cell* cell = cit->first;
...@@ -335,8 +338,15 @@ void Homology::_createPhysicalGroup(std::map<Cell*, int, Less_Cell>& chain, std: ...@@ -335,8 +338,15 @@ void Homology::_createPhysicalGroup(std::map<Cell*, int, Less_Cell>& chain, std:
elements.push_back(ecopy); elements.push_back(ecopy);
} }
//std::vector<double> coeffs (1,1);
std::vector<double> coeffs (1,abs(coeff)); std::vector<double> coeffs (1,abs(coeff));
data[e->getNum()] = coeffs; data[e->getNum()] = coeffs;
/*if(abs(coeff) > 1) {
if(coeff < 0) chain2[cell] = coeff+1;
else chain2[cell] = coeff-1;
}*/
} }
GModel* m = this->getModel(); GModel* m = this->getModel();
...@@ -373,6 +383,7 @@ void Homology::_createPhysicalGroup(std::map<Cell*, int, Less_Cell>& chain, std: ...@@ -373,6 +383,7 @@ void Homology::_createPhysicalGroup(std::map<Cell*, int, Less_Cell>& chain, std:
view->setOptions(opt); view->setOptions(opt);
#endif #endif
} }
//if(!chain2.empty()) _createPhysicalGroup(chain2, name2);
} }
bool Homology::writeBasisMSH(bool binary) bool Homology::writeBasisMSH(bool binary)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment