diff --git a/Geo/Cell.cpp b/Geo/Cell.cpp
index 906262739f5b77b7ea8e471c681d58a7f875e963..6f97ea643498802b85dd4c35cd5e839d320347c3 100755
--- a/Geo/Cell.cpp
+++ b/Geo/Cell.cpp
@@ -103,7 +103,7 @@ int Cell::getFacetOri(std::vector<MVertex*> &v)
void Cell::printCell()
{
- printf("%d-cell %d: \n" , getDim(), getNum());
+ printf("%d-cell: \n" , getDim());
printf("Vertices: ");
for(int i = 0; i < this->getNumVertices(); i++){
printf("%d ", this->getSortedVertex(i));
@@ -205,14 +205,6 @@ bool Cell::hasCoboundary(Cell* cell, bool org)
}
}
-void Cell::makeDualCell()
-{
- std::map<Cell*, int, Less_Cell > temp = _boundary;
- _boundary = _coboundary;
- _coboundary = temp;
- _dim = 3-_dim;
-}
-
void Cell::printBoundary(bool org)
{
for(biter it = firstBoundary(org); it != lastBoundary(org); it++){
diff --git a/Geo/Cell.h b/Geo/Cell.h
index 60d0b1742c201d6e3f3edf9bb0b473e35d577d4f..ccf6213b31d3937356ff2a2993b83016c29bf304 100644
--- a/Geo/Cell.h
+++ b/Geo/Cell.h
@@ -87,10 +87,6 @@ class Cell
virtual int getDim() const { return _dim; };
virtual int getIndex() const { return _index; };
virtual void setIndex(int index) { _index = index; };
- virtual int getNum() const { return _image->getNum(); }
- virtual int getType() const { return _image->getType(); }
- virtual int getTypeForMSH() const { return _image->getTypeForMSH(); }
- virtual int getPartition() const { return _image->getPartition(); }
virtual void setImmune(bool immune) { _immune = immune; };
virtual bool getImmune() const { return _immune; };
virtual bool inSubdomain() const { return _inSubdomain; }
@@ -146,10 +142,7 @@ class Cell
virtual void clearBoundary() { _boundary.clear(); }
virtual void clearCoboundary() { _coboundary.clear(); }
- // algebraic dual of the cell
- virtual void makeDualCell();
-
- // print cell info
+ // print cell debug info
virtual void printCell();
virtual void printBoundary(bool org=false);
virtual void printCoboundary(bool org=false);
@@ -160,14 +153,14 @@ class Cell
virtual void getFacetVertices(const int num, std::vector<MVertex*> &v) const;
// get boundary cell orientation
+ virtual int getFacetOri(std::vector<MVertex*> &v);
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; }
diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index c38226ccad226d9f54462d8711469fae03c676da..49e5b626001d2cfa10883af0279c0881245608ba 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -9,37 +9,21 @@
#if defined(HAVE_KBIPACK)
-CellComplex::CellComplex( std::vector<GEntity*> domain,
- std::vector<GEntity*> subdomain)
+CellComplex::CellComplex( std::vector<MElement*>& domainElements,
+ std::vector<MElement*>& subdomainElements)
{
- _domain = domain;
- _subdomain = subdomain;
-
_dim = 0;
_simplicial = true;
-
- _multidim = false;
- int dim = 0;
- for(unsigned int i = 0; i < domain.size(); i++){
- GEntity* entity = domain.at(i);
- if(i == 0) dim = entity->dim();
- if(dim != entity->dim()){
- _multidim = true;
- //printf("Warning: Domain is not a manifold.");
- break;
- }
- }
-
+
// insert cells into cell complex
- // subdomain need to be inserted first!
- if(!insert_cells(_subdomain, true)){ panic_exit(); return; }
- if(!insert_cells(_domain, false)) { panic_exit(); return; }
+ // subdomain needs to be inserted first!
+ if(!insert_cells(subdomainElements, true)){ panic_exit(); return; }
+ if(!insert_cells(domainElements, false)) { panic_exit(); return; }
for(int i = 0; i < 4; i++){
_ocells[i] = _cells[i];
if(getSize(i) > _dim) _dim = i;
- }
-
+ }
}
void CellComplex::panic_exit(){
@@ -53,68 +37,50 @@ void CellComplex::panic_exit(){
}
for(unsigned int i = 0; i < _newcells.size(); i++) delete _newcells.at(i);
_newcells.clear();
- _domain.clear();
- _subdomain.clear();
}
-bool CellComplex::insert_cells(std::vector<GEntity*>& domain, bool subdomain)
+bool CellComplex::insert_cells(std::vector<MElement*>& elements,
+ bool subdomain)
{
std::vector<MVertex*> vertices;
std::pair<citer, bool> insertInfo;
// add highest dimensional cells
- for(unsigned int j=0; j < domain.size(); j++) {
- for(unsigned int i=0; i < domain.at(j)->getNumMeshElements(); i++){
- vertices.clear();
- MElement* element = domain.at(j)->getMeshElement(i);
-
- for(int k=0; k < element->getNumVertices(); k++){
- MVertex* vertex = element->getVertex(k);
- vertices.push_back(vertex);
- }
-
- int dim = element->getDim();
- int type = element->getTypeForMSH();
+ for(unsigned int i=0; i < elements.size(); i++){
+ vertices.clear();
+ MElement* element = elements.at(i);
+
+ for(int k=0; k < element->getNumVertices(); k++){
+ MVertex* vertex = element->getVertex(k);
+ vertices.push_back(vertex);
+ }
- Cell* cell;
- // simplex types
- if(type == MSH_LIN_2 || type == MSH_TRI_3 || type == MSH_TET_4
- || type == MSH_LIN_3 || type == MSH_TRI_6 || type == MSH_TET_10
- || type == MSH_PNT || type == MSH_TRI_9 || type == MSH_TRI_10
- || type == MSH_TRI_12 || 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 ){
- cell = new Cell(element);
- }
- else if(type == MSH_QUA_4 || type == MSH_QUA_8 || type == MSH_QUA_9){
- cell = new Cell(element);
- _simplicial = false;
- }
- else if(type == MSH_HEX_8 || type == MSH_HEX_27 || type == MSH_HEX_20){
- cell = new Cell(element);
- _simplicial = false;
- }
- else if(type == MSH_PRI_6 || type == MSH_PRI_18 || type == MSH_PRI_15){
- cell = new Cell(element);
- _simplicial = false;
- }/* FIXME: no getFaceInfo methods for these MElements
- else if(type == MSH_PYR_5 || type == MSH_PYR_14 || type == MSH_PYR_13){
- cell = new Cell(element, subdomain, boundary);
- _simplicial = false;
- }*/
- else {
- //printf("Error: mesh element %d not implemented yet! \n", type);
- return false;
- }
- cell->setImmune(false);
- cell->setInSubdomain(subdomain);
- insertInfo = _cells[dim].insert(cell);
- if(!insertInfo.second) delete cell;
+ int dim = element->getDim();
+ int type = element->getTypeForMSH();
+
+ Cell* cell = new Cell(element);
+ // simplex types
+ if( !(type == MSH_PNT
+ || type == MSH_LIN_2 || type == MSH_TRI_3 || type == MSH_TET_4
+ || type == MSH_LIN_3 || type == MSH_TRI_6 || type == MSH_TET_10
+ || type == MSH_PNT || type == MSH_TRI_9 || type == MSH_TRI_10
+ || type == MSH_TRI_12 || 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) ){
+ _simplicial = false;
}
+ /* FIXME: no getFaceInfo methods for these MElements */
+ if(type == MSH_PYR_5 || type == MSH_PYR_14 || type == MSH_PYR_13){
+ //printf("Error: mesh element %d not implemented yet! \n", type);
+ return false;
+ }
+ cell->setImmune(false);
+ cell->setInSubdomain(subdomain);
+ insertInfo = _cells[dim].insert(cell);
+ if(!insertInfo.second) delete cell;
}
-
+
// add lower dimensional cells recursively
MElementFactory factory;
for (int dim = 3; dim > 0; dim--){
@@ -123,7 +89,7 @@ bool CellComplex::insert_cells(std::vector<GEntity*>& domain, bool subdomain)
std::vector<MVertex*> vertices;
for(int i = 0; i < cell->getNumFacets(); i++){
cell->getFacetVertices(i, vertices);
- int type = cell->getTypeForMSH();
+ int type = cell->getImageMElement()->getTypeForMSH();
int newtype = 0;
//FIXME: add missing boundary cell type relations
//FIXME: high order meshes don't work
@@ -131,30 +97,34 @@ bool CellComplex::insert_cells(std::vector<GEntity*>& domain, bool subdomain)
if(type == MSH_TET_4) newtype = MSH_TRI_3;
/*else if(type == MSH_TET_10) newtype = MSH_TRI_6;
else if(type == MSH_TET_20) newtype = MSH_TRI_9;*/
- else if(type == MSH_HEX_8) newtype = MSH_QUA_4;
+ /*else if(type == MSH_HEX_8) newtype = MSH_QUA_4;*/
/*else if(type == MSH_HEX_20) newtype = MSH_QUA_8;
else if(type == MSH_HEX_27) newtype = MSH_QUA_9;*/
- else if(type == MSH_PRI_6 && vertices.size() == 3) newtype = MSH_TRI_3;
- else if(type == MSH_PRI_6 && vertices.size() == 4) newtype = MSH_QUA_4;
+ else if(type == MSH_PRI_6
+ && vertices.size() == 3) newtype = MSH_TRI_3;
+ else if(type == MSH_PRI_6
+ && vertices.size() == 4) newtype = MSH_QUA_4;
}
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(type == MSH_TRI_9) newtype = MSH_LIN_4;
- else if(type == MSH_QUA_9) newtype = MSH_LIN_3;*/
- }
- else if(dim == 1){
- if(type == MSH_LIN_2) newtype = MSH_PNT;
- /*else if(type == MSH_LIN_3 || type == MSH_LIN_4 ||
- type == MSH_LIN_5 || type == MSH_LIN_6) newtype = MSH_PNT;*/
- }
- if(newtype == 0){
+ 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(type == MSH_TRI_9) newtype = MSH_LIN_4;
+ else if(type == MSH_QUA_9) newtype = MSH_LIN_3;*/
+ }
+ else if(dim == 1){
+ if(type == MSH_LIN_2) newtype = MSH_PNT;
+ /*else if(type == MSH_LIN_3 || type == MSH_LIN_4 ||
+ type == MSH_LIN_5 || type == MSH_LIN_6) newtype = MSH_PNT;*/
+ }
+ if(newtype == 0){
//printf("Error: mesh element %d not implemented yet! \n", type);
return false;
}
MElement* element = factory.create(newtype, vertices, 0,
- cell->getPartition());
+ cell->getImageMElement()->
+ getPartition());
Cell* newCell = new Cell(element);
newCell->setImmune(cell->getImmune());
newCell->setInSubdomain(subdomain);
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index 05586ac5ea7d909acf28dc865b08fa8a632ba31f..30df99cfeae572e2fa2cbd8cda364667374560c9 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -20,7 +20,7 @@
#include "Cell.h"
#include "MElement.h"
#include "ChainComplex.h"
-#include "GModel.h"
+//#include "GModel.h"
class Cell;
@@ -28,16 +28,6 @@ class Cell;
class CellComplex
{
private:
-
- // the domain in the model which this cell complex covers
- std::vector<GEntity*> _domain;
- // a subdomain of the given domain
- // used in relative homology computation, may be empty
- std::vector<GEntity*> _subdomain;
-
- // entities on the boundary of the homology computation domain
- std::vector<GEntity*> _boundary;
-
// sorted containers of unique cells in this cell complex
// one for each dimension
std::set<Cell*, Less_Cell> _cells[4];
@@ -56,9 +46,6 @@ class CellComplex
// is the cell complex simplicial
bool _simplicial;
- // does the domain contain entities of different dimensions
- bool _multidim;
-
// enqueue cells in queue if they are not there already
void enqueueCells(std::map<Cell*, int, Less_Cell>& cells,
std::queue<Cell*>& Q, std::set<Cell*, Less_Cell>& Qset);
@@ -66,7 +53,7 @@ class CellComplex
void removeCellQset(Cell* cell, std::set<Cell*, Less_Cell>& Qset);
// for constructor
- bool insert_cells(std::vector<GEntity*>& domain, bool subdomain);
+ bool insert_cells(std::vector<MElement*>& elements, bool subdomain);
void panic_exit();
// insert/remove a cell from this cell complex
@@ -78,13 +65,10 @@ class CellComplex
public:
- CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> subdomain);
- //CellComplex(CellComplex* cellComplex)
+ CellComplex( std::vector<MElement*>& domainElements,
+ std::vector<MElement*>& subdomainElements);
~CellComplex();
- std::vector<GEntity*> getDomain() const { return _domain; }
- std::vector<GEntity*> getSubdomain() const { return _subdomain; }
-
// get the number of certain dimensional cells
int getSize(int dim, bool org=false){
if(!org) return _cells[dim].size();
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index 929e3e838268b5087b7d4b2a2058918b97ed0f88..7b0cd2bab4cb4d5c95be6d34f7ab99c39c27e83a 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -69,8 +69,24 @@ CellComplex* Homology::createCellComplex(std::vector<GEntity*>& domainEntities,
if(domainEntities.empty()) Msg::Error("Domain is empty.");
if(subdomainEntities.empty()) Msg::Info("Subdomain is empty.");
- CellComplex* cellComplex = new CellComplex(domainEntities,
- subdomainEntities);
+ std::vector<MElement*> domainElements;
+ std::vector<MElement*> subdomainElements;
+ for(unsigned int j=0; j < domainEntities.size(); j++) {
+ for(unsigned int i=0; i < domainEntities.at(j)->getNumMeshElements(); i++){
+ MElement* element = domainEntities.at(j)->getMeshElement(i);
+ domainElements.push_back(element);
+ }
+ }
+ for(unsigned int j=0; j < subdomainEntities.size(); j++) {
+ for(unsigned int i=0; i < subdomainEntities.at(j)->getNumMeshElements();
+ i++){
+ MElement* element = subdomainEntities.at(j)->getMeshElement(i);
+ subdomainElements.push_back(element);
+ }
+ }
+
+ CellComplex* cellComplex = new CellComplex(domainElements,
+ subdomainElements);
if(cellComplex->getSize(0) == 0){
Msg::Error("Cell Complex is empty!");
@@ -94,7 +110,8 @@ Homology::~Homology()
void Homology::findGenerators()
{
- CellComplex* cellComplex = createCellComplex(_domainEntities, _subdomainEntities);
+ CellComplex* cellComplex = createCellComplex(_domainEntities,
+ _subdomainEntities);
std::string domainString = getDomainString(_domain, _subdomain);
Msg::Info("Reducing the Cell Complex...");
@@ -320,9 +337,8 @@ void Homology::findDualGenerators()
}
void Homology::findHomSequence(){
-
- CellComplex* cellComplex = new CellComplex(_domainEntities,
- _subdomainEntities);
+ CellComplex* cellComplex = createCellComplex(_domainEntities,
+ _subdomainEntities);
Msg::Info("Reducing the Cell Complex...");
Msg::StatusBar(1, false, "Reducing...");
double t1 = Cpu();
@@ -718,7 +734,7 @@ int Chain::writeChainMSH(const std::string &name)
for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
Cell* cell = (*cit).first;
int coeff = (*cit).second;
- fprintf(fp, "%d %d \n", cell->getNum(), coeff );
+ fprintf(fp, "%d %d \n", cell->getImageMElement()->getNum(), coeff );
}
fprintf(fp, "$EndElementData\n");
diff --git a/Geo/Homology.h b/Geo/Homology.h
index 88a630ea2ccf56b1a961f765b5b92d2076a5cf5e..705a27bfff223d442d52a2aba219363fef3ea40e 100644
--- a/Geo/Homology.h
+++ b/Geo/Homology.h
@@ -12,6 +12,7 @@
#include "CellComplex.h"
#include "ChainComplex.h"
#include "OS.h"
+#include "GModel.h"
#include "PView.h"
#include "Options.h"