diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 892fc24c667f1169656bd9f764d670802775495e..dcbad2b1d6657380c71bfdae3b71848cfdf01071 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -69,13 +69,12 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su
}
}
}
-
// insert cells into cell complex
// subdomain need to be inserted first!
insert_cells(true, true);
insert_cells(false, true);
insert_cells(false, false);
-
+
// find mesh vertices in the domain
for(unsigned int j=0; j < domain.size(); j++) {
for(unsigned int i=0; i < domain.at(j)->getNumMeshElements(); i++){
@@ -85,6 +84,7 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su
}
}
}
+
for(unsigned int j=0; j < subdomain.size(); j++) {
for(unsigned int i=0; i < subdomain.at(j)->getNumMeshElements(); i++){
for(int k=0; k < subdomain.at(j)->getMeshElement(i)->getNumVertices(); k++){
@@ -94,7 +94,7 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su
}
}
-
+
// attach individual tags to cells
int tag = 1;
for(int i = 0; i < 4; i++){
@@ -120,11 +120,12 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){
else if(boundary) domain = _boundary;
else domain = _domain;
- std::vector<int> vertices;
+ std::vector<MVertex*> vertices;
int vertex = 0;
std::pair<citer, bool> insertInfo;
+
// add highest dimensional cells
for(unsigned int j=0; j < domain.size(); j++) {
for(unsigned int i=0; i < domain.at(j)->getNumMeshElements(); i++){
@@ -132,7 +133,7 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){
for(int k=0; k < domain.at(j)->getMeshElement(i)->getNumVertices(); k++){
MVertex* vertex = domain.at(j)->getMeshElement(i)->getVertex(k);
- vertices.push_back(vertex->getNum());
+ vertices.push_back(vertex);
}
int dim = domain.at(j)->getMeshElement(i)->getDim();
@@ -141,23 +142,23 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){
if(dim == 3) cell = new ThreeSimplex(vertices, 0, subdomain, boundary);
else if(dim == 2) cell = new TwoSimplex(vertices, 0, subdomain, boundary);
else if(dim == 1) cell = new OneSimplex(vertices, 0, subdomain, boundary);
- else cell = new ZeroSimplex(vertices.at(0), 0, subdomain, boundary);
+ else cell = new ZeroSimplex(vertices, 0, subdomain, boundary);
insertInfo = _cells[dim].insert(cell);
if(!insertInfo.second) delete cell;
}
}
-
+
// add lower dimensional cells recursively
for (int dim = 3; dim > 0; dim--){
for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
Cell* cell = *cit;
- std::vector<int> vertices;
+ std::vector<MVertex*> vertices;
for(int i = 0; i < cell->getNumFacets(); i++){
cell->getFacetVertices(i, vertices);
Cell* newCell;
if(dim == 3) newCell = new TwoSimplex(vertices, 0, subdomain, boundary);
else if(dim == 2) newCell = new OneSimplex(vertices, 0, subdomain, boundary);
- else if(dim == 1) newCell = new ZeroSimplex(vertices.at(0), 0, subdomain, boundary);
+ else if(dim == 1) newCell = new ZeroSimplex(vertices, 0, subdomain, boundary);
insertInfo = _cells[dim-1].insert(newCell);
if(!insertInfo.second){
delete newCell;
@@ -186,7 +187,7 @@ void CellComplex::insertCell(Cell* cell){
int Simplex::kappa(Cell* tau) const{
for(int i=0; i < tau->getNumVertices(); i++){
- if( !(this->hasVertex(tau->getVertex(i))) ) return 0;
+ if( !(this->hasVertex(tau->getVertex(i)->getNum())) ) return 0;
}
if(this->getDim() - tau->getDim() != 1) return 0;
@@ -249,62 +250,6 @@ int OneSimplex::kappa(Cell* tau) const{
}
*/
-int Quadrangle::kappa(Cell* tau) const{
- for(int i=0; i < tau->getNumVertices(); i++){
- if( !(this->hasVertex(tau->getVertex(i))) ) return 0;
- }
-
- if(tau->getDim() != 1) return 0;
-
- int value=1;
-
- std::vector<int> vTau = tau->getVertexVector();
-
- for(int i = 0; i < this->getNumFacets(); i++){
- std::vector<int> v;
- this->getFacetVertices(i, v);
- value = -1;
- do {
- value = value*-1;
- if(v == vTau) return value;
- }
- while (std::next_permutation(v.begin(), v.end()) );
- }
-
- return value;
-}
-
-/*
-std::set<Cell*, Less_Cell>::iterator CellComplex::findCell(int dim, std::vector<int>& vertices, bool original){
- Cell* cell;
-
- if(dim == 0) cell = new ZeroSimplex(vertices.at(0));
- else if(dim == 1) cell = new OneSimplex(vertices);
- else if(dim == 2) cell = new TwoSimplex(vertices);
- else cell = new ThreeSimplex(vertices);
-
- citer cit;
- if(!original) cit = _cells[dim].find(cell);
- else cit = _originalCells[dim].find(cell);
-
- delete cell;
- return cit;
-}
-
-std::set<Cell*, Less_Cell>::iterator CellComplex::findCell(int dim, int vertex, int dummy){
- Cell* cell;
- if(dim == 0) cell = new ZeroSimplex(vertex);
- else if(dim == 1) cell = new OneSimplex(vertex, dummy);
- else if(dim == 2) cell = new TwoSimplex(vertex, dummy);
- else cell = new ThreeSimplex(vertex, dummy);
-
- citer cit = _cells[dim].lower_bound(cell);
-
- delete cell;
- return cit;
-}
-*/
-
void CellComplex::removeCell(Cell* cell){
_cells[cell->getDim()].erase(cell);
@@ -380,79 +325,6 @@ int CellComplex::coreduction(Cell* generator){
//printf("Coreduction: %d loops with %d coreductions\n", round, coreductions);
return coreductions;
}
-/*
-int CellComplex::coreduction(int dim){
-
- if(dim < 0 || dim > 2) return 0;
-
- std::list<Cell*> bd_c;
- int count = 0;
-
- bool coreduced = true;
- while (coreduced){
- coreduced = false;
- for(citer cit = firstCell(dim+1); cit != lastCell(dim+1); cit++){
- Cell* cell = *cit;
-
- bd_c = cell->getBoundary();
- if(bd_c.size() == 1 && inSameDomain(cell, bd_c.front()) ){
- removeCell(cell);
- removeCell(bd_c.front());
- count++;
- coreduced = true;
- }
-
- }
- }
-
- return count;
-
-}
-*/
-/*
-int CellComplex::reduction(Cell* generator){
-
- int count = 0;
-
- std::queue<Cell*> Q;
- std::set<Cell*, Less_Cell> Qset;
-
- Q.push(generator);
- Qset.insert(generator);
-
-
- std::list<Cell*> cbd_s;
- std::list<Cell*> bd_c;
- Cell* s;
- int round = 0;
- while( !Q.empty() ){
- round++;
- //printf("%d ", round);
-
- s = Q.front();
- Q.pop();
- removeCellQset(s, Qset);
-
- cbd_s = s->getCoboundary();
- if( cbd_s.size() == 1 && inSameDomain(s, cbd_s.front()) ){
- removeCell(s);
- bd_c = cbd_s.front()->getBoundary();
- enqueueCells(bd_c, Q, Qset);
- removeCell(cbd_s.front());
- count++;
- }
- else if(cbd_s.empty()){
- bd_c = s->getBoundary();
- enqueueCells(bd_c, Q, Qset);
- }
-
-
- }
-
- return count;
-}
-*/
-
int CellComplex::reduction(int dim){
if(dim < 1 || dim > 3) return 0;
@@ -483,13 +355,6 @@ void CellComplex::reduceComplex(bool omitHighdim){
printf("Cell complex before reduction: %d volumes, %d faces, %d edges and %d vertices.\n",
getSize(3), getSize(2), getSize(1), getSize(0));
-
- //reduction(3);
- //reduction(2);
- //reduction(1);
-
-
-
int count = 0;
for(int i = 3; i > 0; i--) count = count + reduction(i);
@@ -527,45 +392,7 @@ void CellComplex::reduceComplex(bool omitHighdim){
getSize(3), getSize(2), getSize(1), getSize(0));
}
-/*
-void CellComplex::coreduceComplex(int generatorDim){
-
-
- std::set<Cell*, Less_Cell> generatorCells[4];
-
- printf("Cell complex before coreduction: %d volumes, %d faces, %d edges and %d vertices.\n",
- getSize(3), getSize(2), getSize(1), getSize(0));
- for(int i = 0; i < 4; i++){
- while (getSize(i) != 0){
- //if(generatorDim == i || i == generatorDim+1) break;
- citer cit = firstCell(i);
- Cell* cell = *cit;
- while(!cell->inSubdomain() && cit != lastCell(i)){
- cell = *cit;
- cit++;
- }
- if(!cell->inSubdomain()) generatorCells[i].insert(cell);
- removeCell(cell);
- coreduction(cell);
- }
- if(generatorDim == i) break;
-
- }
-
- for(int i = 0; i < 4; i++){
- for(citer cit = generatorCells[i].begin(); cit != generatorCells[i].end(); cit++){
- Cell* cell = *cit;
- _cells[i].insert(cell);
- }
- }
-
- printf("Cell complex after coreduction: %d volumes, %d faces, %d edges and %d vertices.\n",
- getSize(3), getSize(2), getSize(1), getSize(0));
-
- return;
-}
-*/
void CellComplex::removeSubdomain(){
for(int i = 0; i < 4; i++){
@@ -581,55 +408,7 @@ void CellComplex::removeSubdomain(){
return;
}
-/*
-int CellComplex::coreduction(int dim){
-
- int count = 0;
-
- std::queue<Cell*> Q;
- std::set<Cell*, Less_Cell> Qset;
-
- std::list<Cell*> bd_s;
- std::list<Cell*> cbd_c;
-
- for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
-
- Cell* cell = *cit;
- cbd_c = cell->getCoboundary();
- enqueueCells(cbd_c, Q, Qset);
-
- int round = 0;
- while( !Q.empty() ){
- round++;
- //printf("%d ", round);
-
- Cell* s = Q.front();
- Q.pop();
-
- bd_s = s->getBoundary();
-
- if( bd_s.size() == 1 && inSameDomain(s, bd_s.front()) ){
- removeCell(s);
- cbd_c = bd_s.front()->getCoboundary();
- enqueueCells(cbd_c, Q, Qset);
- removeCell(bd_s.front());
- cit = firstCell(dim);
- count++;
- }
- else if(bd_s.empty()){
- cbd_c = s->getCoboundary();
- enqueueCells(cbd_c, Q, Qset);
- }
-
- removeCellQset(s, Qset);
-
- }
-
-
- }
- return count;
-}
-*/
+
void CellComplex::coreduceComplex(bool omitLowdim){
printf("Cell complex before coreduction: %d volumes, %d faces, %d edges and %d vertices.\n",
@@ -668,24 +447,7 @@ void CellComplex::coreduceComplex(bool omitLowdim){
}
- /*
- std::set<Cell*, Less_Cell> generatorCells;
-
- while (getSize(0) != 0){
- citer cit = firstCell(0);
- Cell* cell = *cit;
- while(!cell->inSubdomain() && cit != lastCell(0)){
- cell = *cit;
- cit++;
- }
- if(!cell->inSubdomain()) generatorCells.insert(cell);
- removeCell(cell);
- coreduction(cell);
- }
- */
-
-
printf("Cell complex after coreduction: %d volumes, %d faces, %d edges and %d vertices.\n",
getSize(3), getSize(2), getSize(1), getSize(0));
@@ -694,18 +456,6 @@ void CellComplex::coreduceComplex(bool omitLowdim){
void CellComplex::computeBettiNumbers(){
- /*
- removeSubdomain();
-
- int count = 0;
- for(int dim = 0; dim < 4; dim++){
- citer cit = firstCell(dim);
- if(cit != lastCell(dim)){
- Cell* cell = *cit;
- count = count + coreduction(cell);
- }
- }
- */
for(int i = 0; i < 4; i++){
printf("Betti number computation process: step %d of 4 \n", i+1);
@@ -913,41 +663,6 @@ int CellComplex::combine(int dim){
return count;
}
-/*
-void CellComplex::repairComplex(int i){
-
- printf("Cell complex before repair: %d volumes, %d faces, %d edges and %d vertices.\n",
- getSize(3), getSize(2), getSize(1), getSize(0));
-
- while(i > 0){
-
- for(citer cit = firstCell(i); cit != lastCell(i); cit++){
- Cell* cell = *cit;
- std::vector<int> vertices = cell->getVertexVector();
-
- for(int j=0; j < vertices.size(); j++){
- std::vector<int> bVertices;
-
- for(int k=0; k < vertices.size(); k++){
- if (k !=j ) bVertices.push_back(vertices.at(k));
- }
- citer cit2 = findCell(i-1, bVertices, true);
- Cell* cell2 = *cit2;
- _cells[i-1].insert(cell2);
-
- }
- }
- i--;
-
- }
-
- printf("Cell complex after repair: %d volumes, %d faces, %d edges and %d vertices.\n",
- getSize(3), getSize(2), getSize(1), getSize(0));
-
- return;
-}
-*/
-
void CellComplex::swapSubdomain(){
for(int i = 0; i < 4; i++){
@@ -1021,7 +736,7 @@ int CellComplex::writeComplexMSH(const std::string &name){
if(vertex->inSubdomain()) physical = 3;
else if(vertex->onDomainBoundary()) physical = 2;
else physical = 1;
- fprintf(fp, "%d %d %d %d %d %d %d\n", vertex->getTag(), 15, 3, 0, 0, physical, vertex->getVertex(0));
+ fprintf(fp, "%d %d %d %d %d %d %d\n", vertex->getTag(), 15, 3, 0, 0, physical, vertex->getVertex(0)->getNum());
}
std::list< std::pair<int, Cell*> > cells;
@@ -1033,7 +748,7 @@ int CellComplex::writeComplexMSH(const std::string &name){
cells = edge->getCells();
for(std::list< std::pair<int, Cell*> >::iterator it = cells.begin(); it != cells.end(); it++){
Cell* cell = (*it).second;
- fprintf(fp, "%d %d %d %d %d %d %d %d\n", cell->getTag(), 1, 3, 0, 0, physical, cell->getVertex(0), cell->getVertex(1));
+ fprintf(fp, "%d %d %d %d %d %d %d %d\n", cell->getTag(), 1, 3, 0, 0, physical, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum());
}
}
@@ -1045,7 +760,7 @@ int CellComplex::writeComplexMSH(const std::string &name){
cells = face->getCells();
for(std::list< std::pair<int, Cell*> >::iterator it = cells.begin(); it != cells.end(); it++){
Cell* cell = (*it).second;
- fprintf(fp, "%d %d %d %d %d %d %d %d %d\n", cell->getTag(), 2, 3, 0, 0, physical, cell->getVertex(0), cell->getVertex(1), cell->getVertex(2));
+ fprintf(fp, "%d %d %d %d %d %d %d %d %d\n", cell->getTag(), 2, 3, 0, 0, physical, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum());
}
}
for(citer cit = firstCell(3); cit != lastCell(3); cit++) {
@@ -1056,7 +771,7 @@ int CellComplex::writeComplexMSH(const std::string &name){
cells = volume->getCells();
for(std::list< std::pair<int, Cell*> >::iterator it = cells.begin(); it != cells.end(); it++){
Cell* cell = (*it).second;
- fprintf(fp, "%d %d %d %d %d %d %d %d %d %d\n", cell->getTag(), 4, 3, 0, 0, physical, cell->getVertex(0), cell->getVertex(1), cell->getVertex(2), cell->getVertex(3));
+ fprintf(fp, "%d %d %d %d %d %d %d %d %d %d\n", cell->getTag(), 4, 3, 0, 0, physical, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum());
}
}
@@ -1072,20 +787,10 @@ void CellComplex::printComplex(int dim){
for (citer cit = firstCell(dim); cit != lastCell(dim); cit++){
Cell* cell = *cit;
for(int i = 0; i < cell->getNumVertices(); i++){
- printf("%d ", cell->getVertex(i));
+ printf("%d ", cell->getVertex(i)->getNum());
}
printf("\n");
}
}
-
-DualCellComplex::DualCellComplex(CellComplex* cellComplex){
-
- for(int i = 0; i < 4; i++){
- _cells[i] = cellComplex->getCells(i);
- }
-
-
-
-}
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index 5eec690c5ac830d2c988829a3f7a4257f16a0800..8961d5cbb12c31382953a7362532c1f6fdee195a 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -15,12 +15,18 @@
#include <queue>
#include "GmshConfig.h"
#include "MElement.h"
+#include "MPoint.h"
+#include "MLine.h"
+#include "MTriangle.h"
+#include "MTetrahedron.h"
#include "GModel.h"
#include "GEntity.h"
#include "GRegion.h"
#include "GFace.h"
#include "GVertex.h"
+
+
// Abstract class representing an elemtary cell of a cell complex.
class Cell
{
@@ -28,11 +34,7 @@ class Cell
// cell dimension
int _dim;
- // maximum number of lower dimensional cells on the boundary of this cell
- int _bdMaxSize;
- // maximum number of higher dimensional cells on the coboundary of this cell
- int _cbdMaxSize;
-
+
// whether this cell belongs to a subdomain
// used in relative homology computation
bool _inSubdomain;
@@ -40,6 +42,7 @@ class Cell
// whether this cell belongs to the boundary of a cell complex
bool _onDomainBoundary;
+ // whether this cell a combinded cell of elemetary cells
bool _combined;
// unique tag for each cell
@@ -63,10 +66,10 @@ class Cell
// get the number of vertices this cell has
- virtual int getNumVertices() const = 0; //{return _vertices.size();}
- virtual int getVertex(int vertex) const = 0; //{return _vertices.at(vertex);}
+ virtual int getNumVertices() const = 0;
+ virtual MVertex* getVertex(int vertex) const = 0; //{return _vertices.at(vertex);}
virtual int getSortedVertex(int vertex) const = 0;
- virtual std::vector<int> getVertexVector() const = 0;
+ virtual std::vector<MVertex*> getVertexVector() const = 0;
// returns 1 or -1 if lower dimensional cell tau is on the boundary of this cell
// otherwise returns 0
@@ -74,11 +77,8 @@ class Cell
virtual int kappa(Cell* tau) const = 0;
// true if this cell has given vertex
- virtual bool hasVertex(int vertex) const =0;
+ virtual bool hasVertex(int vertex) const = 0;
- virtual unsigned int getBdMaxSize() const { return _bdMaxSize; };
- virtual unsigned int getCbdMaxSize() const { return _cbdMaxSize; };
-
virtual int getBoundarySize() { return _boundary.size(); }
virtual int getCoboundarySize() { return _coboundary.size(); }
@@ -146,13 +146,8 @@ class Cell
// print cell vertices
virtual void printCell() const = 0;
- // return the coordinates of this cell (only used on zero-cells if anywhere)
- virtual inline double x() const { return 0; }
- virtual inline double y() const { return 0; }
- virtual inline double z() const { return 0; }
-
virtual int getNumFacets() const { return 0; }
- virtual void getFacetVertices(const int num, std::vector<int> &v) const {};
+ virtual void getFacetVertices(const int num, std::vector<MVertex*> &v) const {};
virtual bool combined() { 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; }
@@ -202,304 +197,159 @@ class Simplex : public Cell
};
// Zero simplex cell type.
-class ZeroSimplex : public Simplex
+class ZeroSimplex : public Simplex, public MPoint
{
private:
- // number of the vertex
- // same as the corresponding vertex in the finite element mesh
- int _v;
- // coordinates of this zero simplex, if needed
- double _x, _y, _z;
+
public:
- ZeroSimplex(int vertex, int tag=0, bool subdomain=false, bool boundary=false, double x=0, double y=0, double z=0) : Simplex(){
- _v = vertex;
- _tag = tag;
- _dim = 0;
- _bdMaxSize = 0;
- _cbdMaxSize = 1000; // big number
- _x = x;
- _y = y;
- _z = z;
- _inSubdomain = subdomain;
- _onDomainBoundary = boundary;
- }
+ ZeroSimplex(std::vector<MVertex*> v, int tag=0, bool subdomain=false, bool boundary=false, int num=0, int part=0)
+ : MPoint(v, num, part){
+ _tag = tag;
+ _inSubdomain = subdomain;
+ _onDomainBoundary = boundary;
+ }
~ZeroSimplex(){}
- //int getDim() const { return 0; }
+ int getDim() const { return 0; }
int getNumVertices() const { return 1; }
- int getVertex(int vertex) const {return _v; }
- int getSortedVertex(int vertex) const {return _v; }
- bool hasVertex(int vertex) const {return (_v == vertex); }
+ MVertex* getVertex(int vertex) const {return _v[0]; }
+ int getSortedVertex(int vertex) const {return _v[0]->getNum(); }
+ bool hasVertex(int vertex) const {return (_v[0]->getNum() == vertex); }
- std::vector<int> getVertexVector() const { return std::vector<int>(1,_v); }
- std::vector<int> getSortedVertexVector() const { return std::vector<int>(1,_v); }
+ std::vector<MVertex*> getVertexVector() const { std::vector<MVertex*> v; v.push_back(_v[0]); return v; }
+ std::vector<int> getSortedVertexVector() const { return std::vector<int>(1,_v[0]->getNum()); }
- inline double x() const { return _x; }
- inline double y() const { return _y; }
- inline double z() const { return _z; }
-
- void printCell() const { printf("Vertices: %d, in subdomain: %d \n", _v, _inSubdomain); }
+ void printCell() const { printf("Vertices: %d, in subdomain: %d \n", _v[0]->getNum(), _inSubdomain); }
};
// One simplex cell type.
-class OneSimplex : public Simplex
+class OneSimplex : public Simplex, public MLine
{
private:
- // numbers of the vertices of this simplex
- // same as the corresponding vertices in the finite element mesh
- int _v[2];
int _vs[2];
public:
- OneSimplex(std::vector<int> vertices, int tag=0, bool subdomain=false, bool boundary=false) : Simplex(){
- _v[0] = vertices.at(0);
- _v[1] = vertices.at(1);
- sort(vertices.begin(), vertices.end());
- _vs[0] = vertices.at(0);
- _vs[1] = vertices.at(1);
- _tag = tag;
- _dim = 1;
- _bdMaxSize = 2;
- _cbdMaxSize = 1000;
- _inSubdomain = subdomain;
- _onDomainBoundary = boundary;
- }
-
-
-
- // constructor for the dummy one simplex
- // used to find another definite one simplex from the cell complex
- // first vertex gives the lower bound from where to look
- OneSimplex(int vertex, int dummy){
- _vs[0] = vertex;
- _vs[1] = dummy;
- }
+
+ OneSimplex(std::vector<MVertex*> v, int tag=0, bool subdomain=false, bool boundary=false, int num=0, int part=0)
+ : MLine(v, num, part){
+ _tag = tag;
+ _inSubdomain = subdomain;
+ _onDomainBoundary = boundary;
+ _vs[0] = v.at(0)->getNum();
+ _vs[1] = v.at(1)->getNum();
+ std::sort(_vs, _vs+2);
+ }
~OneSimplex(){}
- //int getDim() const { return 1; }
+ int getDim() const { return 1; }
int getNumVertices() const { return 2; }
int getNumFacets() const { return 2; }
- int getVertex(int vertex) const {return _v[vertex]; }
+ MVertex* getVertex(int vertex) const {return _v[vertex]; }
int getSortedVertex(int vertex) const {return _vs[vertex]; }
- bool hasVertex(int vertex) const {return (_v[0] == vertex || _v[1] == vertex); }
+ bool hasVertex(int vertex) const {return (_v[0]->getNum() == vertex || _v[1]->getNum() == vertex); }
- std::vector<int> getVertexVector() const {
- return std::vector<int>(_v, _v + sizeof(_v)/sizeof(int) ); }
+ std::vector<MVertex*> getVertexVector() const {
+ return std::vector<MVertex*>(_v, _v + sizeof(_v)/sizeof(MVertex*) ); }
std::vector<int> getSortedVertexVector() const {
return std::vector<int>(_vs, _vs + sizeof(_vs)/sizeof(int) );
}
- void getFacetVertices(const int num, std::vector<int> &v) const {
+ void getFacetVertices(const int num, std::vector<MVertex*> &v) const {
v.resize(1);
v[0] = _v[num];
}
//int kappa(Cell* tau) const;
- void printCell() const { printf("Vertices: %d %d, in subdomain: %d \n", _v[0], _v[1], _inSubdomain); }
+ void printCell() const { printf("Vertices: %d %d, in subdomain: %d \n", _v[0]->getNum(), _v[1]->getNum(), _inSubdomain); }
};
// Two simplex cell type.
-class TwoSimplex : public Simplex
+class TwoSimplex : public Simplex, public MTriangle
{
private:
- // numbers of the vertices of this simplex
- // same as the corresponding vertices in the finite element mesh
- int _v[3];
+
int _vs[3];
- int edges_tri(const int edge, const int vert) const{
- static const int e[3][2] = {
- {0, 1},
- {1, 2},
- {2, 0}
- };
- return e[edge][vert];
- }
-
public:
-
- TwoSimplex(std::vector<int> vertices, int tag=0, bool subdomain=false, bool boundary=false) : Simplex(){
- _v[0] = vertices.at(0);
- _v[1] = vertices.at(1);
- _v[2] = vertices.at(2);
- sort(vertices.begin(), vertices.end());
- _vs[0] = vertices.at(0);
- _vs[1] = vertices.at(1);
- _vs[2] = vertices.at(2);
- _tag = tag;
- _dim = 2;
- _bdMaxSize = 3;
- _cbdMaxSize = 2;
- _inSubdomain = subdomain;
- _onDomainBoundary = boundary;
- }
- // constructor for the dummy one simplex
- TwoSimplex(int vertex, int dummy){
- _v[0] = vertex;
- _v[1] = dummy;
- _v[2] = dummy;
- }
+
+ TwoSimplex(std::vector<MVertex*> v, int tag=0, bool subdomain=false, bool boundary=false, int num=0, int part=0)
+ : MTriangle(v, num, part){
+ _tag = tag;
+ _inSubdomain = subdomain;
+ _onDomainBoundary = boundary;
+ _vs[0] = v.at(0)->getNum();
+ _vs[1] = v.at(1)->getNum();
+ _vs[2] = v.at(2)->getNum();
+ std::sort(_vs, _vs+3);
+ }
~TwoSimplex(){}
- //int getDim() const { return 2; }
+ int getDim() const { return 2; }
int getNumVertices() const { return 3; }
int getNumFacets() const { return 3; }
- int getVertex(int vertex) const {return _v[vertex]; }
+ MVertex* getVertex(int vertex) const {return _v[vertex]; }
int getSortedVertex(int vertex) const {return _vs[vertex]; }
bool hasVertex(int vertex) const {return
- (_v[0] == vertex || _v[1] == vertex || _v[2] == vertex); }
- std::vector<int> getVertexVector() const {
- return std::vector<int>(_v, _v + sizeof(_v)/sizeof(int) ); }
+ (_v[0]->getNum() == vertex || _v[1]->getNum() == vertex || _v[2]->getNum() == vertex); }
+
+ std::vector<MVertex*> getVertexVector() const {
+ return std::vector<MVertex*>(_v, _v + sizeof(_v)/sizeof(MVertex*) ); }
std::vector<int> getSortedVertexVector() const {
return std::vector<int>(_vs, _vs + sizeof(_vs)/sizeof(int) );
}
- void getFacetVertices(const int num, std::vector<int> &v) const {
- v.resize(2);
- v[0] = _v[edges_tri(num, 0)];
- v[1] = _v[edges_tri(num, 1)];
+ void getFacetVertices(const int num, std::vector<MVertex*> &v) const {
+ MTriangle::getEdgeVertices(num, v);
}
- void printCell() const { printf("Vertices: %d %d %d, in subdomain: %d\n", _v[0], _v[1], _v[2], _inSubdomain); }
+ void printCell() const { printf("Vertices: %d %d %d, in subdomain: %d\n", _v[0]->getNum(), _v[1]->getNum(), _v[2]->getNum(), _inSubdomain); }
};
// Three simplex cell type.
-class ThreeSimplex : public Simplex
+class ThreeSimplex : public Simplex, public MTetrahedron
{
private:
- // numbers of the vertices of this simplex
- // same as the corresponding vertices in the finite element mesh
- int _v[4];
int _vs[4];
-
- int faces_tetra(const int face, const int vert) const{
- static const int f[4][3] = {
- {0, 2, 1},
- {0, 1, 3},
- {0, 3, 2},
- {3, 1, 2}
- };
- return f[face][vert];
- }
-
+
public:
-
- ThreeSimplex(std::vector<int> vertices, int tag=0, bool subdomain=false, bool boundary=false) : Simplex(){
- _v[0] = vertices.at(0);
- _v[1] = vertices.at(1);
- _v[2] = vertices.at(2);
- _v[3] = vertices.at(3);
- sort(vertices.begin(), vertices.end());
- _vs[0] = vertices.at(0);
- _vs[1] = vertices.at(1);
- _vs[2] = vertices.at(2);
- _vs[3] = vertices.at(3);
- _tag = tag;
- _dim = 3;
- _bdMaxSize = 4;
- _cbdMaxSize = 0;
- _inSubdomain = subdomain;
- _onDomainBoundary = boundary;
- }
- // constructor for the dummy one simplex
- ThreeSimplex(int vertex, int dummy){
- _v[0] = vertex;
- _v[1] = dummy;
- _v[2] = dummy;
- _v[3] = dummy;
- }
+ ThreeSimplex(std::vector<MVertex*> v, int tag=0, bool subdomain=false, bool boundary=false, int num=0, int part=0)
+ : MTetrahedron(v, num, part){
+ _tag = tag;
+ _inSubdomain = subdomain;
+ _onDomainBoundary = boundary;
+ _vs[0] = v.at(0)->getNum();
+ _vs[1] = v.at(1)->getNum();
+ _vs[2] = v.at(2)->getNum();
+ _vs[3] = v.at(3)->getNum();
+ std::sort(_vs, _vs+4);
+ }
~ThreeSimplex(){}
- //int getDim() const { return 3; }
+ int getDim() const { return 3; }
int getNumVertices() const { return 4; }
int getNumFacets() const { return 4; }
- int getVertex(int vertex) const {return _v[vertex]; }
+ MVertex* getVertex(int vertex) const {return _v[vertex]; }
int getSortedVertex(int vertex) const {return _vs[vertex]; }
bool hasVertex(int vertex) const {return
- (_v[0] == vertex || _v[1] == vertex || _v[2] == vertex || _v[3] == vertex); }
- std::vector<int> getVertexVector() const {
- return std::vector<int>(_v, _v + sizeof(_v)/sizeof(int) ); }
+ (_v[0]->getNum() == vertex || _v[1]->getNum() == vertex || _v[2]->getNum() == vertex || _v[3]->getNum() == vertex); }
+ std::vector<MVertex*> getVertexVector() const {
+ return std::vector<MVertex*>(_v, _v + sizeof(_v)/sizeof(MVertex*) ); }
std::vector<int> getSortedVertexVector() const {
return std::vector<int>(_vs, _vs + sizeof(_vs)/sizeof(int) ); }
- void getFacetVertices(const int num, std::vector<int> &v) const {
- v.resize(3);
- v[0] = _v[faces_tetra(num, 0)];
- v[1] = _v[faces_tetra(num, 1)];
- v[2] = _v[faces_tetra(num, 2)];
+ void getFacetVertices(const int num, std::vector<MVertex*> &v) const {
+ MTetrahedron::getFaceVertices(num, v);
}
- virtual void printCell() const { printf("Vertices: %d %d %d %d, in subdomain: %d \n", _v[0], _v[1], _v[2], _v[3], _inSubdomain); }
+ virtual void printCell() const { printf("Vertices: %d %d %d %d, in subdomain: %d \n", _v[0]->getNum(), _v[1]->getNum(), _v[2]->getNum(), _v[3]->getNum(), _inSubdomain); }
};
-
-class Quadrangle : public Cell
-{
- private:
-
- int _v[4];
- int _vs[4];
-
- int edges_quad(const int edge, const int vert) const{
- static const int e[4][2] = {
- {0, 1},
- {1, 2},
- {2, 3},
- {3, 0}
- };
- return e[edge][vert];
- }
-
- public:
-
- Quadrangle(std::vector<int> vertices, int tag=0, bool subdomain=false, bool boundary=false) : Cell(){
- _v[0] = vertices.at(0);
- _v[1] = vertices.at(1);
- _v[2] = vertices.at(2);
- _v[3] = vertices.at(3);
- sort(vertices.begin(), vertices.end());
- _vs[0] = vertices.at(0);
- _vs[1] = vertices.at(1);
- _vs[2] = vertices.at(2);
- _vs[3] = vertices.at(3);
- _tag = tag;
- _dim = 2;
- _bdMaxSize = 4;
- _cbdMaxSize = 2;
- _inSubdomain = subdomain;
- _onDomainBoundary = boundary;
- }
- ~Quadrangle(){}
-
- //int getDim() const { return 2; }
- int getNumVertices() const { return 4; }
- int getNumFacets() const { return 4; }
- int getVertex(int vertex) const {return _v[vertex]; }
- int getSortedVertex(int vertex) const { return _vs[vertex]; }
- int kappa(Cell* tau) const;
-
- bool hasVertex(int vertex) const {return
- (_v[0] == vertex || _v[1] == vertex || _v[2] == vertex || _v[3] == vertex); }
- std::vector<int> getVertexVector() const {
- return std::vector<int>(_v, _v + sizeof(_v)/sizeof(int) ); }
- std::vector<int> getSortedVertexVector() const {
- return std::vector<int>(_vs, _vs + sizeof(_vs)/sizeof(int) ); }
-
- void getFacetVertices(const int num, std::vector<int> &v) const {
- v.resize(2);
- v[0] = _v[edges_quad(num, 0)];
- v[1] = _v[edges_quad(num, 1)];
- }
-
-};
-
-
// Ordering for the cells.
class Less_Cell{
public:
@@ -519,18 +369,7 @@ class Less_Cell{
if(c1->getSortedVertex(i) < c2->getSortedVertex(i)) return true;
else if (c1->getSortedVertex(i) > c2->getSortedVertex(i)) return false;
}
-
- /*
- std::vector<int> c1v = c1->getVertexVector();
- std::vector<int> c2v = c2->getVertexVector();
- std::sort(c1v.begin(), c1v.end());
- std::sort(c2v.begin(), c2v.end());
-
- for(int i=0; i < c1v.size();i++){
- if(c1v.at(i) < c2v.at(i)) return true;
- else if (c1v.at(i) > c2v.at(i)) return false;
- }
- */
+
return false;
}
};
@@ -563,7 +402,7 @@ class Less_MVertex{
class CombinedCell : public Cell{
private:
- std::vector<int> _v;
+ std::vector<MVertex*> _v;
std::vector<int> _vs;
std::list< std::pair<int, Cell*> > _cells;
@@ -572,18 +411,16 @@ class CombinedCell : public Cell{
CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co=false) : Cell(){
_tag = c1->getTag();
_dim = c1->getDim();
- _bdMaxSize = 1000000;
- _cbdMaxSize = 1000000;
_inSubdomain = c1->inSubdomain();
_onDomainBoundary = c1->onDomainBoundary();
_combined = true;
_v = c1->getVertexVector();
for(int i = 0; i < c2->getNumVertices(); i++){
- if(!this->hasVertex(c2->getVertex(i))) _v.push_back(c2->getVertex(i));
+ if(!this->hasVertex(c2->getVertex(i)->getNum())) _v.push_back(c2->getVertex(i));
}
- _vs = _v;
+ for(int i = 0; i < _v.size(); i++) _vs.push_back(_v.at(i)->getNum());
std::sort(_vs.begin(), _vs.end());
_cells = c1->getCells();
@@ -630,16 +467,16 @@ class CombinedCell : public Cell{
~CombinedCell(){}
int getNumVertices() const { return _v.size(); }
- int getVertex(int vertex) const { return _v.at(vertex); }
+ MVertex* getVertex(int vertex) const { return _v.at(vertex); }
int getSortedVertex(int vertex) const { return _vs.at(vertex); }
- std::vector<int> getVertexVector() const { return _v; }
+ std::vector<MVertex*> getVertexVector() const { return _v; }
int kappa(Cell* tau) const { return 0; }
// true if this cell has given vertex
bool hasVertex(int vertex) const {
for(int i = 0; i < _v.size(); i++){
- if(_v.at(i) == vertex) return true;
+ if(_v.at(i)->getNum() == vertex) return true;
}
return false;
}
@@ -647,7 +484,7 @@ class CombinedCell : public Cell{
virtual void printCell() const {
printf("Vertices: ");
for(int i = 0; i < this->getNumVertices(); i++){
- printf("%d ", this->getVertex(i));
+ printf("%d ", this->getVertex(i)->getNum());
}
printf(", in subdomain: %d\n", _inSubdomain);
}
@@ -748,10 +585,6 @@ class CellComplex
virtual citer firstCell(int dim) {return _cells[dim].begin(); }
virtual citer lastCell(int dim) {return _cells[dim].end(); }
- // find a cell in this cell complex
- //virtual std::set<Cell*, Less_Cell>::iterator findCell(int dim, std::vector<int>& vertices, bool original=false);
- //virtual std::set<Cell*, Less_Cell>::iterator findCell(int dim, int vertex, int dummy=0);
-
// kappa for two cells of this cell complex
// implementation will vary depending on cell type
virtual inline int kappa(Cell* sigma, Cell* tau) const { return sigma->kappa(tau); }
@@ -767,27 +600,17 @@ class CellComplex
virtual bool inSameDomain(Cell* c1, Cell* c2) const { return
( (!c1->inSubdomain() && !c2->inSubdomain()) || (c1->inSubdomain() && c2->inSubdomain()) ); }
- // coreduction of this cell complex
- // removes corection pairs of cells of dimension dim and dim+1
- //virtual int coreduction(int dim);
- //virtual int coreduction();
-
- // stores removed cells
// reduction of this cell complex
// removes reduction pairs of cell of dimension dim and dim-1
virtual int reduction(int dim);
- //virtual int reduction(Cell* generator);
+
// useful functions for (co)reduction of cell complex
virtual void reduceComplex(bool omitHighdim = false);
virtual void coreduceComplex(bool omitlowDim = false);
- //virtual void coreduceComplex(int generatorDim);
- //virtual void coreduceComplex();
-
// queued coreduction presented in Mrozek's paper
- // slower, but produces cleaner result
virtual int coreduction(Cell* generator);
// add every volume, face and edge its missing boundary cells
@@ -827,15 +650,4 @@ class CellComplex
};
-
-class DualCellComplex : public CellComplex
-{
-
- public:
-
- DualCellComplex(CellComplex* cellComplex);
- ~DualCellComplex(){}
-
-};
-
#endif