From d7f6bb25c579a14436d01da96b64af1fce320a48 Mon Sep 17 00:00:00 2001
From: Matti Pellika <matti.pellikka@tut.fi>
Date: Wed, 22 Apr 2009 07:07:53 +0000
Subject: [PATCH] Added Chain class in order to visualize the homology
generator chains. Modified utils/misc/celldriver.cpp example to demonstrate
this.
---
Geo/CellComplex.cpp | 38 ++++++++++------
Geo/CellComplex.h | 23 +++++-----
Geo/ChainComplex.cpp | 92 +++++++++++++++++++++++++++++++++++++--
Geo/ChainComplex.h | 40 ++++++++++++++++-
utils/misc/celldriver.cpp | 56 ++++++++++++++++++++----
5 files changed, 212 insertions(+), 37 deletions(-)
diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 5d17436d71..b7da2fc688 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -15,6 +15,15 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su
// subdomain need to be inserted first!
insertCells(true);
insertCells(false);
+
+ 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++;
+ }
+ }
}
void CellComplex::insertCells(bool subdomain){
@@ -493,13 +502,18 @@ 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());
printf("Unable to open file.");
return 0;
}
- fprintf(fp, "$NOD\n");
+
+
+ fprintf(fp, "$MeshFormat\n2.0 0 8\n$EndMeshFormat\n");
+
+ fprintf(fp, "$Nodes\n");
std::set<MVertex*, Less_MVertex> domainVertices;
getDomainVertices(domainVertices, true);
@@ -513,38 +527,34 @@ int CellComplex::writeComplexMSH(const std::string &name){
}
- fprintf(fp, "$ENDNOD\n");
- fprintf(fp, "$ELM\n");
+ fprintf(fp, "$EndNodes\n");
+ fprintf(fp, "$Elements\n");
fprintf(fp, "%d\n", _cells[0].size() + _cells[1].size() + _cells[2].size() + _cells[3].size());
- int index = 1;
-
for(citer cit = firstCell(0); cit != lastCell(0); cit++) {
Cell* vertex = *cit;
- fprintf(fp, "%d %d %d %d %d %d\n", index, 15, 0, 1, 1, vertex->getVertex(0));
- index++;
+ fprintf(fp, "%d %d %d %d %d %d %d\n", vertex->getTag(), 15, 3, 0, 0, 0, vertex->getVertex(0));
}
for(citer cit = firstCell(1); cit != lastCell(1); cit++) {
Cell* edge = *cit;
- fprintf(fp, "%d %d %d %d %d %d %d\n", index, 1, 0, 1, 2, edge->getVertex(0), edge->getVertex(1));
- index++;
+ fprintf(fp, "%d %d %d %d %d %d %d %d\n", edge->getTag(), 1, 3, 0, 0, 0, edge->getVertex(0), edge->getVertex(1));
}
for(citer cit = firstCell(2); cit != lastCell(2); cit++) {
Cell* face = *cit;
- fprintf(fp, "%d %d %d %d %d %d %d %d\n", index, 2, 0, 1, 3, face->getVertex(0), face->getVertex(1), face->getVertex(2));
- index++;
+ fprintf(fp, "%d %d %d %d %d %d %d %d %d\n", face->getTag(), 2, 3, 0, 0, 0, face->getVertex(0), face->getVertex(1), face->getVertex(2));
}
for(citer cit = firstCell(3); cit != lastCell(3); cit++) {
Cell* volume = *cit;
- fprintf(fp, "%d %d %d %d %d %d %d %d %d\n", index, 4, 0, 1, 4, volume->getVertex(0), volume->getVertex(1), volume->getVertex(2), volume->getVertex(3));
- index++;
+ fprintf(fp, "%d %d %d %d %d %d %d %d %d %d\n", volume->getTag(), 4, 3, 0, 0, 0, volume->getVertex(0), volume->getVertex(1), volume->getVertex(2), volume->getVertex(3));
}
- fprintf(fp, "$ENDELM\n");
+ fprintf(fp, "$EndElements\n");
+
+ fclose(fp);
return 1;
}
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index 9cf34a2a59..ee70c278c8 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -35,17 +35,20 @@ class Cell
int _bdSize;
int _cbdSize;
-
+
// whether this cell belongs to a subdomain
// used in relative homology computation
bool _inSubdomain;
+ int _tag;
+
public:
Cell(){}
virtual ~Cell(){}
virtual int getDim() const { return _dim; };
- //virtual int getTag() const { return _tag; };
+ virtual int getTag() const { return _tag; };
+ virtual void setTag(int tag) { _tag = tag; };
// get the number of vertices this cell has
virtual int getNumVertices() const = 0; //{return _vertices.size();}
@@ -68,7 +71,7 @@ class Cell
virtual void setBdSize(int size) { _bdSize = size; }
virtual int getCbdSize() { return _cbdSize; }
virtual void setCbdSize(int size) { _cbdSize = size; }
-
+
virtual bool inSubdomain() const { return _inSubdomain; }
// print cell vertices
@@ -112,7 +115,7 @@ class Simplex : public Cell
{
public:
- Simplex(){
+ Simplex() : Cell() {
_cbdSize = 1000; // big number
_bdSize = 1000;
}
@@ -344,11 +347,7 @@ class CellComplex
// remove cell from the queue set
virtual void removeCellQset(Cell*& cell, std::set<Cell*, Less_Cell>& Qset);
-
- // get all the finite element mesh vertices in the domain of this cell complex
- virtual void getDomainVertices(std::set<MVertex*, Less_MVertex>& domainVertices,
- bool subdomain);
-
+
// insert cells into this cell complex
virtual void insertCells(bool subdomain);
@@ -369,6 +368,9 @@ class CellComplex
virtual std::set<Cell*, Less_Cell> getCells(int dim){ return _cells[dim]; }
+ // get all the finite element mesh vertices in the domain of this cell complex
+ virtual void getDomainVertices(std::set<MVertex*, Less_MVertex>& domainVertices, bool subdomain);
+
// iterators to the first and last cells of certain dimension
virtual citer firstCell(int dim) {return _cells[dim].begin(); }
virtual citer lastCell(int dim) {return _cells[dim].end(); }
@@ -377,7 +379,6 @@ class CellComplex
virtual std::set<Cell*, Less_Cell>::iterator findCell(int dim, std::vector<int>& vertices);
virtual std::set<Cell*, Less_Cell>::iterator findCell(int dim, int vertex, int dummy=0);
-
// kappa for two cells of this cell complex
// implementation will vary depending on cell type
virtual int kappa(Cell* sigma, Cell* tau) const { return sigma->kappa(tau); }
@@ -424,7 +425,7 @@ class CellComplex
// print the vertices of cells of certain dimension
virtual void printComplex(int dim, bool subcomplex);
- // write this cell complex in legacy .msh format
+ // write this cell complex in .msh format
virtual int writeComplexMSH(const std::string &name);
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index 5bdcd8f517..7bb70d2e6a 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -12,6 +12,11 @@
ChainComplex::ChainComplex(CellComplex* cellComplex){
_HMatrix[0] = NULL;
+ _kerH[0] = NULL;
+ _codH[0] = NULL;
+ _JMatrix[0] = NULL;
+ _QMatrix[0] = NULL;
+ _Hbasis[0] = NULL;
for(int dim = 1; dim < 4; dim++){
unsigned int cols = cellComplex->getSize(dim);
@@ -256,15 +261,17 @@ void ChainComplex::computeHomology(){
if(getCodHMatrix(i+1) == NULL){
_Hbasis[i] = copy_gmp_matrix(getKerHMatrix(i), 1, 1,
gmp_matrix_rows(getKerHMatrix(i)), gmp_matrix_cols(getKerHMatrix(i)));
- }
+ }
else if(getJMatrix(i) == NULL || getQMatrix(i) == NULL){
_Hbasis[i] = NULL;
- }
+ }
else{
_Hbasis[i] = copy_gmp_matrix(getKerHMatrix(i), 1, 1,
gmp_matrix_rows(getKerHMatrix(i)), gmp_matrix_cols(getKerHMatrix(i)));
+
gmp_matrix_right_mult(_Hbasis[i], getQMatrix(i));
- }
+ }
+
}
@@ -303,5 +310,84 @@ void ChainComplex::matrixTest(){
return;
}
+std::vector<int> ChainComplex::getCoeffVector(int dim, int chainNumber){
+
+ std::vector<int> coeffVector;
+
+ if(dim < 0 || dim > 3) return coeffVector;
+ if(_Hbasis[dim] == NULL || gmp_matrix_cols(_Hbasis[dim]) < chainNumber) return coeffVector;
+
+ int rows = gmp_matrix_rows(_Hbasis[dim]);
+
+ int elemi;
+ long int elemli;
+ mpz_t elem;
+ mpz_init(elem);
+
+ for(int i = 1; i <= rows; i++){
+ gmp_matrix_get_elem(elem, i, chainNumber, _Hbasis[dim]);
+ elemli = mpz_get_si(elem);
+ elemi = elemli;
+ coeffVector.push_back(elemi);
+ //printf("coeff: %d \n", coeffVector.at(i-1));
+ }
+
+ mpz_clear(elem);
+ return coeffVector;
+
+}
#endif
+
+Chain::Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs, CellComplex* cellComplex, std::string name){
+
+ int i = 0;
+ for(std::set<Cell*, Less_Cell>::iterator cit = cells.begin(); cit != cells.end(); cit++){
+ Cell* cell = *cit;
+ if(!cell->inSubdomain() && coeffs.size() > i){
+ if(coeffs.at(i) != 0) _cells.push_back( std::make_pair(cell, coeffs.at(i)) );
+ i++;
+ }
+ }
+
+ _name = name;
+ _cellComplex = cellComplex;
+
+}
+
+int Chain::writeChainMSH(const std::string &name){
+
+ //_cellComplex->writeComplexMSH(name);
+
+ FILE *fp = fopen(name.c_str(), "a");
+ if(!fp){
+ Msg::Error("Unable to open file '%s'", name.c_str());
+ printf("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(int i = 0; i < _cells.size(); i++){
+
+ fprintf(fp, "%d %d \n", getCell(i)->getTag(), getCoeff(i));
+
+ }
+
+ fprintf(fp, "$EndElementData\n");
+
+ fclose(fp);
+
+ return 1;
+
+}
diff --git a/Geo/ChainComplex.h b/Geo/ChainComplex.h
index 89b505a900..4e37193107 100644
--- a/Geo/ChainComplex.h
+++ b/Geo/ChainComplex.h
@@ -49,6 +49,8 @@ class ChainComplex{
// bases for homology groups
gmp_matrix* _Hbasis[4];
+
+
public:
ChainComplex(CellComplex* cellComplex);
@@ -85,6 +87,9 @@ class ChainComplex{
// Compute bases for the homology groups of this chain complex
virtual void computeHomology();
+ virtual std::vector<int> getCoeffVector(int dim, int chainNumber);
+ virtual int getBasisSize(int dim) { if(dim > -1 || dim < 4) return gmp_matrix_cols(_Hbasis[dim]); else return 0; }
+
virtual int printMatrix(gmp_matrix* matrix){
printf("%d rows and %d columns\n", gmp_matrix_rows(matrix), gmp_matrix_cols(matrix));
return gmp_matrix_printf(matrix); }
@@ -93,7 +98,40 @@ class ChainComplex{
virtual void matrixTest();
};
-
#endif
+// A class representing a chain.
+// Used to visualize generators of the homology groups.
+class Chain{
+
+ private:
+ // cells and their coefficients in this chain
+ std::vector< std::pair <Cell*, int> > _cells;
+ // name of the chain (optional)
+ std::string _name;
+ // cell complex this chain belongs to
+ CellComplex* _cellComplex;
+
+ public:
+ Chain(){}
+ Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs, CellComplex* cellComplex, std::string name="Chain");
+ ~Chain(){}
+
+ // get i:th cell of this chain
+ virtual Cell* getCell(int i) { return _cells.at(i).first; }
+ // get coeffcient of i:th cell of this chain
+ virtual int getCoeff(int i) { return _cells.at(i).second; }
+
+ // number of cells in this chain
+ virtual int getSize() { return _cells.size(); }
+
+ // get/set chain name
+ virtual std::string getName() { return _name; }
+ virtual void setName(std::string name) { _name=name; }
+
+ // append this chain to a 2.0 .msh file as $ElementData
+ virtual int writeChainMSH(const std::string &name);
+
+};
+
#endif
diff --git a/utils/misc/celldriver.cpp b/utils/misc/celldriver.cpp
index 278eff7a47..a515b2c819 100755
--- a/utils/misc/celldriver.cpp
+++ b/utils/misc/celldriver.cpp
@@ -5,15 +5,26 @@
//
// Then compile this driver with "g++ celldriver.cpp -lGmsh -llapack -lblas -lgmp"
//
-// This program creates .msh files reduced_complex.msh and coreduced_complex.msh
-// of an two port transformer model to represent its relative homology groups.
+//
+// This program creates .msh file chains.msh which represents the generators of
+// an two port transformer model's homology groups.
#include <stdio.h>
+#include <sstream>
#include <gmsh/Gmsh.h>
#include <gmsh/GModel.h>
#include <gmsh/MElement.h>
#include <gmsh/CellComplex.h>
+#include <gmsh/ChainComplex.h>
+
+template <class TTypeA, class TTypeB>
+bool convert(const TTypeA& input, TTypeB& output ){
+ std::stringstream stream;
+ stream << input;
+ stream >> output;
+ return stream.good();
+}
int main(int argc, char **argv)
@@ -46,16 +57,45 @@ int main(int argc, char **argv)
s= *(++sub);
subdomain.push_back(s);
- CellComplex complex = CellComplex(domain, subdomain);
+ // create a cell complex
+ CellComplex* complex = new CellComplex(domain, subdomain);
printf("Cell complex of this model has: %d volumes, %d faces, %d edges and %d vertices\n",
- complex.getSize(3), complex.getSize(2), complex.getSize(1), complex.getSize(0));
+ complex->getSize(3), complex->getSize(2), complex->getSize(1), complex->getSize(0));
+
+ // reduce the complex in order to ease homology computation
+ complex->reduceComplex();
+
+ // create a chain complex of a cell complex (construct boundary operator matrices)
+ ChainComplex* chains = new ChainComplex(complex);
+ // compute the homology using the boundary operator matrices
+ chains->computeHomology();
+
+ // write the reduced cell complex to a .msh file
+ complex->writeComplexMSH("chains.msh");
+
+ // Append the homology generators to the .msh file as post-processing views.
+ // To visualize the result, open chains.msh with Gmsh GUI and deselect all mesh
+ // entities from Tools->Options->Mesh->Visibility.
+ for(int j = 0; j < 4; j++){
+ for(int i = 1; i <= chains->getBasisSize(j); i++){
+
+ std::string generator;
+ std::string dimension;
+ convert(i, generator);
+ convert(j, dimension);
+ std::string name = dimension + "D Generator " + generator;
+
+ Chain* chain = new Chain(complex->getCells(j), chains->getCoeffVector(j,i), complex, name);
+ chain->writeChainMSH("chains.msh");
+ delete chain;
+
+ }
+ }
- complex.reduceComplex();
- complex.writeComplexMSH("reduced_complex.msh");
- complex.coreduceComplex();
- complex.writeComplexMSH("coreduced_complex.msh");
+ delete chains;
+ delete complex;
delete m;
GmshFinalize();
--
GitLab