diff --git a/Geo/Cell.h b/Geo/Cell.h
index 94e771b480d8d459af4b8c5dcc3457698c5c1642..fe02f850c076b178404f4e8779e2915f7595df5b 100644
--- a/Geo/Cell.h
+++ b/Geo/Cell.h
@@ -16,10 +16,12 @@
#include <string>
#include <algorithm>
#include <set>
+#include <list>
+#include <map>
#include <queue>
#include "CellComplex.h"
#include "MElement.h"
-#include "GModel.h"
+//#include "GModel.h"
class Cell;
diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index d6e0918e503a0b9318b62d7be05fd04e60b19ee4..7d99b7f51b9c5173f60c1bb4e8c0fbd890910fbb 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -6,13 +6,9 @@
// Contributed by Matti Pellikka <matti.pellikka@tut.fi>.
#include "CellComplex.h"
-#include "OS.h"
-#include "Context.h"
#if defined(HAVE_KBIPACK)
-
-
CellComplex::CellComplex( std::vector<GEntity*> domain,
std::vector<GEntity*> subdomain)
{
@@ -140,13 +136,6 @@ void CellComplex::insert_cells(std::vector<GEntity*>& domain,
int dim = element->getDim();
int type = element->getTypeForMSH();
- if(CTX::instance()->mesh.radiusSup){
- double max = CTX::instance()->mesh.radiusSup;
- double min = CTX::instance()->mesh.radiusInf;
- double r = element->maxEdge();
- if(r < min || r > max) continue;
- }
-
Cell* cell;
// simplex types
if(type == MSH_LIN_2 || type == MSH_TRI_3 || type == MSH_TET_4
@@ -204,6 +193,8 @@ void CellComplex::insert_cells(std::vector<GEntity*>& domain,
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(dim == 2){
if(type == MSH_TRI_3 || type == MSH_QUA_4) newtype = MSH_LIN_2;
@@ -260,21 +251,6 @@ CellComplex::~CellComplex()
_newcells.clear();
}
-/*
-CellComplex::CellComplex(CellComplex* cellComplex)
-{
- _domain = cellComplex->getDomain();
- _subdomain = cellComplex->getSubdomain();
- _boundary = cellComplex->getBoundary();
-
- for(int i = 0; i < 4; i++){
- _betti[i] = cellComplex->getBetti(i);
- _ocells[i] = cellComplex->getOcells(i);
- }
- _dim = cellComplex->getDim();
- restoreComplex();
- }*/
-
void CellComplex::insertCell(Cell* cell)
{
_newcells.push_back(cell);
@@ -377,9 +353,7 @@ int CellComplex::reduction(int dim, int omitted)
while(cit != lastCell(dim-1)){
Cell* cell = *cit;
if( cell->getCoboundarySize() == 1
- && inSameDomain(cell, cell->firstCoboundary()->first)
- && !cell->getImmune()
- && !cell->firstCoboundary()->first->getImmune()){
+ && inSameDomain(cell, cell->firstCoboundary()->first)){
++cit;
if(dim == getDim() && omitted > 0){
_store.at(omitted-1).insert(cell->firstCoboundary()->first);
@@ -477,30 +451,6 @@ void CellComplex::removeSubdomain()
}
for(unsigned int i = 0; i < toRemove.size(); i++) removeCell(toRemove[i]);
}
-/*
-void CellComplex::setSubdomain(std::vector<GEntity*>& subdomain){
- for(int i = 3; i > -1; i--){
- for(citer cit = firstCell(i); cit != lastCell(i); cit++){
- Cell* cell = *cit;
- for(unsigned int j = 0; j < subdomain.size(); j++){
- if(cell->getGEntity()->tag() == subdomain.at(j)->tag()){
- cell->setInSubdomain(true);
- }
- }
- }
- }
-}*/
-
-void CellComplex::deImmuneCells(bool subdomain)
-{
- for(int i = 0; i < 4; i++){
- for(citer cit = firstCell(i); cit != lastCell(i); ++cit){
- Cell *cell = *cit;
- if(cell->getImmune() && subdomain) cell->setInSubdomain(true);
- cell->setImmune(false);
- }
- }
-}
int CellComplex::coreduceComplex()
{
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index daef3684dcc7cc5077e57215de03666c5752b93d..8a77c8cc38622217ee1b35139cbd1ce0e0bf8a18 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -20,6 +20,8 @@
#include "Cell.h"
#include "MElement.h"
#include "ChainComplex.h"
+#include "OS.h"
+#include "GModel.h"
class Cell;
@@ -177,6 +179,7 @@ class CellComplex
// store cells of dimension dim to the GModel as a
// physical group of MElements
+ // for debugging purposes
void storeCells(int dim);
};
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index 6960108497b6d906e0fa980c4cb1821a622269fd..1dd04bb090ac6cf978e76af155d740b37c86fc5e 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -7,9 +7,6 @@
// Contributed by Matti Pellikka <matti.pellikka@tut.fi>.
#include "ChainComplex.h"
-#include "OS.h"
-#include "PView.h"
-#include "Options.h"
#if defined(HAVE_KBIPACK)
@@ -42,6 +39,7 @@ ChainComplex::ChainComplex(CellComplex* cellComplex, int domain)
index--;
cols--;
}
+ else _cellIndices[dim].insert( std::make_pair(index, cell));
}
index = 1;
if(dim > 0){
@@ -55,6 +53,7 @@ ChainComplex::ChainComplex(CellComplex* cellComplex, int domain)
index--;
rows--;
}
+ else _cellIndices[dim-1].insert( std::make_pair(index, cell));
}
}
@@ -434,351 +433,39 @@ std::vector<int> ChainComplex::getCoeffVector(int dim, int chainNumber)
return coeffVector;
}
-int ChainComplex::getTorsion(int dim, int chainNumber)
-{
- if(dim < 0 || dim > 4) return 0;
- if(_Hbasis[dim] == NULL
- || (int)gmp_matrix_cols(_Hbasis[dim]) < chainNumber) return 0;
- if(_torsion[dim].empty()
- || (int)_torsion[dim].size() < chainNumber) return 1;
- else return _torsion[dim].at(chainNumber-1);
-}
-
-Chain::Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs,
- CellComplex* cellComplex, GModel* model,
- std::string name, int torsion)
-{
- int i = 0;
- for(std::set<Cell*, Less_Cell>::iterator cit = cells.begin();
- cit != cells.end(); cit++){
- Cell* cell = *cit;
- _dim = cell->getDim();
- if((int)coeffs.size() > i){
- if(coeffs.at(i) != 0){
- std::list< std::pair<int, Cell*> > subCells;
- cell->getCells(subCells);
- for(std::list< std::pair<int, Cell*> >::iterator it =
- subCells.begin(); it != subCells.end(); it++){
- Cell* subCell = (*it).second;
- int coeff = (*it).first;
- _cells.insert( std::make_pair(subCell, coeffs.at(i)*coeff));
- subCell->setDeleteWithCellComplex(false);
- }
- }
- i++;
- }
-
- }
- _name = name;
- _cellComplex = cellComplex;
- _torsion = torsion;
- _model = model;
-
-}
-
-bool Chain::deform(std::map<Cell*, int, Less_Cell>& cellsInChain,
- std::map<Cell*, int, Less_Cell>& cellsNotInChain)
-{
- std::vector<int> cc;
- std::vector<int> bc;
-
- for(citer cit = cellsInChain.begin(); cit != cellsInChain.end(); cit++){
- Cell* c = (*cit).first;
- c->setImmune(false);
- if(!c->inSubdomain()) {
- cc.push_back(getCoeff(c));
- bc.push_back((*cit).second);
- }
- }
-
- if(cc.empty() || (getDim() == 2 && cc.size() < 2) ) return false;
- int inout = cc[0]*bc[0];
- for(unsigned int i = 0; i < cc.size(); i++){
- if(cc[i]*bc[i] != inout) return false;
- }
-
- for(citer cit = cellsInChain.begin(); cit != cellsInChain.end(); cit++){
- removeCell((*cit).first);
- }
-
- int n = 1;
- for(citer cit = cellsNotInChain.begin(); cit != cellsNotInChain.end();
- cit++){
- Cell* c = (*cit).first;
- if(n == 2) c->setImmune(true);
- else c->setImmune(false);
- int coeff = -1*inout*(*cit).second;
- addCell(c, coeff);
- n++;
- }
-
- return true;
-}
-
-bool Chain::deformChain(std::pair<Cell*, int> cell, bool bend)
-{
- Cell* c1 = cell.first;
- for(citer cit = c1->firstCoboundary(true); cit != c1->lastCoboundary(true);
- cit++){
-
- std::map<Cell*, int, Less_Cell> cellsInChain;
- std::map<Cell*, int, Less_Cell> cellsNotInChain;
- Cell* c1CbdCell = (*cit).first;
-
- for(citer cit2 = c1CbdCell->firstBoundary(true);
- cit2 != c1CbdCell->lastBoundary(true); cit2++){
- Cell* c1CbdBdCell = (*cit2).first;
- int coeff = (*cit2).second;
- if( (hasCell(c1CbdBdCell) && getCoeff(c1CbdBdCell) != 0)
- || c1CbdBdCell->inSubdomain()){
- cellsInChain.insert(std::make_pair(c1CbdBdCell, coeff));
- }
- else cellsNotInChain.insert(std::make_pair(c1CbdBdCell, coeff));
- }
-
- bool next = false;
-
- for(citer cit2 = cellsInChain.begin(); cit2 != cellsInChain.end(); cit2++){
- Cell* c = (*cit2).first;
- int coeff = getCoeff(c);
- if(c->getImmune()) next = true;
- if(c->inSubdomain()) bend = false;
- if(coeff > 1 || coeff < -1) next = true;
- }
-
- for(citer cit2 = cellsNotInChain.begin(); cit2 != cellsNotInChain.end();
- cit2++){
- Cell* c = (*cit2).first;
- if(c->inSubdomain()) next = true;
- }
- if(next) continue;
-
- if( (getDim() == 1 && cellsInChain.size() == 2
- && cellsNotInChain.size() == 1) ||
- (getDim() == 2 && cellsInChain.size() == 3
- && cellsNotInChain.size() == 1)){
- //printf("straighten \n");
- return deform(cellsInChain, cellsNotInChain);
- }
- else if ( (getDim() == 1 && cellsInChain.size() == 1
- && cellsNotInChain.size() == 2 && bend) ||
- (getDim() == 2 && cellsInChain.size() == 2
- && cellsNotInChain.size() == 2 && bend)){
- //printf("bend \n");
- return deform(cellsInChain, cellsNotInChain);
- }
- else if ((getDim() == 1 && cellsInChain.size() == 3
- && cellsNotInChain.size() == 0) ||
- (getDim() == 2 && cellsInChain.size() == 4
- && cellsNotInChain.size() == 0)){
- //printf("remove boundary \n");
- return deform(cellsInChain, cellsNotInChain);
- }
- }
-
- return false;
-}
-
-void Chain::smoothenChain()
+void ChainComplex::getChain(std::map<Cell*, int, Less_Cell>& chain, int dim, int chainNumber)
{
- if(!_cellComplex->simplicial()) return;
-
- int start = getSize();
- double t1 = Cpu();
-
- int useless = 0;
- for(int i = 0; i < 20; i++){
- int size = getSize();
- for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
- //if(!deformChain(*cit, false) && getDim() == 2) deformChain(*cit, true);
- if(getDim() == 2) deformChain(*cit, true);
- deformChain(*cit, false);
- }
- deImmuneCells();
- eraseNullCells();
- if (size >= getSize()) useless++;
- else useless = 0;
- if (useless > 5) break;
- }
-
- deImmuneCells();
- for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
- deformChain(*cit, false);
- }
-
- eraseNullCells();
- double t2 = Cpu();
- Msg::Debug("Smoothened a %d-chain from %d cells to %d cells (%g s).\n",
- getDim(), start, getSize(), t2-t1);
- return;
-}
-
-
-int Chain::writeChainMSH(const std::string &name)
-{
- if(getSize() == 0) return 1;
-
- FILE *fp = fopen(name.c_str(), "a");
- if(!fp){
- Msg::Error("Unable to open file '%s'", name.c_str());
- Msg::Debug("Unable to open file.");
- return 0;
- }
-
- fprintf(fp, "\n$ElementData\n");
-
- fprintf(fp, "1 \n");
- fprintf(fp, "\"%s\" \n", getName().c_str());
- fprintf(fp, "1 \n");
- fprintf(fp, "0.0 \n");
- fprintf(fp, "4 \n");
- fprintf(fp, "0 \n");
- fprintf(fp, "1 \n");
- fprintf(fp, "%d \n", getSize());
- fprintf(fp, "0 \n");
-
- 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, "$EndElementData\n");
- fclose(fp);
-
- return 1;
-}
-
-int Chain::createPGroup()
-{
- std::vector<MElement*> elements;
- std::map<int, std::vector<double> > data;
- MElementFactory factory;
-
- for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
- Cell* cell = (*cit).first;
- int coeff = (*cit).second;
-
- MElement* e = cell->getImageMElement();
- std::vector<MVertex*> v;
- for(int i = 0; i < e->getNumVertices(); i++){
- v.push_back(e->getVertex(i));
- }
- MElement* ne = factory.create(e->getTypeForMSH(), v, 0, e->getPartition());
-
- if(cell->getDim() > 0 && coeff < 0) ne->revert(); // flip orientation
- for(int i = 0; i < abs(coeff); i++) elements.push_back(ne);
-
- std::vector<double> coeffs (1,abs(coeff));
- data[ne->getNum()] = coeffs;
- }
-
- int max[4];
- for(int i = 0; i < 4; i++) max[i] = _model->getMaxElementaryNumber(i);
- int entityNum = *std::max_element(max,max+4) + 1;
- for(int i = 0; i < 4; i++) max[i] = _model->getMaxPhysicalNumber(i);
- int physicalNum = *std::max_element(max,max+4) + 1;
- setNum(physicalNum);
-
- std::map<int, std::vector<MElement*> > entityMap;
- entityMap[entityNum] = elements;
- std::map<int, std::map<int, std::string> > physicalMap;
- std::map<int, std::string> physicalInfo;
- physicalInfo[physicalNum]=getName();
- physicalMap[entityNum] = physicalInfo;
-
- // hide mesh
- opt_mesh_points(0, GMSH_SET, 0);
- opt_mesh_lines(0, GMSH_SET, 0);
- opt_mesh_triangles(0, GMSH_SET, 0);
- opt_mesh_quadrangles(0, GMSH_SET, 0);
- opt_mesh_tetrahedra(0, GMSH_SET, 0);
- opt_mesh_hexahedra(0, GMSH_SET, 0);
- opt_mesh_prisms(0, GMSH_SET, 0);
- opt_mesh_pyramids(0, GMSH_SET, 0);
-
- // show post-processing normals, tangents and element boundaries
- //opt_view_normals(0, GMSH_SET, 20);
- //opt_view_tangents(0, GMSH_SET, 20);
- opt_view_show_element(0, GMSH_SET, 1);
-
- if(!data.empty()){
- _model->storeChain(getDim(), entityMap, physicalMap);
- _model->setPhysicalName(getName(), getDim(), physicalNum);
-
- // create PView for visualization
- PView* chain = new PView(getName(), "ElementData", getGModel(),
- data, 0, 1);
- }
-
- return physicalNum;
-}
+ chain.clear();
+ if(dim < 0 || dim > 4) return;
+ if(_Hbasis[dim] == NULL
+ || (int)gmp_matrix_cols(_Hbasis[dim]) < chainNumber) return;
+ int rows = gmp_matrix_rows(_Hbasis[dim]);
-void Chain::removeCell(Cell* cell)
-{
- citer it = _cells.find(cell);
- if(it != _cells.end()){
- (*it).second = 0;
- }
- return;
-}
+ int elemi;
+ long int elemli;
+ mpz_t elem;
+ mpz_init(elem);
-void Chain::addCell(Cell* cell, int coeff)
-{
- std::pair<citer,bool> insert = _cells.insert( std::make_pair( cell, coeff));
- if(!insert.second && (*insert.first).second == 0){
- (*insert.first).second = coeff;
- cell->setDeleteWithCellComplex(false);
- }
- else if (!insert.second && (*insert.first).second != 0){
- Msg::Debug("Error: invalid chain smoothening add! \n");
+ for(int i = 1; i <= rows; i++){
+ gmp_matrix_get_elem(elem, i, chainNumber, _Hbasis[dim]);
+ elemli = mpz_get_si(elem);
+ elemi = elemli;
+ Cell* cell = _cellIndices[dim][i];
+ chain[cell] = elemi;
}
- return;
-}
-
-bool Chain::hasCell(Cell* c)
-{
- citer it = _cells.find(c);
- if(it != _cells.end() && (*it).second != 0) return true;
- return false;
-}
-
-Cell* Chain::findCell(Cell* c)
-{
- citer it = _cells.find(c);
- if(it != _cells.end() && (*it).second != 0) return (*it).first;
- return NULL;
-}
-
-int Chain::getCoeff(Cell* c)
-{
- citer it = _cells.find(c);
- if(it != _cells.end()) return (*it).second;
- return 0;
-}
-void Chain::eraseNullCells()
-{
- std::vector<Cell*> toRemove;
- for(int i = 0; i < 4; i++){
- for(citer cit = _cells.begin(); cit != _cells.end(); ++cit){
- if((*cit).second == 0) toRemove.push_back((*cit).first);
- }
- }
- for(unsigned int i = 0; i < toRemove.size(); i++){
- _cells.erase(toRemove[i]);
- toRemove[i]->setDeleteWithCellComplex(true);
- }
- return;
+ mpz_clear(elem);
}
-void Chain::deImmuneCells()
+int ChainComplex::getTorsion(int dim, int chainNumber)
{
- for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
- Cell* cell = (*cit).first;
- cell->setImmune(false);
- }
+ if(dim < 0 || dim > 4) return 0;
+ if(_Hbasis[dim] == NULL
+ || (int)gmp_matrix_cols(_Hbasis[dim]) < chainNumber) return 0;
+ if(_torsion[dim].empty()
+ || (int)_torsion[dim].size() < chainNumber) return 1;
+ else return _torsion[dim].at(chainNumber-1);
}
#endif
diff --git a/Geo/ChainComplex.h b/Geo/ChainComplex.h
index 38f6197a3461cc3ffe406b39572dde260c2a9434..fce3e0a9502cef9f506bad4509547d6552dbb518 100644
--- a/Geo/ChainComplex.h
+++ b/Geo/ChainComplex.h
@@ -17,14 +17,8 @@
#include <algorithm>
#include <set>
#include <queue>
-#include "MElement.h"
-#include "GModel.h"
-#include "GEntity.h"
-#include "GRegion.h"
-#include "GFace.h"
-#include "GVertex.h"
#include "CellComplex.h"
-
+#include "OS.h"
#if defined(HAVE_GMP)
#include "gmp.h"
@@ -60,6 +54,8 @@ class ChainComplex{
int _dim;
CellComplex* _cellComplex;
+
+ std::map<int, Cell*> _cellIndices[4];
// set the matrices
void setHMatrix(int dim, gmp_matrix* matrix) {
@@ -141,6 +137,7 @@ class ChainComplex{
// get coefficient vector for dim-dimensional chain chainNumber
// (columns of _Hbasis[dim])
std::vector<int> getCoeffVector(int dim, int chainNumber);
+ void getChain(std::map<Cell*, int, Less_Cell>& chain, int dim, int chainNumber);
// torsion coefficient for dim-dimensional chain chainNumber
int getTorsion(int dim, int chainNumber);
// get rank of homology group of dimension dim
@@ -157,97 +154,6 @@ class ChainComplex{
void matrixTest();
};
-// A class representing a chain.
-// Used to visualize generators of the homology groups.
-class Chain{
-
- private:
- // cells and their coefficients in this chain
- std::map< Cell*, int, Less_Cell > _cells;
- // name of the chain (optional)
- std::string _name;
- // physical group number of the chain
- int _num;
- // cell complex this chain belongs to
- CellComplex* _cellComplex;
- GModel* _model;
-
- // torsion coefficient
- int _torsion;
-
- int _dim;
-
- bool deform(std::map<Cell*, int, Less_Cell> &cellsInChain,
- std::map<Cell*, int, Less_Cell> &cellsNotInChain);
- bool deformChain(std::pair<Cell*, int> cell, bool bend);
-
-
- public:
- Chain() {}
- Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs,
- CellComplex* cellComplex, GModel* model,
- std::string name="Chain", int torsion=0);
- Chain(Chain* chain){
- _cells = chain->getCells();
- _torsion = chain->getTorsion();
- _name = chain->getName();
- _cellComplex = chain->getCellComplex();
- _dim = chain->getDim();
- _model = chain->getGModel();
- _num = chain->getNum();
- }
- typedef std::map<Cell*, int, Less_Cell>::iterator citer;
- citer firstCell() {return _cells.begin(); }
- citer lastCell() {return _cells.end(); }
-
- ~Chain() {}
-
- // remove a cell from this chain
- void removeCell(Cell* cell);
-
- // add a cell to this chain
- void addCell(Cell* cell, int coeff);
-
- bool hasCell(Cell* c);
- Cell* findCell(Cell* c);
- int getCoeff(Cell* c);
-
-
- int getTorsion() { return _torsion; }
- int getDim() { return _dim; }
- CellComplex* getCellComplex() { return _cellComplex; }
- GModel* getGModel() { return _model; }
- std::map<Cell*, int, Less_Cell> getCells() { return _cells; }
-
- // erase cells from the chain with zero coefficient
- void eraseNullCells();
-
- void deImmuneCells();
-
- // number of cells in this chain
- int getSize() { return _cells.size();}
-
- // get/set chain name
- std::string getName() { return _name; }
- void setName(std::string name) { _name=name; }
- // get/set physical group number
- int getNum() { return _num; }
- void setNum(int num) { _num=num; }
-
- // make local deformations to the chain to make it smoother and smaller
- // (for primary complex chains only, not for dual chains represented
- // by primary cells (yet).)
- void smoothenChain();
-
- // append this chain to a MSH ASCII file as $ElementData
- // for debugging only
- int writeChainMSH(const std::string &name);
-
- // create a physical group of this chain.
- int createPGroup();
-
-};
-
#endif
#endif
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index ee17d26ba84bade56e8069f4f94cddf2f13d17c7..85f3a9d58c69800e9cefe1e544c696b6252fd77b 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -7,8 +7,6 @@
#include "Homology.h"
-#include "ChainComplex.h"
-#include "OS.h"
#if defined(HAVE_KBIPACK)
@@ -94,17 +92,10 @@ Homology::~Homology()
}
-void Homology::findGenerators(CellComplex* cellComplex, bool omit,
- std::string domainString)
+void Homology::findGenerators()
{
- bool newComplex = false;
- if(cellComplex == NULL){
- cellComplex = createCellComplex(_domainEntities, _subdomainEntities);
- newComplex = true;
- }
- if(domainString == ""){
- domainString = getDomainString(_domain, _subdomain);
- }
+ CellComplex* cellComplex = createCellComplex(_domainEntities, _subdomainEntities);
+ std::string domainString = getDomainString(_domain, _subdomain);
Msg::Info("Reducing the Cell Complex...");
Msg::StatusBar(1, false, "Reducing...");
@@ -114,7 +105,7 @@ void Homology::findGenerators(CellComplex* cellComplex, bool omit,
cellComplex->getSize(3), cellComplex->getSize(2),
cellComplex->getSize(1), cellComplex->getSize(0));
- int omitted = cellComplex->reduceComplex(omit);
+ int omitted = cellComplex->reduceComplex();
printf(" %d volumes, %d faces, %d edges and %d vertices. \n",
cellComplex->getSize(3), cellComplex->getSize(2),
@@ -168,7 +159,7 @@ void Homology::findGenerators(CellComplex* cellComplex, bool omit,
chains->getTorsion(j,i));
t1 = Cpu();
int start = chain->getSize();
- //chain->smoothenChain();
+ chain->smoothenChain();
t2 = Cpu();
Msg::Info("Smoothened H%d %d from %d cells to %d cells (%g s).",
j, i, start, chain->getSize(), t2 - t1);
@@ -197,7 +188,7 @@ void Homology::findGenerators(CellComplex* cellComplex, bool omit,
if(_fileName != "") writeGeneratorsMSH();
- if(newComplex) delete cellComplex;
+ delete cellComplex;
delete chains;
Msg::Info("Ranks of homology spaces for primal cell complex:");
@@ -499,5 +490,355 @@ bool Homology::writeGeneratorsMSH(bool binary)
Msg::Info("Wrote homology computation results to %s.", _fileName.c_str());
return true;
}
+Chain::Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs,
+ CellComplex* cellComplex, GModel* model,
+ std::string name, int torsion)
+{
+ int i = 0;
+ for(std::set<Cell*, Less_Cell>::iterator cit = cells.begin();
+ cit != cells.end(); cit++){
+ Cell* cell = *cit;
+ _dim = cell->getDim();
+ if((int)coeffs.size() > i){
+ if(coeffs.at(i) != 0){
+ std::list< std::pair<int, Cell*> > subCells;
+ cell->getCells(subCells);
+ for(std::list< std::pair<int, Cell*> >::iterator it =
+ subCells.begin(); it != subCells.end(); it++){
+ Cell* subCell = (*it).second;
+ int coeff = (*it).first;
+ _cells.insert( std::make_pair(subCell, coeffs.at(i)*coeff));
+ subCell->setDeleteWithCellComplex(false);
+ }
+ }
+ i++;
+ }
+
+ }
+ _name = name;
+ _cellComplex = cellComplex;
+ _torsion = torsion;
+ _model = model;
+
+}
+
+Chain::Chain(std::map<Cell*, int, Less_Cell>& chain,
+ CellComplex* cellComplex, GModel* model,
+ std::string name, int torsion)
+{
+
+ _cells = chain;
+
+ _name = name;
+ _cellComplex = cellComplex;
+ _torsion = torsion;
+ _model = model;
+}
+
+bool Chain::deform(std::map<Cell*, int, Less_Cell>& cellsInChain,
+ std::map<Cell*, int, Less_Cell>& cellsNotInChain)
+{
+ std::vector<int> cc;
+ std::vector<int> bc;
+
+ for(citer cit = cellsInChain.begin(); cit != cellsInChain.end(); cit++){
+ Cell* c = (*cit).first;
+ c->setImmune(false);
+ if(!c->inSubdomain()) {
+ cc.push_back(getCoeff(c));
+ bc.push_back((*cit).second);
+ }
+ }
+
+ if(cc.empty() || (getDim() == 2 && cc.size() < 2) ) return false;
+ int inout = cc[0]*bc[0];
+ for(unsigned int i = 0; i < cc.size(); i++){
+ if(cc[i]*bc[i] != inout) return false;
+ }
+
+ for(citer cit = cellsInChain.begin(); cit != cellsInChain.end(); cit++){
+ removeCell((*cit).first);
+ }
+
+ int n = 1;
+ for(citer cit = cellsNotInChain.begin(); cit != cellsNotInChain.end();
+ cit++){
+ Cell* c = (*cit).first;
+ if(n == 2) c->setImmune(true);
+ else c->setImmune(false);
+ int coeff = -1*inout*(*cit).second;
+ addCell(c, coeff);
+ n++;
+ }
+
+ return true;
+}
+
+bool Chain::deformChain(std::pair<Cell*, int> cell, bool bend)
+{
+ Cell* c1 = cell.first;
+ for(citer cit = c1->firstCoboundary(true); cit != c1->lastCoboundary(true);
+ cit++){
+
+ std::map<Cell*, int, Less_Cell> cellsInChain;
+ std::map<Cell*, int, Less_Cell> cellsNotInChain;
+ Cell* c1CbdCell = (*cit).first;
+
+ for(citer cit2 = c1CbdCell->firstBoundary(true);
+ cit2 != c1CbdCell->lastBoundary(true); cit2++){
+ Cell* c1CbdBdCell = (*cit2).first;
+ int coeff = (*cit2).second;
+ if( (hasCell(c1CbdBdCell) && getCoeff(c1CbdBdCell) != 0)
+ || c1CbdBdCell->inSubdomain()){
+ cellsInChain.insert(std::make_pair(c1CbdBdCell, coeff));
+ }
+ else cellsNotInChain.insert(std::make_pair(c1CbdBdCell, coeff));
+ }
+
+ bool next = false;
+
+ for(citer cit2 = cellsInChain.begin(); cit2 != cellsInChain.end(); cit2++){
+ Cell* c = (*cit2).first;
+ int coeff = getCoeff(c);
+ if(c->getImmune()) next = true;
+ if(c->inSubdomain()) bend = false;
+ if(coeff > 1 || coeff < -1) next = true;
+ }
+
+ for(citer cit2 = cellsNotInChain.begin(); cit2 != cellsNotInChain.end();
+ cit2++){
+ Cell* c = (*cit2).first;
+ if(c->inSubdomain()) next = true;
+ }
+ if(next) continue;
+
+ if( (getDim() == 1 && cellsInChain.size() == 2
+ && cellsNotInChain.size() == 1) ||
+ (getDim() == 2 && cellsInChain.size() == 3
+ && cellsNotInChain.size() == 1)){
+ //printf("straighten \n");
+ return deform(cellsInChain, cellsNotInChain);
+ }
+ else if ( (getDim() == 1 && cellsInChain.size() == 1
+ && cellsNotInChain.size() == 2 && bend) ||
+ (getDim() == 2 && cellsInChain.size() == 2
+ && cellsNotInChain.size() == 2 && bend)){
+ //printf("bend \n");
+ return deform(cellsInChain, cellsNotInChain);
+ }
+ else if ((getDim() == 1 && cellsInChain.size() == 3
+ && cellsNotInChain.size() == 0) ||
+ (getDim() == 2 && cellsInChain.size() == 4
+ && cellsNotInChain.size() == 0)){
+ //printf("remove boundary \n");
+ return deform(cellsInChain, cellsNotInChain);
+ }
+ }
+
+ return false;
+}
+
+void Chain::smoothenChain()
+{
+ if(!_cellComplex->simplicial()) return;
+
+ int start = getSize();
+ double t1 = Cpu();
+
+ int useless = 0;
+ for(int i = 0; i < 20; i++){
+ int size = getSize();
+ for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
+ //if(!deformChain(*cit, false) && getDim() == 2) deformChain(*cit, true);
+ if(getDim() == 2) deformChain(*cit, true);
+ deformChain(*cit, false);
+ }
+ deImmuneCells();
+ eraseNullCells();
+ if (size >= getSize()) useless++;
+ else useless = 0;
+ if (useless > 5) break;
+ }
+
+ deImmuneCells();
+ for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
+ deformChain(*cit, false);
+ }
+
+ eraseNullCells();
+ double t2 = Cpu();
+ Msg::Debug("Smoothened a %d-chain from %d cells to %d cells (%g s).\n",
+ getDim(), start, getSize(), t2-t1);
+ return;
+}
+
+
+int Chain::writeChainMSH(const std::string &name)
+{
+ if(getSize() == 0) return 1;
+
+ FILE *fp = fopen(name.c_str(), "a");
+ if(!fp){
+ Msg::Error("Unable to open file '%s'", name.c_str());
+ Msg::Debug("Unable to open file.");
+ return 0;
+ }
+
+ fprintf(fp, "\n$ElementData\n");
+
+ fprintf(fp, "1 \n");
+ fprintf(fp, "\"%s\" \n", getName().c_str());
+ fprintf(fp, "1 \n");
+ fprintf(fp, "0.0 \n");
+ fprintf(fp, "4 \n");
+ fprintf(fp, "0 \n");
+ fprintf(fp, "1 \n");
+ fprintf(fp, "%d \n", getSize());
+ fprintf(fp, "0 \n");
+
+ 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, "$EndElementData\n");
+ fclose(fp);
+
+ return 1;
+}
+
+int Chain::createPGroup()
+{
+ std::vector<MElement*> elements;
+ std::map<int, std::vector<double> > data;
+ MElementFactory factory;
+
+ for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
+ Cell* cell = (*cit).first;
+ int coeff = (*cit).second;
+
+ MElement* e = cell->getImageMElement();
+ std::vector<MVertex*> v;
+ for(int i = 0; i < e->getNumVertices(); i++){
+ v.push_back(e->getVertex(i));
+ }
+ MElement* ne = factory.create(e->getTypeForMSH(), v, 0, e->getPartition());
+
+ if(cell->getDim() > 0 && coeff < 0) ne->revert(); // flip orientation
+ for(int i = 0; i < abs(coeff); i++) elements.push_back(ne);
+
+ std::vector<double> coeffs (1,abs(coeff));
+ data[ne->getNum()] = coeffs;
+ }
+
+ int max[4];
+ for(int i = 0; i < 4; i++) max[i] = _model->getMaxElementaryNumber(i);
+ int entityNum = *std::max_element(max,max+4) + 1;
+ for(int i = 0; i < 4; i++) max[i] = _model->getMaxPhysicalNumber(i);
+ int physicalNum = *std::max_element(max,max+4) + 1;
+ setNum(physicalNum);
+
+ std::map<int, std::vector<MElement*> > entityMap;
+ entityMap[entityNum] = elements;
+ std::map<int, std::map<int, std::string> > physicalMap;
+ std::map<int, std::string> physicalInfo;
+ physicalInfo[physicalNum]=getName();
+ physicalMap[entityNum] = physicalInfo;
+
+ // hide mesh
+ opt_mesh_points(0, GMSH_SET, 0);
+ opt_mesh_lines(0, GMSH_SET, 0);
+ opt_mesh_triangles(0, GMSH_SET, 0);
+ opt_mesh_quadrangles(0, GMSH_SET, 0);
+ opt_mesh_tetrahedra(0, GMSH_SET, 0);
+ opt_mesh_hexahedra(0, GMSH_SET, 0);
+ opt_mesh_prisms(0, GMSH_SET, 0);
+ opt_mesh_pyramids(0, GMSH_SET, 0);
+
+ // show post-processing normals, tangents and element boundaries
+ //opt_view_normals(0, GMSH_SET, 20);
+ //opt_view_tangents(0, GMSH_SET, 20);
+ opt_view_show_element(0, GMSH_SET, 1);
+
+ if(!data.empty()){
+ _model->storeChain(getDim(), entityMap, physicalMap);
+ _model->setPhysicalName(getName(), getDim(), physicalNum);
+
+ // create PView for visualization
+ PView* chain = new PView(getName(), "ElementData", getGModel(),
+ data, 0, 1);
+ }
+
+ return physicalNum;
+}
+
+
+void Chain::removeCell(Cell* cell)
+{
+ citer it = _cells.find(cell);
+ if(it != _cells.end()){
+ (*it).second = 0;
+ }
+ return;
+}
+
+void Chain::addCell(Cell* cell, int coeff)
+{
+ std::pair<citer,bool> insert = _cells.insert( std::make_pair( cell, coeff));
+ if(!insert.second && (*insert.first).second == 0){
+ (*insert.first).second = coeff;
+ cell->setDeleteWithCellComplex(false);
+ }
+ else if (!insert.second && (*insert.first).second != 0){
+ Msg::Debug("Error: invalid chain smoothening add! \n");
+ }
+ return;
+}
+
+bool Chain::hasCell(Cell* c)
+{
+ citer it = _cells.find(c);
+ if(it != _cells.end() && (*it).second != 0) return true;
+ return false;
+}
+
+Cell* Chain::findCell(Cell* c)
+{
+ citer it = _cells.find(c);
+ if(it != _cells.end() && (*it).second != 0) return (*it).first;
+ return NULL;
+}
+
+int Chain::getCoeff(Cell* c)
+{
+ citer it = _cells.find(c);
+ if(it != _cells.end()) return (*it).second;
+ return 0;
+}
+
+void Chain::eraseNullCells()
+{
+ std::vector<Cell*> toRemove;
+ for(int i = 0; i < 4; i++){
+ for(citer cit = _cells.begin(); cit != _cells.end(); ++cit){
+ if((*cit).second == 0) toRemove.push_back((*cit).first);
+ }
+ }
+ for(unsigned int i = 0; i < toRemove.size(); i++){
+ _cells.erase(toRemove[i]);
+ toRemove[i]->setDeleteWithCellComplex(true);
+ }
+ return;
+}
+
+void Chain::deImmuneCells()
+{
+ for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
+ Cell* cell = (*cit).first;
+ cell->setImmune(false);
+ }
+}
#endif
+
diff --git a/Geo/Homology.h b/Geo/Homology.h
index 62d1fd224171664d0293380bb020d3bb6a2b3445..49f80aa0aa993769c617f815dc3df86112b76121 100644
--- a/Geo/Homology.h
+++ b/Geo/Homology.h
@@ -11,6 +11,9 @@
#include <sstream>
#include "CellComplex.h"
#include "ChainComplex.h"
+#include "OS.h"
+#include "PView.h"
+#include "Options.h"
#if defined(HAVE_KBIPACK)
@@ -22,7 +25,9 @@ template <class TTypeA, class TTypeB>
return stream.good();
}
+class Chain;
+// Interface class for homology computation in Gmsh
class Homology
{
private:
@@ -57,8 +62,7 @@ class Homology
// Find the generators/duals of homology spaces,
// or just compute the ranks of homology spaces
- void findGenerators(CellComplex* cellComplex=NULL,
- bool omit=true, std::string domainString="");
+ void findGenerators();
void findDualGenerators();
void computeBettiNumbers();
@@ -77,6 +81,100 @@ class Homology
bool writeGeneratorsMSH(bool binary=false);
};
+// A class representing a chain.
+// Used to visualize generators of the homology groups.
+class Chain{
+
+ private:
+ // cells and their coefficients in this chain
+ std::map< Cell*, int, Less_Cell > _cells;
+ // name of the chain (optional)
+ std::string _name;
+ // physical group number of the chain
+ int _num;
+ // cell complex this chain belongs to
+ CellComplex* _cellComplex;
+ GModel* _model;
+
+ // torsion coefficient
+ int _torsion;
+
+ int _dim;
+
+ bool deform(std::map<Cell*, int, Less_Cell> &cellsInChain,
+ std::map<Cell*, int, Less_Cell> &cellsNotInChain);
+ bool deformChain(std::pair<Cell*, int> cell, bool bend);
+
+
+ public:
+ Chain() {}
+ Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs,
+ CellComplex* cellComplex, GModel* model,
+ std::string name="Chain", int torsion=0);
+ Chain(std::map<Cell*, int, Less_Cell>& chain,
+ CellComplex* cellComplex, GModel* model,
+ std::string name="Chain", int torsion=0);
+ Chain(Chain* chain){
+ _cells = chain->getCells();
+ _torsion = chain->getTorsion();
+ _name = chain->getName();
+ _cellComplex = chain->getCellComplex();
+ _dim = chain->getDim();
+ _model = chain->getGModel();
+ _num = chain->getNum();
+ }
+ typedef std::map<Cell*, int, Less_Cell>::iterator citer;
+ citer firstCell() {return _cells.begin(); }
+ citer lastCell() {return _cells.end(); }
+
+ ~Chain() {}
+
+ // remove a cell from this chain
+ void removeCell(Cell* cell);
+
+ // add a cell to this chain
+ void addCell(Cell* cell, int coeff);
+
+ bool hasCell(Cell* c);
+ Cell* findCell(Cell* c);
+ int getCoeff(Cell* c);
+
+
+ int getTorsion() { return _torsion; }
+ int getDim() { return _dim; }
+ CellComplex* getCellComplex() { return _cellComplex; }
+ GModel* getGModel() { return _model; }
+ std::map<Cell*, int, Less_Cell> getCells() { return _cells; }
+
+ // erase cells from the chain with zero coefficient
+ void eraseNullCells();
+
+ void deImmuneCells();
+
+ // number of cells in this chain
+ int getSize() { return _cells.size();}
+
+ // get/set chain name
+ std::string getName() { return _name; }
+ void setName(std::string name) { _name=name; }
+ // get/set physical group number
+ int getNum() { return _num; }
+ void setNum(int num) { _num=num; }
+
+ // make local deformations to the chain to make it smoother and smaller
+ // (for primary complex chains only, not for dual chains represented
+ // by primary cells (yet).)
+ void smoothenChain();
+
+ // append this chain to a MSH ASCII file as $ElementData
+ // for debugging only
+ int writeChainMSH(const std::string &name);
+
+ // create a physical group of this chain.
+ int createPGroup();
+
+};
+
#endif
#endif
diff --git a/Solver/dgRungeKuttaMultirate.cpp b/Solver/dgRungeKuttaMultirate.cpp
index ac99bd9d2777daf808ea6ae546d2b90ec969e508..ed4c6fc4a9812e7621250324dd142d6d6159599f 100644
--- a/Solver/dgRungeKuttaMultirate.cpp
+++ b/Solver/dgRungeKuttaMultirate.cpp
@@ -5,6 +5,7 @@
#include "dgDofContainer.h"
#include "dgGroupOfElements.h"
#include <algorithm>
+#include <limits.h>
dgRungeKuttaMultirate::dgRungeKuttaMultirate(dgGroupCollection* gc,dgConservationLaw *law, int nStages){
_law=law;
diff --git a/tutorial/t10.geo b/tutorial/t10.geo
index 5c32ee9b3692fae7b19c59e9447b5503b3b863de..cc9ccfa31c163d06fb4907ff85b581fbad6e5185 100644
--- a/tutorial/t10.geo
+++ b/tutorial/t10.geo
@@ -144,10 +144,6 @@ HomGen("t10_hom.msh") = {{69}, {75}};
// Save the cut chains to t10_hom.msh.
HomCut("t10_hom.msh") = {{69}, {70, 71, 72, 73}};
-// Only find and save the ranks of the relative homology spaces
-// (Betti numbers) to t10_homrank.txt. Does not find generators.
-HomRank("t10_hom.txt") = {{69},{70, 71, 72, 73}};
-
// More examples (uncomment):
// HomGen("t10_hom.msh") = {{69}, {}};
// HomGen("t10_hom.msh") = {{}, {}};
diff --git a/utils/api_demos/mainHomology.cpp b/utils/api_demos/mainHomology.cpp
index cc1fcbe01b0289638a3c292a0a4229dab7d3d7d4..6b6aecc25229581eee025548bdbd4faf0ccbfa3c 100644
--- a/utils/api_demos/mainHomology.cpp
+++ b/utils/api_demos/mainHomology.cpp
@@ -31,18 +31,15 @@ int main(int argc, char **argv)
Homology* homology1 = new Homology(m, domain, subdomain);
std::string fileName = "homology.msh";
- homology->findGenerators(fileName);
+ homology->findGenerators();
+ homology->setFileName(fileName);
delete homology1;
Homology* homology2 = new Homology(m, domain, subdomain);
- fileName = "homology.msh";
- homology->findDualGenerators(fileName);
+ homology->findDualGenerators();
+ homology->setFileName(fileName);
delete homology2;
- Homology* homology3 = new Homology(m, domain, subdomain);
- homology->computeBettiNumbers();
- delete homology3;
-
delete m;
GmshFinalize();