From 3453e93e192b682fc2bfe598ffeac25a81c3fbd3 Mon Sep 17 00:00:00 2001
From: Matti Pellika <matti.pellikka@tut.fi>
Date: Sun, 17 Jan 2010 15:24:15 +0000
Subject: [PATCH] A major renovation of the Cell class. Rather than imitating
MElements, each Cell maps to an existing or newly created MElement and uses
MElement's services to their extend. This brings simplicity and memory
efficiency.
---
Geo/Cell.cpp | 146 +++++-----
Geo/Cell.h | 247 ++++++++++++++++
Geo/CellComplex.cpp | 282 ++++---------------
Geo/CellComplex.h | 657 +------------------------------------------
Geo/ChainComplex.cpp | 23 +-
Geo/Homology.cpp | 3 -
tutorial/t10.geo | 2 +-
7 files changed, 382 insertions(+), 978 deletions(-)
create mode 100644 Geo/Cell.h
diff --git a/Geo/Cell.cpp b/Geo/Cell.cpp
index 5cdaaba5cf..962a53b3ad 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 0000000000..3fb2439e83
--- /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 8b49f25cee..cb2b779a5a 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 ff2bd3b08c..49c9ffc072 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 73c36f04ba..231bee797b 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 0a8d7bc04c..fe8850d0db 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 a97eb68816..6e8606054f 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
--
GitLab