diff --git a/Geo/Cell.cpp b/Geo/Cell.cpp
index 5cdaaba5cf84065da9abffbef17bde3de94d524d..962a53b3ad9c2239fea66986b1416cfd863c1d34 100755
--- a/Geo/Cell.cpp
+++ b/Geo/Cell.cpp
@@ -6,6 +6,7 @@
// Contributed by Matti Pellikka <matti.pellikka@tut.fi>.
#include "CellComplex.h"
+#include "Cell.h"
#if defined(HAVE_KBIPACK)
@@ -27,12 +28,69 @@ bool Less_Cell::operator()(const Cell* c1, const Cell* c2) const {
return false;
}
+bool Cell::hasVertex(int vertex) const {
+ std::vector<int>::const_iterator it = std::find(_vs.begin(), _vs.end(), vertex);
+ if (it != _vs.end()) return true;
+ else return false;
+}
+
+int Cell::getNumFacets() const {
+ if(getDim() == 0) return 0;
+ else if(getDim() == 1) return 2;
+ else if(getDim() == 2) return _image->getNumEdges();
+ else if(getDim() == 3) return _image->getNumFaces();
+ else return 0;
+}
+void Cell::getFacetVertices(const int num, std::vector<MVertex*> &v) const {
+ if(getDim() == 0) return;
+ else if(getDim() == 1) { v.resize(1); v[0] = getVertex(num); }
+ else if(getDim() == 2) _image->getEdgeVertices(num, v);
+ else if(getDim() == 3) _image->getFaceVertices(num, v);
+ return;
+}
+
+
+int Cell::getFacetOri(std::vector<MVertex*> &v) {
+ if(getDim() == 0) return 0;
+ else if(getDim() == 1){
+ if(v.size() != 1) return 0;
+ else if(v[0] == getVertex(0)) return -1;
+ else if(v[0] == getVertex(1)) return 1;
+ else return 0;
+ }
+ else if(getDim() == 2){
+ if(v.size() != 2) return 0;
+ MEdge facet = MEdge(v[0], v[1]);
+ int ithFacet = 0;
+ int sign = 0;
+ _image->getEdgeInfo(facet, ithFacet, sign);
+ return sign;
+ }
+ else if(getDim() == 3){
+ if(v.size() != 3) return 0;
+ MFace facet = MFace(v);
+ int ithFacet = 0;
+ int sign = 0;
+ int rot = 0;
+ _image->getFaceInfo(facet, ithFacet, sign, rot);
+ return sign;
+ }
+ else return 0;
+}
+
+void Cell::printCell() {
+ printf("%d-cell %d: \n" , getDim(), getNum());
+ printf("Vertices: ");
+ for(int i = 0; i < this->getNumVertices(); i++){
+ printf("%d ", this->getSortedVertex(i));
+ }
+ printf(", in subdomain: %d\n", _inSubdomain);
+ printf("Combined: %d \n" , isCombined() );
+};
void Cell::restoreCell(){
_boundary = _obd;
_coboundary = _ocbd;
- _bdSize = _boundary.size();
- _cbdSize = _coboundary.size();
_combined = false;
_index = 0;
_immune = false;
@@ -58,12 +116,11 @@ std::list< Cell* > Cell::getCoboundary() {
bool Cell::addBoundaryCell(int orientation, Cell* cell) {
- _bdSize++;
biter it = _boundary.find(cell);
if(it != _boundary.end()){
(*it).second = (*it).second + orientation;
if((*it).second == 0) {
- _boundary.erase(it); _bdSize--;
+ _boundary.erase(it);
(*it).first->removeCoboundaryCell(this,false);
return false;
}
@@ -74,12 +131,11 @@ bool Cell::addBoundaryCell(int orientation, Cell* cell) {
}
bool Cell::addCoboundaryCell(int orientation, Cell* cell) {
- _cbdSize++;
biter it = _coboundary.find(cell);
if(it != _coboundary.end()){
(*it).second = (*it).second + orientation;
if((*it).second == 0) {
- _coboundary.erase(it); _cbdSize--;
+ _coboundary.erase(it);
(*it).first->removeBoundaryCell(this,false);
return false;
}
@@ -94,7 +150,6 @@ int Cell::removeBoundaryCell(Cell* cell, bool other) {
if(it != _boundary.end()){
_boundary.erase(it);
if(other) (*it).first->removeCoboundaryCell(this, false);
- _bdSize--;
return (*it).second;
}
@@ -105,7 +160,6 @@ int Cell::removeCoboundaryCell(Cell* cell, bool other) {
if(it != _coboundary.end()){
_coboundary.erase(it);
if(other) (*it).first->removeBoundaryCell(this, false);
- _cbdSize--;
return (*it).second;
}
return 0;
@@ -189,26 +243,27 @@ CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell() {
}
_index = c1->getIndex();
- _tag = c1->getTag();
+ //_tag = c1->getTag();
_dim = c1->getDim();
_num = c1->getNum();
_inSubdomain = c1->inSubdomain();
_onDomainBoundary = c1->onDomainBoundary();
_combined = true;
+ _image = NULL;
_v.reserve(c1->getNumVertices() + c2->getNumVertices());
_vs.reserve(c1->getNumVertices() + c2->getNumVertices());
-
- _v = c1->getVertexVector();
+
+ for(int i = 0; i < c1->getNumVertices(); i++) _v.push_back(c1->getVertex(i));
for(int i = 0; i < c2->getNumVertices(); i++){
- if(!this->hasVertex(c2->getVertex(i)->getNum())) _v.push_back(c2->getVertex(i));
+ if(!this->hasVertex(c2->getVertex(i)->getNum())) _v.push_back(c2->getVertex(i));
}
-
+
// sorted vertices
for(unsigned int i = 0; i < _v.size(); i++) _vs.push_back(_v.at(i)->getNum());
std::sort(_vs.begin(), _vs.end());
-
+
// cells
_cells = c1->getCells();
std::list< std::pair<int, Cell*> > c2Cells = c2->getCells();
@@ -216,7 +271,7 @@ CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell() {
if(!orMatch) (*it).first = -1*(*it).first;
_cells.push_back(*it);
}
-
+
// boundary cells
std::map< Cell*, int, Less_Cell > c1Boundary = c1->getOrientedBoundary();
std::map< Cell*, int, Less_Cell > c2Boundary = c2->getOrientedBoundary();
@@ -242,7 +297,7 @@ CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell() {
if(this->addBoundaryCell(ori, cell)) cell->addCoboundaryCell(ori, this);
}
}
-
+
// coboundary cells
std::map<Cell*, int, Less_Cell > c1Coboundary = c1->getOrientedCoboundary();
std::map<Cell*, int, Less_Cell > c2Coboundary = c2->getOrientedCoboundary();
@@ -268,7 +323,7 @@ CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell() {
if(this->addCoboundaryCell(ori, cell)) cell->addBoundaryCell(ori, this);
}
}
-
+
}
@@ -280,61 +335,4 @@ CombinedCell::~CombinedCell(){
}
}
-bool CombinedCell::hasVertex(int vertex) const {
- /*std::vector<int>::const_iterator it = std::find(_vs.begin(), _vs.end(), vertex);
- if (it != _vs.end()) return true;
- else return false;*/
- for(unsigned int i = 0; i < _v.size(); i++){
- if(_v.at(i)->getNum() == vertex) return true;
- }
- return false;
-}
-
-void CombinedCell::printCell() const {
- printf("Cell dimension: %d, ", getDim() );
- printf("Vertices: ");
- for(int i = 0; i < this->getNumVertices(); i++){
- printf("%d ", this->getSortedVertex(i));
- }
- printf(", in subdomain: %d\n", _inSubdomain);
-}
-
-/*
- int Cell::incidence(Cell* tau) const{
- for(int i=0; i < tau->getNumVertices(); i++){
- if( !(this->hasVertex(tau->getVertex(i)->getNum())) ) return 0;
- }
-
- if(this->getDim() - tau->getDim() != 1) return 0;
-
- int value=1;
-
- for(int i = 0; i < this->getNumFacets(); i++){
- std::vector<MVertex*> vTau = tau->getVertexVector();
- std::vector<MVertex*> v;
- this->getFacetVertices(i, v);
- value = -1;
-
- if(v.size() != vTau.size()) printf("Error: invalid facet!");
-
- do {
- value = value*-1;
- if(v == vTau) return value;
- }
- while (std::next_permutation(vTau.begin(), vTau.end()) );
-
- vTau = tau->getVertexVector();
- value = -1;
- do {
- value = value*-1;
- if(v == vTau) return value;
- }
- while (std::prev_permutation(vTau.begin(), vTau.end()) );
-
-
- }
-
- return 0;
-}*/
-
#endif
diff --git a/Geo/Cell.h b/Geo/Cell.h
new file mode 100644
index 0000000000000000000000000000000000000000..3fb2439e8300f8ffb4c4b76ab0f65a5634216d6a
--- /dev/null
+++ b/Geo/Cell.h
@@ -0,0 +1,247 @@
+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+//
+// Contributed by Matti Pellikka <matti.pellikka@tut.fi>.
+
+#ifndef _CELL_H_
+#define _CELL_H_
+
+#include "GmshConfig.h"
+
+#if defined(HAVE_KBIPACK)
+
+#include <stdio.h>
+#include <string>
+#include <algorithm>
+#include <set>
+#include <queue>
+#include "CellComplex.h"
+#include "MElement.h"
+#include "MPoint.h"
+#include "MLine.h"
+#include "MTriangle.h"
+#include "MQuadrangle.h"
+#include "MHexahedron.h"
+#include "MPrism.h"
+#include "MPyramid.h"
+#include "MTetrahedron.h"
+#include "GModel.h"
+
+class Cell;
+
+// Ordering of the cells
+class Less_Cell{
+ public:
+ bool operator()(const Cell* c1, const Cell* c2) const;
+};
+
+// Class representing an elementary cell of a cell complex.
+class Cell
+{
+ protected:
+
+ // cell dimension
+ int _dim;
+
+ // whether this cell belongs to a subdomain, immutable
+ // used in relative homology computation
+ bool _inSubdomain;
+
+ // whether this cell belongs to the boundary of a cell complex
+ bool _onDomainBoundary;
+
+ // whether this cell a combinded cell of elemetary cells
+ bool _combined;
+
+ // mutable index for each cell (used to create boundary operator matrices)
+ int _index;
+
+ // for some algorithms to omit this cell
+ bool _immune;
+
+ // mutable list of cells on the boundary and on the coboundary of this cell
+ std::map< Cell*, int, Less_Cell > _boundary;
+ std::map< Cell*, int, Less_Cell > _coboundary;
+
+ // immutable original boundary and coboundary before the reduction of the cell complex
+ std::map<Cell*, int, Less_Cell > _obd;
+ std::map<Cell*, int, Less_Cell > _ocbd;
+
+ // The mesh element that is the image of this cell
+ MElement* _image;
+ // Whether to delete the mesh element when done (created for homology computation only)
+ bool _deleteImage;
+
+ // sorted vertices of this cell (used for ordering of the cells)
+ std::vector<int> _vs;
+
+ public:
+ Cell(){
+ _combined = false;
+ _index = 0;
+ _immune = false;
+ _image = NULL;
+ _deleteImage = false;
+ }
+
+ Cell(MElement* image, bool subdomain, bool boundary){
+ _combined = false;
+ _index = 0;
+ _immune = false;
+ _image = image;
+ _onDomainBoundary = boundary;
+ _inSubdomain = subdomain;
+ _dim = image->getDim();
+ _deleteImage = false;
+ for(int i = 0; i < image->getNumVertices(); i++) _vs.push_back(image->getVertex(i)->getNum());
+ std::sort(_vs.begin(), _vs.end());
+ }
+ virtual ~Cell() {if(_deleteImage) delete _image; }
+
+ virtual MElement* getImageMElement() const { return _image; };
+
+ virtual int getDim() const { return _dim; };
+ virtual int getIndex() const { return _index; };
+ virtual void setIndex(int index) { _index = index; };
+ virtual int getNum() { return _image->getNum(); }
+ virtual int getType() { return _image->getType(); }
+ virtual int getTypeForMSH() { return _image->getTypeForMSH(); }
+ virtual int getPartition() { return _image->getPartition(); }
+ virtual void setImmune(bool immune) { _immune = immune; };
+ virtual bool getImmune() const { return _immune; };
+ virtual void setDeleteImage(bool deleteImage) { _deleteImage = deleteImage; };
+ virtual bool getDeleteImage() const { return _deleteImage; };
+
+ // get the number of vertices this cell has
+ virtual int getNumVertices() const { return _image->getNumVertices(); }
+ virtual MVertex* getVertex(int vertex) const { return _image->getVertex(vertex); }
+ virtual int getSortedVertex(int vertex) const { return _vs.at(vertex); }
+
+ // restores the cell information to its original state before reduction
+ virtual void restoreCell();
+
+ // true if this cell has given vertex
+ virtual bool hasVertex(int vertex) const;
+
+ // (co)boundary cell iterator
+ typedef std::map<Cell*, int, Less_Cell>::iterator biter;
+
+ virtual int getBoundarySize() { return _boundary.size(); }
+ virtual int getCoboundarySize() { return _coboundary.size(); }
+
+ // get the cell boundary
+ virtual std::map<Cell*, int, Less_Cell > getOrientedBoundary() { return _boundary; }
+ virtual std::list< Cell* > getBoundary();
+ virtual std::map<Cell*, int, Less_Cell > getOrientedCoboundary() { return _coboundary; }
+ virtual std::list< Cell* > getCoboundary();
+
+ // add (co)boundary cell
+ virtual bool addBoundaryCell(int orientation, Cell* cell);
+ virtual bool addCoboundaryCell(int orientation, Cell* cell);
+
+ // remove (co)boundary cell
+ virtual int removeBoundaryCell(Cell* cell, bool other=true);
+ virtual int removeCoboundaryCell(Cell* cell, bool other=true);
+
+ // true if has given cell on (original) (co)boundary
+ virtual bool hasBoundary(Cell* cell, bool org=false);
+ virtual bool hasCoboundary(Cell* cell, bool org=false);
+
+ virtual void clearBoundary() { _boundary.clear(); }
+ virtual void clearCoboundary() { _coboundary.clear(); }
+
+ // algebraic dual of the cell
+ virtual void makeDualCell();
+
+ // print cell info
+ virtual void printCell();
+ virtual void printBoundary();
+ virtual void printCoboundary();
+
+ // original (co)boundary
+ virtual void addOrgBdCell(int orientation, Cell* cell) { _obd.insert( std::make_pair(cell, orientation ) ); };
+ virtual void addOrgCbdCell(int orientation, Cell* cell) { _ocbd.insert( std::make_pair(cell, orientation ) ); };
+ virtual std::map<Cell*, int, Less_Cell > getOrgBd() { return _obd; }
+ virtual std::map<Cell*, int, Less_Cell > getOrgCbd() { return _ocbd; }
+ virtual void printOrgBd();
+ virtual void printOrgCbd();
+
+ virtual bool inSubdomain() const { return _inSubdomain; }
+ //virtual void setInSubdomain(bool subdomain) { _inSubdomain = subdomain; }
+
+ virtual bool onDomainBoundary() const { return _onDomainBoundary; }
+ virtual void setOnDomainBoundary(bool domainboundary) { _onDomainBoundary = domainboundary; }
+
+ // get the number of facets of this cell
+ virtual int getNumFacets() const;
+ // get the vertices on a facet of this cell
+ virtual void getFacetVertices(const int num, std::vector<MVertex*> &v) const;
+
+ // get boundary cell orientation
+ virtual int getFacetOri(Cell* cell) {
+ std::vector<MVertex*> v;
+ for(int i = 0; i < cell->getNumVertices(); i++) v.push_back(cell->getVertex(i));
+ return getFacetOri(v);
+ }
+ virtual int getFacetOri(std::vector<MVertex*> &v);
+
+ // tools for combined cells
+ virtual bool isCombined() { return _combined; }
+ virtual std::list< std::pair<int, Cell*> > getCells() {
+ std::list< std::pair<int, Cell*> >cells;
+ cells.push_back( std::make_pair(1, this));
+ return cells;
+ }
+ virtual int getNumCells() {return 1;}
+
+ // equivalence
+ bool operator==(const Cell& c2) const {
+ if(this->getNumVertices() != c2.getNumVertices()){
+ return false;
+ }
+ for(int i=0; i < this->getNumVertices();i++){
+ if(this->getSortedVertex(i) != c2.getSortedVertex(i)){
+ return false;
+ }
+ }
+ return true;
+ }
+};
+
+// A cell that is a combination of cells of same dimension
+class CombinedCell : public Cell{
+
+ private:
+ // vertices
+ std::vector<MVertex*> _v;
+ // sorted list of all vertices
+ std::vector<int> _vs;
+ // list of cells this cell is a combination of
+ std::list< std::pair<int, Cell*> > _cells;
+ int _num;
+
+ public:
+
+ CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co=false);
+ ~CombinedCell();
+
+ int getNum() { return _num; }
+ int getType() { return 0; }
+ int getTypeForMSH() { return 0; }
+ int getPartition() { return 0; }
+
+ int getNumVertices() const { return _v.size(); }
+ MVertex* getVertex(int vertex) const { return _v.at(vertex); }
+ int getSortedVertex(int vertex) const { return _vs.at(vertex); }
+
+ std::list< std::pair<int, Cell*> > getCells() { return _cells; }
+ int getNumCells() {return _cells.size();}
+
+};
+
+
+#endif
+
+#endif
diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 8b49f25cee53fdaa0bf5285a8e4d8e4188807a4b..cb2b779a5a8bb1914192f955c3044561acfe870f 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -26,19 +26,9 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su
insert_cells(true, true);
insert_cells(false, true);
insert_cells(false, false);
-
- // find vertices in domain
- find_vertices(domain, subdomain);
- //int tag = 1;
for(int i = 0; i < 4; i++){
- /*
- for(citer cit = firstCell(i); cit != lastCell(i); cit++){
- Cell* cell = *cit;
- cell->setTag(tag);
- tag++;
- }*/
_ocells[i] = _cells[i];
_betti[i] = 0;
if(getSize(i) > _dim) _dim = i;
@@ -47,27 +37,6 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su
}
-void CellComplex::find_vertices(std::vector<GEntity*>& domain, std::vector<GEntity*>& subdomain){
- // 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++){
- for(int k=0; k < domain.at(j)->getMeshElement(i)->getNumVertices(); k++){
- MVertex* vertex = domain.at(j)->getMeshElement(i)->getVertex(k);
- _domainVertices.insert(vertex);
- }
- }
- }
-
- 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++){
- MVertex* vertex = subdomain.at(j)->getMeshElement(i)->getVertex(k);
- _domainVertices.insert(vertex);
- }
- }
- }
-}
-
void CellComplex::find_boundary(std::vector<GEntity*>& domain, std::vector<GEntity*>& subdomain){
// determine mesh entities on boundary of the domain
@@ -140,7 +109,6 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){
std::vector<MVertex*> vertices;
std::pair<citer, bool> insertInfo;
- int tag = 1;
// add highest dimensional cells
for(unsigned int j=0; j < domain.size(); j++) {
@@ -161,15 +129,12 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){
type == MSH_TRI_15 || type == MSH_TRI_15I || type == MSH_TRI_21 || type == MSH_LIN_4 ||
type == MSH_LIN_5 || type == MSH_LIN_6 || type == MSH_TET_20 || type == MSH_TET_35 || type == MSH_TET_56
|| type == MSH_TET_34 || type == MSH_TET_52 ){
- if(dim == 3) cell = new ThreeSimplex(vertices, tag, subdomain, boundary);
- else if(dim == 2) cell = new TwoSimplex(vertices, tag, subdomain, boundary);
- else if(dim == 1) cell = new OneSimplex(vertices, tag, subdomain, boundary);
- else cell = new ZeroSimplex(vertices, tag, subdomain, boundary);
+ cell = new Cell(domain.at(j)->getMeshElement(i), subdomain, boundary);
}
else if(type == MSH_QUA_4 || type == MSH_QUA_8 || type == MSH_QUA_9){
- cell = new CQuadrangle(vertices, tag, subdomain, boundary);
+ cell = new Cell(domain.at(j)->getMeshElement(i), subdomain, boundary);
_simplicial = false;
- }/* FIXME: reduction doesn't work for these
+ }/* FIXME: reduction doesn't work for these (but should).
else if(type == MSH_HEX_8 || type == MSH_HEX_27 || type == MSH_HEX_20){
cell = new CHexahedron(vertices, tag, subdomain, boundary);
_simplicial = false;
@@ -182,10 +147,7 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){
cell = new CPyramid(vertices, tag, subdomain, boundary);
_simplicial = false;
}*/
- else Msg::Error("Error: mesh element %d not implemented yet! \n", type);
- tag++;
- //if(domain.at(j)->getMeshElement(i)->getNum() == 8196) cell->setImmune(true);
- //else
+ else Msg::Error("Error: mesh element %d not implemented yet! \n", type);
cell->setImmune(false);
insertInfo = _cells[dim].insert(cell);
if(!insertInfo.second) delete cell;
@@ -193,28 +155,40 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){
}
// add lower dimensional cells recursively
+ MElementFactory factory;
for (int dim = 3; dim > 0; dim--){
for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
Cell* cell = *cit;
std::vector<MVertex*> vertices;
for(int i = 0; i < cell->getNumFacets(); i++){
cell->getFacetVertices(i, vertices);
- Cell* newCell;
+ int type = cell->getTypeForMSH();
+
+ int newtype = 0;
+ //FIXME: add missing boundary cell type relations
if(dim == 3){
- if(vertices.size() == 3) newCell = new TwoSimplex(vertices, tag, subdomain, boundary);
- else if(vertices.size() == 4) newCell = new CQuadrangle(vertices, tag, subdomain, boundary);
- else Msg::Debug("Error: invalid face! \n");
+ if(type == MSH_TET_4) newtype = MSH_TRI_3;
+ else if(type == MSH_TET_10) newtype = MSH_TRI_6;
+ }
+ else if(dim == 2){
+ if(type == MSH_TRI_3 || type == MSH_QUA_4) newtype = MSH_LIN_2;
+ else if(type == MSH_TRI_6 || type == MSH_QUA_8) newtype = MSH_LIN_3;
}
- else if(dim == 2) newCell = new OneSimplex(vertices, tag, subdomain, boundary);
- else if(dim == 1) newCell = new ZeroSimplex(vertices, tag, subdomain, boundary);
+ else if(dim == 1){
+ if(type == MSH_LIN_2 || type == MSH_LIN_3 || type == MSH_LIN_4 ||
+ type == MSH_LIN_5 || type == MSH_LIN_6) newtype = MSH_PNT;
+ }
+ if(newtype == 0) Msg::Error("Error: mesh element %d not implemented yet! \n", type);
+
+ MElement* element = factory.create(newtype, vertices, 0, cell->getPartition());
+ Cell* newCell = new Cell(element, subdomain, boundary);
newCell->setImmune(cell->getImmune());
- tag++;
+ newCell->setDeleteImage(true);
insertInfo = _cells[dim-1].insert(newCell);
if(!insertInfo.second){
- delete newCell;
+ delete newCell;
Cell* oldCell = *(insertInfo.first);
if(!subdomain && !boundary){
- //int ori = cell->kappa(oldCell);
int ori = cell->getFacetOri(oldCell);
oldCell->addCoboundaryCell( ori, cell );
oldCell->addOrgCbdCell( ori, cell );
@@ -223,7 +197,6 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){
}
}
else if(!subdomain && !boundary) {
- //int ori = cell->kappa(newCell);
int ori = cell->getFacetOri(vertices);
cell->addBoundaryCell( ori, newCell );
cell->addOrgBdCell( ori, newCell );
@@ -237,36 +210,35 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){
}
CellComplex::~CellComplex(){
- for(int i = 0; i < 4; i++){
- _cells[i].clear();
-
- for(citer cit = _ocells[i].begin(); cit != _ocells[i].end(); cit++){
- Cell* cell = *cit;
- delete cell;
- }
- _ocells[i].clear();
-
- }
- //emptyTrash();
- }
+ for(int i = 0; i < 4; i++){
+ _cells[i].clear();
+
+ for(citer cit = _ocells[i].begin(); cit != _ocells[i].end(); cit++){
+ Cell* cell = *cit;
+ delete cell;
+ }
+ _ocells[i].clear();
+
+ }
+}
- /*
- CellComplex::CellComplex(CellComplex* cellComplex){
-
- _domain = cellComplex->getDomain();
- _subdomain = cellComplex->getSubdomain;
- _boundary = cellComplex->getBoundary;
+/*
+ CellComplex::CellComplex(CellComplex* cellComplex){
+ *
+ _domain = cellComplex->getDomain();
+ _subdomain = cellComplex->getSubdomain;
+ _boundary = cellComplex->getBoundary;
_domainVertices = cellComplex->getDomainVertices;
-
- for(int i = 0; i < 4; i++){
- _betti[i] = cellComplex->getBetti(i);
- _cells[i] = ocells[i];
- _ocells[i] = ocells[i];
- }
-
- _dim = cellComplex->getDim();
+ *
+ for(int i = 0; i < 4; i++){
+ _betti[i] = cellComplex->getBetti(i);
+ _cells[i] = ocells[i];
+ _ocells[i] = ocells[i];
+ }
+ *
+ _dim = cellComplex->getDim();
}
- */
+ */
void CellComplex::removeCell(Cell* cell, bool other){
@@ -693,12 +665,13 @@ int CellComplex::combine(int dim){
std::map<Cell*, int, Less_Cell> cbd_c;
std::list<Cell*> bd_c;
int count = 0;
-
+
for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
Cell* cell = *cit;
bd_c = cell->getBoundary();
enqueueCells(bd_c, Q, Qset);
+
while(Q.size() != 0){
Cell* s = Q.front();
@@ -725,10 +698,11 @@ int CellComplex::combine(int dim){
enqueueCells(bd_c, Q, Qset);
bd_c = c2->getBoundary();
enqueueCells(bd_c, Q, Qset);
-
+
CombinedCell* newCell = new CombinedCell(c1, c2, (or1 != or2) );
removeCell(c1);
removeCell(c2);
+
std::pair<citer, bool> insertInfo = _cells[dim].insert(newCell);
if(!insertInfo.second) Msg::Debug("Warning: Combined cell not inserted! \n");
cit = firstCell(dim);
@@ -827,152 +801,6 @@ void CellComplex::swapSubdomain(){
}*/
-int CellComplex::writeComplexMSH(const std::string &name){
-
- FILE *fp = fopen(name.c_str(), "w");
-
- if(!fp){
- Msg::Error("Unable to open file '%s'", name.c_str());
- Msg::Debug("Unable to open file.");
- return 0;
- }
-
- fprintf(fp, "$MeshFormat\n2.1 0 8\n$EndMeshFormat\n");
- fprintf(fp, "$Nodes\n");
- fprintf(fp, "%d\n", (int)_domainVertices.size());
-
- for(std::set<MVertex*, Less_MVertex>::iterator vit = _domainVertices.begin(); vit != _domainVertices.end(); vit++){
- MVertex* vertex = *vit;
- fprintf(fp, "%d %.16g %.16g %.16g\n", vertex->getNum(), vertex->x(), vertex->y(), vertex->z() );
- }
-
-
- fprintf(fp, "$EndNodes\n");
- fprintf(fp, "$Elements\n");
-
- int count = 0;
- for(int i = 0; i < 4; i++){
- for(citer cit = firstCell(i); cit != lastCell(i); cit++){
- Cell* cell = *cit;
- count = count + cell->getNumCells();
- }
- }
- for(int i = 0; i < _store.size(); i++){
- count = count + _store.at(i).size();
- }
-
- fprintf(fp, "%d\n", count);
-
- int type = -1;
- int physical = 0;
- int elementary = 0;
- int partition = 0;
-
- for(citer cit = firstCell(0); cit != lastCell(0); cit++) {
- Cell* vertex = *cit;
- if(vertex->inSubdomain()) partition = 3;
- else if(vertex->onDomainBoundary()) partition = 2;
- else partition = 1;
- type = vertex->getTypeForMSH();
- //vertex->writeMSH(fp, 2.1, false, 0, elementary, physical);
- fprintf(fp, "%d %d %d %d %d %d %d\n", vertex->getNum(), type, 3, physical, elementary, partition, vertex->getVertex(0)->getNum());
- }
-
- std::list< std::pair<int, Cell*> > cells;
- for(citer cit = firstCell(1); cit != lastCell(1); cit++) {
- Cell* edge = *cit;
- if(edge->inSubdomain()) partition = 3;
- else if(edge->onDomainBoundary()) partition = 2;
- else partition = 1;
- cells = edge->getCells();
- for(std::list< std::pair<int, Cell*> >::iterator it = cells.begin(); it != cells.end(); it++){
- Cell* cell = (*it).second;
- type = cell->getTypeForMSH();
- //cell->writeMSH(fp, 2.1, false, 0, elementary, physical);
- fprintf(fp, "%d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum());
- }
- }
-
- for(citer cit = firstCell(2); cit != lastCell(2); cit++) {
- Cell* face = *cit;
- if(face->inSubdomain()) partition = 3;
- else if(face->onDomainBoundary()) partition = 2;
- else partition = 1;
- cells = face->getCells();
- for(std::list< std::pair<int, Cell*> >::iterator it = cells.begin(); it != cells.end(); it++){
- Cell* cell = (*it).second;
- type = cell->getTypeForMSH();
- if(cell->getNumVertices() == 3) fprintf(fp, "%d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum());
- else if (cell->getNumVertices() == 4) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum());
- }
- }
- for(citer cit = firstCell(3); cit != lastCell(3); cit++) {
- Cell* volume = *cit;
- if(volume->inSubdomain()) partition = 3;
- else if(volume->onDomainBoundary()) partition = 2;
- else partition = 1;
- cells = volume->getCells();
- for(std::list< std::pair<int, Cell*> >::iterator it = cells.begin(); it != cells.end(); it++){
- Cell* cell = (*it).second;
- type = cell->getTypeForMSH();
- if(cell->getNumVertices() == 4) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum());
- if(cell->getNumVertices() == 8) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(), cell->getVertex(4)->getNum(), cell->getVertex(5)->getNum(), cell->getVertex(6)->getNum(), cell->getVertex(7)->getNum() );
- if(cell->getNumVertices() == 6) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(),cell->getVertex(4)->getNum(), cell->getVertex(5)->getNum());
- if(cell->getNumVertices() == 5) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(), cell->getVertex(4)->getNum());
- }
- }
-
- for(int i = 0; i < _store.size(); i++){
- for(citer cit = _store.at(i).begin(); cit != _store.at(i).end(); cit++){
- Cell* cell = *cit;
- type = cell->getTypeForMSH();
- if(cell->inSubdomain()) partition = 3;
- else if(cell->onDomainBoundary()) partition = 2;
- else partition = 1;
- if(cell->getDim() == 0) fprintf(fp, "%d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum());
- if(cell->getDim() == 1) fprintf(fp, "%d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum());
- if(cell->getDim() == 2){
- if(cell->getNumVertices() == 3) fprintf(fp, "%d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum());
- else if (cell->getNumVertices() == 4) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum());
- }
- if(cell->getDim() == 3){
- if(cell->getNumVertices() == 4) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum());
- if(cell->getNumVertices() == 8) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(), cell->getVertex(4)->getNum(), cell->getVertex(5)->getNum(), cell->getVertex(6)->getNum(), cell->getVertex(7)->getNum() );
- if(cell->getNumVertices() == 6) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(),cell->getVertex(4)->getNum(), cell->getVertex(5)->getNum());
- if(cell->getNumVertices() == 5) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(), cell->getVertex(4)->getNum());
- }
- }
- }
-
- /*
- for(citer cit = _ecells.begin(); cit != _ecells.end(); cit++){
- Cell* cell = *cit;
- type = cell->getTypeForMSH();
- if(cell->inSubdomain()) partition = 3;
- else if(cell->onDomainBoundary()) partition = 2;
- else partition = 1;
- if(cell->getDim() == 0) fprintf(fp, "%d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum());
- if(cell->getDim() == 1) fprintf(fp, "%d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum());
- if(cell->getDim() == 2){
- if(cell->getNumVertices() == 3) fprintf(fp, "%d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum());
- else if (cell->getNumVertices() == 4) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum());
- }
- if(cell->getDim() == 3){
- if(cell->getNumVertices() == 4) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum());
- if(cell->getNumVertices() == 8) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(), cell->getVertex(4)->getNum(), cell->getVertex(5)->getNum(), cell->getVertex(6)->getNum(), cell->getVertex(7)->getNum() );
- if(cell->getNumVertices() == 6) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(),cell->getVertex(4)->getNum(), cell->getVertex(5)->getNum());
- if(cell->getNumVertices() == 5) fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d\n", cell->getNum(), type, 3, physical, elementary, partition, cell->getVertex(0)->getNum(), cell->getVertex(1)->getNum(), cell->getVertex(2)->getNum(), cell->getVertex(3)->getNum(), cell->getVertex(4)->getNum());
- }
- }*/
-
- fprintf(fp, "$EndElements\n");
-
- fclose(fp);
-
- return 1;
-}
-
-
void CellComplex::printComplex(int dim){
for (citer cit = firstCell(dim); cit != lastCell(dim); cit++){
Cell* cell = *cit;
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index ff2bd3b08cc81892497804f68e15cb19341dd649..49c9ffc07270df86eb1c188ed93d9d21e6157c3d 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -17,15 +17,8 @@
#include <algorithm>
#include <set>
#include <queue>
+#include "Cell.h"
#include "MElement.h"
-#include "MPoint.h"
-#include "MLine.h"
-#include "MTriangle.h"
-#include "MQuadrangle.h"
-#include "MHexahedron.h"
-#include "MPrism.h"
-#include "MPyramid.h"
-#include "MTetrahedron.h"
#include "GModel.h"
#include "GEntity.h"
#include "GRegion.h"
@@ -34,648 +27,6 @@
class Cell;
-// Ordering of the cells
-class Less_Cell{
- public:
- bool operator()(const Cell* c1, const Cell* c2) const;
-};
-
-// Abstract class representing an elementary cell of a cell complex.
-class Cell
-{
- protected:
-
- // cell dimension
- int _dim;
-
- // whether this cell belongs to a subdomain, immutable
- // used in relative homology computation
- bool _inSubdomain;
-
- // whether this cell belongs to the boundary of a cell complex
- bool _onDomainBoundary;
-
- // whether this cell a combinded cell of elemetary cells
- bool _combined;
-
- // unused
- int _tag;
-
- // mutable index for each cell
- int _index;
-
- // for some algorithms to omit this cell
- bool _immune;
-
- // mutable list of cells on the boundary and on the coboundary of this cell
- std::map< Cell*, int, Less_Cell > _boundary;
- std::map< Cell*, int, Less_Cell > _coboundary;
- int _bdSize;
- int _cbdSize;
-
- // immutable original boundary and coboundary before the reduction of the cell complex
- std::map<Cell*, int, Less_Cell > _obd;
- std::map<Cell*, int, Less_Cell > _ocbd;
-
- public:
- Cell(){
- _bdSize = 0;
- _cbdSize = 0;
- _combined = false;
- _index = 0;
- _immune = false;
- }
- virtual ~Cell(){}
-
- virtual int getDim() const { return _dim; };
- virtual int getTag() const { return _tag; };
- virtual void setTag(int tag) { _tag = tag; };
- virtual int getIndex() const { return _index; };
- virtual void setIndex(int index) { _index = index; };
- virtual int getNum() { return -1; }
- virtual int getType() { return -1; }
- virtual int getTypeForMSH() { return -1; }
- virtual int getPartition() { return -1; }
- virtual void setPartition(int num){}
- virtual void setImmune(bool immune) { _immune = immune; };
- virtual bool getImmune() const { return _immune; };
-
- // get the number of vertices this cell has
- virtual int getNumVertices() const = 0;
- virtual MVertex* getVertex(int vertex) const = 0; //{return _vertices.at(vertex);}
- virtual int getSortedVertex(int vertex) const = 0;
- virtual std::vector<MVertex*> getVertexVector() const = 0;
-
- // restores the cell information to its original state before reduction
- virtual void restoreCell();
-
- // true if this cell has given vertex
- virtual bool hasVertex(int vertex) const = 0;
-
- // (co)boundary cell iterator
- typedef std::map<Cell*, int, Less_Cell>::iterator biter;
-
- virtual int getBoundarySize() { return _boundary.size(); }
- virtual int getCoboundarySize() { return _coboundary.size(); }
-
- // get the cell boundary
- virtual std::map<Cell*, int, Less_Cell > getOrientedBoundary() { return _boundary; }
- virtual std::list< Cell* > getBoundary();
- virtual std::map<Cell*, int, Less_Cell > getOrientedCoboundary() { return _coboundary; }
- virtual std::list< Cell* > getCoboundary();
-
- // add (co)boundary cell
- virtual bool addBoundaryCell(int orientation, Cell* cell);
- virtual bool addCoboundaryCell(int orientation, Cell* cell);
-
- // remove (co)boundary cell
- virtual int removeBoundaryCell(Cell* cell, bool other=true);
- virtual int removeCoboundaryCell(Cell* cell, bool other=true);
-
- // true if has given cell on (original) (co)boundary
- virtual bool hasBoundary(Cell* cell, bool org=false);
- virtual bool hasCoboundary(Cell* cell, bool org=false);
-
- virtual void clearBoundary() { _boundary.clear(); }
- virtual void clearCoboundary() { _coboundary.clear(); }
-
- // algebraic dual of the cell
- virtual void makeDualCell();
-
- virtual void printBoundary();
- virtual void printCoboundary();
-
- // original (co)boundary
- virtual void addOrgBdCell(int orientation, Cell* cell) { _obd.insert( std::make_pair(cell, orientation ) ); };
- virtual void addOrgCbdCell(int orientation, Cell* cell) { _ocbd.insert( std::make_pair(cell, orientation ) ); };
- virtual std::map<Cell*, int, Less_Cell > getOrgBd() { return _obd; }
- virtual std::map<Cell*, int, Less_Cell > getOrgCbd() { return _ocbd; }
- virtual void printOrgBd();
- virtual void printOrgCbd();
-
- virtual bool inSubdomain() const { return _inSubdomain; }
- //virtual void setInSubdomain(bool subdomain) { _inSubdomain = subdomain; }
-
- virtual bool onDomainBoundary() const { return _onDomainBoundary; }
- virtual void setOnDomainBoundary(bool domainboundary) { _onDomainBoundary = domainboundary; }
-
- // print cell vertices
- virtual void printCell() const = 0;
-
- virtual int getNumFacets() const { return 0; }
- virtual void getFacetVertices(const int num, std::vector<MVertex*> &v) const {};
-
- // get boundary cell orientation
- virtual int getFacetOri(Cell* cell) { std::vector<MVertex*> v = cell->getVertexVector(); return getFacetOri(v); }
- virtual int getFacetOri(std::vector<MVertex*> &v) { return 0; }
-
- // tools for combined cells
- virtual bool isCombined() { return _combined; }
- virtual std::list< std::pair<int, Cell*> > getCells() { std::list< std::pair<int, Cell*> >cells; cells.push_back( std::make_pair(1, this)); return cells; }
- virtual int getNumCells() {return 1;}
-
- bool operator==(const Cell& c2){
- if(this->getNumVertices() != c2.getNumVertices()){
- return false;
- }
- for(int i=0; i < this->getNumVertices();i++){
- if(this->getSortedVertex(i) != c2.getSortedVertex(i)){
- return false;
- }
- }
- return true;
- }
-
- bool operator<(const Cell& c2) const {
-
- if(this->getNumVertices() != c2.getNumVertices()){
- return (this->getNumVertices() < c2.getNumVertices());
- }
- for(int i=0; i < this->getNumVertices();i++){
- if(this->getSortedVertex(i) < c2.getSortedVertex(i)) return true;
- else if (this->getSortedVertex(i) > c2.getSortedVertex(i)) return false;
- }
- return false;
- }
-
- //virtual void writeMSH(FILE *fp, double version=1.0, bool binary=false, int num=0, int elementary=1, int physical=1, int parentNum=0) {}
-};
-
-
-// Simplicial cell type.
-class Simplex : public Cell
-{
- public:
- Simplex() : Cell() {}
- ~Simplex(){}
-};
-
-// Zero simplex cell type.
-class ZeroSimplex : public Simplex, public MPoint
-{
- private:
-
- public:
-
- 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 getNum() { return MPoint::getNum(); }
- int getType() { return MPoint::getType(); }
- int getTypeForMSH() { return MPoint::getTypeForMSH(); }
- int getPartition() { return MPoint::getPartition(); }
- int getNumVertices() const { return 1; }
- 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<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()); }
-
- void printCell() const { printf("Cell dimension: %d, Vertices: %d, in subdomain: %d \n", getDim(), _v[0]->getNum(), _inSubdomain); }
-
-};
-
-// One simplex cell type.
-class OneSimplex : public Simplex, public MLine
-{
- private:
- // sorted list of vertices
- int _vs[2];
- public:
- 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 getNum() { return MLine::getNum(); }
- int getType() { return MLine::getType(); }
- int getTypeForMSH() { return MLine::getTypeForMSH(); }
- int getPartition() { return MLine::getPartition(); }
- int getNumVertices() const { return 2; }
- int getNumFacets() const { return 2; }
- MVertex* getVertex(int vertex) const {return _v[vertex]; }
- int getSortedVertex(int vertex) const {return _vs[vertex]; }
- bool hasVertex(int vertex) const {return (_v[0]->getNum() == vertex || _v[1]->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<MVertex*> &v) const {
- v.resize(1);
- v[0] = _v[num];
- }
- int getFacetOri(std::vector<MVertex*> &v){
- if(v.size() != 1) return 0;
- else if(v[0] == _v[0]) return -1;
- else if(v[0] == _v[1]) return 1;
- else return 0;
- }
-
- void printCell() const { printf("Cell dimension: %d, Vertices: %d %d, in subdomain: %d \n", getDim(), _v[0]->getNum(), _v[1]->getNum(), _inSubdomain); }
-
-};
-
-// Two simplex cell type.
-class TwoSimplex : public Simplex, public MTriangle
-{
- private:
- // sorted list of vertices
- int _vs[3];
-
- public:
-
- 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 getNum() { return MTriangle::getNum(); }
- int getType() { return MTriangle::getType(); }
- int getTypeForMSH() { return MTriangle::getTypeForMSH(); }
- int getPartition() { return MTriangle::getPartition(); }
- int getNumVertices() const { return 3; }
- int getNumFacets() const { return 3; }
- MVertex* getVertex(int vertex) const {return _v[vertex]; }
- int getSortedVertex(int vertex) const {return _vs[vertex]; }
- bool hasVertex(int vertex) const {return
- (_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<MVertex*> &v) const {
- MTriangle::getEdgeVertices(num, v);
- }
- int getFacetOri(std::vector<MVertex*> &v){
- if(v.size() != 2) return 0;
- MEdge facet = MEdge(v[0], v[1]);
- int ithFacet = 0;
- int sign = 0;
- MTriangle::getEdgeInfo(facet, ithFacet, sign);
- return sign;
- }
-
- void printCell() const { printf("Cell dimension: %d, Vertices: %d %d %d, in subdomain: %d\n", getDim(), _v[0]->getNum(), _v[1]->getNum(), _v[2]->getNum(), _inSubdomain); }
-};
-
-// Quadrangle cell type.
-class CQuadrangle : public Cell, public MQuadrangle
-{
- private:
- // sorted list of vertices
- int _vs[4];
-
- public:
-
- CQuadrangle(std::vector<MVertex*> v, int tag=0, bool subdomain=false, bool boundary=false, int num=0, int part=0)
- : MQuadrangle(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+3);
- }
-
- ~CQuadrangle(){}
-
- int getDim() const { return 2; }
- int getNum() { return MQuadrangle::getNum(); }
- int getType() { return MQuadrangle::getType(); }
- int getTypeForMSH() { return MQuadrangle::getTypeForMSH(); }
- int getPartition() { return MQuadrangle::getPartition(); }
- int getNumVertices() const { return 4; }
- int getNumFacets() const { return 4; }
- MVertex* getVertex(int vertex) const {return _v[vertex]; }
- int getSortedVertex(int vertex) const {return _vs[vertex]; }
- bool hasVertex(int vertex) const {return
- (_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<MVertex*> &v) const {
- MQuadrangle::getEdgeVertices(num, v);
- }
- int getFacetOri(std::vector<MVertex*> &v){
- if(v.size() != 2) return 0;
- MEdge facet = MEdge(v[0], v[1]);
- int ithFacet = 0;
- int sign = 0;
- MQuadrangle::getEdgeInfo(facet, ithFacet, sign);
- return sign;
- }
-
- void printCell() const { printf("Cell dimension: %d, Vertices: %d %d %d %d, in subdomain: %d\n",
- getDim(), _v[0]->getNum(), _v[1]->getNum(), _v[2]->getNum(),
- _v[3]->getNum(), _inSubdomain); }
-};
-
-
-// ThreeSimplex cell type.
-class ThreeSimplex : public Simplex, public MTetrahedron
-{
- private:
- // sorted list of vertices
- int _vs[4];
-
- public:
- 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 getNum() { return MTetrahedron::getNum(); }
- int getType() { return MTetrahedron::getType(); }
- int getTypeForMSH() { return MTetrahedron::getTypeForMSH(); }
- int getPartition() { return MTetrahedron::getPartition(); }
- int getNumVertices() const { return 4; }
- int getNumFacets() const { return 4; }
- MVertex* getVertex(int vertex) const {return _v[vertex]; }
- int getSortedVertex(int vertex) const {return _vs[vertex]; }
- bool hasVertex(int vertex) const {return
- (_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<MVertex*> &v) const {
- MTetrahedron::getFaceVertices(num, v);
- }
- int getFacetOri(std::vector<MVertex*> &v){
- if(v.size() != 3) return 0;
- MFace facet = MFace(v);
- int ithFacet = 0;
- int sign = 0;
- int rot = 0;
- MTetrahedron::getFaceInfo(facet, ithFacet, sign, rot);
- return sign;
- }
-
- virtual void printCell() const { printf("Cell dimension: %d, Vertices: %d %d %d %d, in subdomain: %d \n", getDim(), _v[0]->getNum(), _v[1]->getNum(), _v[2]->getNum(), _v[3]->getNum(), _inSubdomain); }
-};
-
-// Hexahedron cell type.
-class CHexahedron : public Cell, public MHexahedron
-{
- private:
- // sorted list of vertices
- int _vs[8];
-
- public:
- CHexahedron(std::vector<MVertex*> v, int tag=0, bool subdomain=false, bool boundary=false, int num=0, int part=0)
- : MHexahedron(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();
- _vs[4] = v.at(4)->getNum();
- _vs[5] = v.at(5)->getNum();
- _vs[6] = v.at(6)->getNum();
- _vs[7] = v.at(7)->getNum();
- std::sort(_vs, _vs+8);
- }
-
- ~CHexahedron(){}
-
- int getDim() const { return 3; }
- int getNum() { return MHexahedron::getNum(); }
- int getType() { return MHexahedron::getType(); }
- int getTypeForMSH() { return MHexahedron::getTypeForMSH(); }
- int getPartition() { return MHexahedron::getPartition(); }
- int getNumVertices() const { return 8; }
- int getNumFacets() const { return 6; }
- MVertex* getVertex(int vertex) const {return _v[vertex]; }
- int getSortedVertex(int vertex) const {return _vs[vertex]; }
- bool hasVertex(int vertex) const {return
- (_v[0]->getNum() == vertex || _v[1]->getNum() == vertex || _v[2]->getNum() == vertex || _v[3]->getNum() == vertex ||
- _v[4]->getNum() == vertex || _v[5]->getNum() == vertex || _v[6]->getNum() == vertex || _v[7]->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<MVertex*> &v) const {
- MHexahedron::getFaceVertices(num, v);
- }
- int getFacetOri(std::vector<MVertex*> &v){
- if(v.size() != 4) return 0;
- MFace facet = MFace(v);
- int ithFacet = 0;
- int sign = 0;
- int rot = 0;
- MHexahedron::getFaceInfo(facet, ithFacet, sign, rot);
- return sign;
- }
-
- virtual void printCell() const { printf("Cell dimension: %d, Vertices: %d %d %d %d %d %d %d %d, in subdomain: %d \n", getDim(), _v[0]->getNum(), _v[1]->getNum(), _v[2]->getNum(), _v[3]->getNum(), _v[4]->getNum(), _v[5]->getNum(), _v[6]->getNum(), _v[7]->getNum(), _inSubdomain); }
-};
-
-
-// Prism cell type.
-class CPrism : public Cell, public MPrism
-{
- private:
- // sorted list of vertices
- int _vs[6];
-
- public:
- CPrism(std::vector<MVertex*> v, int tag=0, bool subdomain=false, bool boundary=false, int num=0, int part=0)
- : MPrism(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();
- _vs[4] = v.at(4)->getNum();
- _vs[5] = v.at(5)->getNum();
- std::sort(_vs, _vs+6);
- }
-
- ~CPrism(){}
-
- int getDim() const { return 3; }
- int getNum() { return MPrism::getNum(); }
- int getType() { return MPrism::getType(); }
- int getTypeForMSH() { return MPrism::getTypeForMSH(); }
- int getPartition() { return MPrism::getPartition(); }
- int getNumVertices() const { return 6; }
- int getNumFacets() const { return 5; }
- MVertex* getVertex(int vertex) const {return _v[vertex]; }
- int getSortedVertex(int vertex) const {return _vs[vertex]; }
- bool hasVertex(int vertex) const {return
- (_v[0]->getNum() == vertex || _v[1]->getNum() == vertex || _v[2]->getNum() == vertex || _v[3]->getNum() == vertex ||
- _v[4]->getNum() == vertex || _v[5]->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<MVertex*> &v) const {
- MPrism::getFaceVertices(num, v);
- }
- int getFacetOri(std::vector<MVertex*> &v){
- if(v.size() != 4 && v.size() != 3) return 0;
- MFace facet = MFace(v);
- int ithFacet = 0;
- int sign = 0;
- int rot = 0;
- MPrism::getFaceInfo(facet, ithFacet, sign, rot);
- return sign;
- }
-
- virtual void printCell() const { printf("Cell dimension: %d, Vertices: %d %d %d %d %d %d, in subdomain: %d \n", getDim(), _v[0]->getNum(), _v[1]->getNum(), _v[2]->getNum(), _v[3]->getNum(), _v[4]->getNum(), _v[5]->getNum(), _inSubdomain); }
-};
-
-
-// Pyramid cell type.
-class CPyramid : public Cell, public MPyramid
-{
- private:
- // sorted list of vertices
- int _vs[5];
-
- public:
- CPyramid(std::vector<MVertex*> v, int tag=0, bool subdomain=false, bool boundary=false, int num=0, int part=0)
- : MPyramid(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();
- _vs[4] = v.at(4)->getNum();
- std::sort(_vs, _vs+5);
- }
-
- ~CPyramid(){}
-
- int getDim() const { return 3; }
- int getNum() { return MPyramid::getNum(); }
- int getType() { return MPyramid::getType(); }
- int getTypeForMSH() { return MPyramid::getTypeForMSH(); }
- int getPartition() { return MPyramid::getPartition(); }
- int getNumVertices() const { return 5; }
- int getNumFacets() const { return 5; }
- MVertex* getVertex(int vertex) const {return _v[vertex]; }
- int getSortedVertex(int vertex) const {return _vs[vertex]; }
- bool hasVertex(int vertex) const {return
- (_v[0]->getNum() == vertex || _v[1]->getNum() == vertex || _v[2]->getNum() == vertex || _v[3]->getNum() == vertex ||
- _v[4]->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<MVertex*> &v) const {
- MPyramid::getFaceVertices(num, v);
- }
- int getFacetOri(std::vector<MVertex*> &v){
- if(v.size() != 4 && v.size() != 3) return 0;
- MFace facet = MFace(v);
- int ithFacet = 0;
- int sign = 0;
- int rot = 0;
- MPyramid::getFaceInfo(facet, ithFacet, sign, rot);
- return sign;
- }
-
- virtual void printCell() const { printf("Cell dimension: %d, Vertices: %d %d %d %d %d, in subdomain: %d \n", getDim(), _v[0]->getNum(), _v[1]->getNum(), _v[2]->getNum(), _v[3]->getNum(), _v[4]->getNum(), _inSubdomain); }
-};
-
-// Ordering for the finite element mesh vertices.
-class Less_MVertex{
- public:
- bool operator()(const MVertex* v1, const MVertex* v2) const {
- return (v1->getNum() < v2->getNum());
- }
-};
-
-// A cell that is a combination of cells of same dimension
-class CombinedCell : public Cell{
-
- private:
- // vertices
- std::vector<MVertex*> _v;
- // sorted list of all vertices
- std::vector<int> _vs;
- // list of cells this cell is a combination of
- std::list< std::pair<int, Cell*> > _cells;
- int _num;
-
- public:
-
- CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co=false);
- ~CombinedCell();
-
- int getNum() { return _num; }
- int getNumVertices() const { return _v.size(); }
- MVertex* getVertex(int vertex) const { return _v.at(vertex); }
- int getSortedVertex(int vertex) const { return _vs.at(vertex); }
- std::vector<MVertex*> getVertexVector() const { return _v; }
-
- // true if this cell has given vertex
- bool hasVertex(int vertex) const;
-
- virtual void printCell() const;
-
- std::list< std::pair<int, Cell*> > getCells() { return _cells; }
- int getNumCells() {return _cells.size();}
-
-};
-
// A class representing a cell complex made out of a finite element mesh.
class CellComplex
{
@@ -690,9 +41,6 @@ class CellComplex
// entities on the boundary of the homology computation domain
std::vector<GEntity*> _boundary;
- // mesh vertices in this domain
- std::set<MVertex*, Less_MVertex> _domainVertices;
-
// sorted containers of unique cells in this cell complex
// one for each dimension
std::set<Cell*, Less_Cell> _cells[4];
@@ -720,7 +68,6 @@ class CellComplex
// for constructor
void insert_cells(bool subdomain, bool boundary);
void find_boundary(std::vector<GEntity*>& domain, std::vector<GEntity*>& subdomain);
- void find_vertices(std::vector<GEntity*>& domain, std::vector<GEntity*>& subdomain);
// remove a cell from this cell complex
void removeCell(Cell* cell, bool other=true);
@@ -784,7 +131,7 @@ class CellComplex
// write this cell complex in 2.0 MSH ASCII format
// for debugging only
- int writeComplexMSH(const std::string &name);
+ //int writeComplexMSH(const std::string &name);
// Cell combining for reduction and coreduction
int combine(int dim);
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index 73c36f04bacfe481ab406900ef330798a0e8624d..231bee797bbc7685b6601040e13f9ec92fc195cf 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -566,8 +566,6 @@ void Chain::smoothenChain(){
int Chain::writeChainMSH(const std::string &name){
- //_cellComplex->writeComplexMSH(name);
-
if(getSize() == 0) return 1;
FILE *fp = fopen(name.c_str(), "a");
@@ -610,20 +608,11 @@ void Chain::createPView(){
for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
Cell* cell = (*cit).first;
int coeff = (*cit).second;
- std::vector<MVertex*> v = cell->getVertexVector();
- if(cell->getDim() > 0 && coeff < 0){ // flip orientation
- if(getDim() != 1){
- MVertex* temp = v[1];
- v[1] = v[v.size()-1];
- v[v.size()-1] = temp;
- }
- else{
- MVertex* temp = v[0];
- v[0] = v[1];
- v[1] = temp;
- }
- }
- MElement *e = factory.create(cell->getTypeForMSH(), v, cell->getNum(), cell->getPartition());
+
+ MElement *e = cell->getImageMElement();
+ cell->setDeleteImage(false);
+
+ if(cell->getDim() > 0 && coeff < 0) e->revert(); // flip orientation
for(int i = 0; i < abs(coeff); i++) elements.push_back(e);
std::vector<double> coeffs (1,abs(coeff));
@@ -638,11 +627,9 @@ void Chain::createPView(){
setNum(physicalNum);
std::map<int, std::vector<MElement*> > entityMap;
- //int entityNum = _model->getMaxElementaryNumber(getDim())+1;
entityMap[entityNum] = elements;
std::map<int, std::map<int, std::string> > physicalMap;
std::map<int, std::string> physicalInfo;
- //int physicalNum = _model->getMaxPhysicalNumber(getDim())+1;
physicalInfo[physicalNum]=getName();
physicalMap[entityNum] = physicalInfo;
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index 0a8d7bc04cd617d86c78c779add13b08f33219c7..fe8850d0dbfc146d6727a6737a6fe8a59d8f2411 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -94,9 +94,7 @@ void Homology::findGenerators(std::string fileName){
int omitted = _cellComplex->reduceComplex();
//_cellComplex->printEuler();
- //_cellComplex->emptyTrash();
-
_cellComplex->combine(3);
_cellComplex->reduction(2);
_cellComplex->combine(2);
@@ -106,7 +104,6 @@ void Homology::findGenerators(std::string fileName){
_cellComplex->checkCoherence();
//_cellComplex->printEuler();
- //_cellComplex->emptyTrash();
double t2 = Cpu();
Msg::Info("Cell Complex reduction complete (%g s).", t2 - t1);
diff --git a/tutorial/t10.geo b/tutorial/t10.geo
index a97eb68816713085c9b1090adb2142686598260b..6e8606054f0a3495490e4da8ee6056bc0c835745 100644
--- a/tutorial/t10.geo
+++ b/tutorial/t10.geo
@@ -8,7 +8,7 @@
// Homology computation in Gmsh finds representative chains of (relative) homology spaces using a mesh of a model. Those representatives generate the (relative) homology spaces of the model. Alternatively, Gmsh can only look for the ranks of the (relative) homology spaces, the Betti numbers of the model.
-// The generators chains are stored in a given .msh-file as physical groups, whose mesh elements are oriented such that their coeffecients are 1 in the generator chain.
+// The generators chains are stored in a given .msh-file as physical groups, whose mesh elements are oriented such that their coefficients are 1 in the generator chain.
// Create an example geometry