From aa3c5c65603d37034a145408de7ed6b5949b33fe Mon Sep 17 00:00:00 2001
From: Matti Pellika <matti.pellikka@tut.fi>
Date: Tue, 31 Mar 2009 14:22:02 +0000
Subject: [PATCH] Bugfixes in CellComplex and in the examples.
Begun writing ChainComplex class. Will have own file some day.
---
Geo/CellComplex.cpp | 111 ++++++++++++++++++++++++++++++++-----
Geo/CellComplex.h | 79 +++++++++++++++++++-------
utils/misc/celldriver.cpp | 9 ++-
utils/misc/transformer.geo | 14 ++---
4 files changed, 167 insertions(+), 46 deletions(-)
diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 25b8cb1285..2980645832 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -3,7 +3,7 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
//
-// Contributed by Matti Pellikka, 16.3.2009.
+// Contributed by Matti Pellikka.
#include "CellComplex.h"
@@ -18,7 +18,6 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su
insertCells(true);
insertCells(false);
-
}
void CellComplex::insertCells(bool subdomain){
@@ -79,7 +78,7 @@ int Simplex::kappa(Cell* tau) const{
value = value*-1;
}
- return 1;
+ return value;
}
std::set<Cell*, Less_Cell>::iterator CellComplex::findCell(int dim, std::vector<int>& vertices){
if(dim == 0) return _cells[dim].find(new ZeroSimplex(vertices.at(0)));
@@ -289,7 +288,6 @@ int CellComplex::coreduction(int dim){
bd_c = bd(cell);
if(bd_c.size() == 1 && inSameDomain(cell, bd_c.at(0)) ){
removeCell(cell);
- //cit = firstCell(dim-1);
removeCell(bd_c.at(0));
count++;
coreduced =true;
@@ -314,11 +312,9 @@ int CellComplex::reduction(int dim){
reduced = false;
for(citer cit = firstCell(dim-1); cit != lastCell(dim-1); cit++){
Cell* cell = *cit;
- //if(reductionMrozek(cell) !=0) cit = firstCell(dim-1);
cbd_c = cbd(cell);
if(cbd_c.size() == 1 && inSameDomain(cell, cbd_c.at(0)) ){
removeCell(cell);
- //cit = firstCell(dim-1);
removeCell(cbd_c.at(0));
count++;
reduced =true;
@@ -355,49 +351,117 @@ void CellComplex::coreduceComplex(){
printf("Cell complex after coreduction: %d volumes, %d faces, %d edges and %d vertices.\n",
getSize(3), getSize(2), getSize(1), getSize(0));
}
-
+/*
void CellComplex::constructHMatrix(int dim){
- // h: C_dim -> C_(dim-1)
+ // h_dim: C_dim -> C_(dim-1)
+
if(dim > 3 || dim < 1){
return;
}
+ destroy_gmp_matrix(_HMatrix[dim]);
+
if( _cells[dim].size() == 0 ){
_HMatrix[dim] = create_gmp_matrix_zero(1, 1);
+ //gmp_matrix_printf(_HMatrix[dim]);
return;
}
unsigned int cols = _cells[dim].size();
if( _cells[dim-1].size() == 0){
_HMatrix[dim] = create_gmp_matrix_zero(1, cols);
+ //gmp_matrix_printf(_HMatrix[dim]);
return;
}
unsigned int rows = _cells[dim-1].size();
mpz_t elems[rows*cols];
- mpz_array_init( elems[0], rows*cols, 512 );
-
+ mpz_array_init( elems[0], rows*cols, sizeof(mpz_t) );
+
citer high = firstCell(dim);
citer low = firstCell(dim-1);
- //printf("laa %d %d %d \n", rows, cols, rows*cols);
+ /*
+ printf("laa %d %d %d \n", rows, cols, rows*cols);
for(unsigned int i=0; i < rows*cols; i++){
- //printf("%d, ", i);
+ printf(" %d, ", i);
if(low == lastCell(dim-1)){
- //printf("rowfull %d", i);
+ printf("rowfull %d ", i);
high++;
low = firstCell(dim-1);
}
mpz_set_ui(elems[i], kappa(*high, *low));
low++;
}
+ */
+ /*
+ unsigned int i = 0;
+ while(i < rows*cols){
+ while(low != lastCell(dim-1)){
+ mpz_set_si(elems[i], kappa(*high, *low));
+ //printf(" %d %d, ",i, kappa(*high, *low));
+ i++;
+ low++;
+ }
+ low = firstCell(dim-1);
+ high++;
+ }
+
+ _HMatrix[dim] = create_gmp_matrix(rows, cols,(const mpz_t *) elems);
+
+ //gmp_matrix_printf(_HMatrix[dim]);
- _HMatrix[dim] = create_gmp_matrix(rows, cols, elems);
-
return;
}
+*/
+std::vector<gmp_matrix*> CellComplex::constructHMatrices(){
+
+ // h_dim: C_dim -> C_(dim-1)
+
+ std::vector<gmp_matrix*> HMatrix;
+
+ HMatrix.push_back(create_gmp_matrix_zero(1, 1));
+
+ for(int dim = 1; dim < 4; dim++){
+ if( _cells[dim].size() == 0 ){
+ HMatrix.push_back( create_gmp_matrix_zero(1, 1));
+ //gmp_matrix_printf(HMatrix[dim]);
+ break;
+ }
+ unsigned int cols = _cells[dim].size();
+
+ if( _cells[dim-1].size() == 0){
+ HMatrix.push_back(create_gmp_matrix_zero(1, cols));
+ //gmp_matrix_printf(HMatrix[dim]);
+ break;
+ }
+ unsigned int rows = _cells[dim-1].size();
+ mpz_t elems[rows*cols];
+ mpz_array_init( elems[0], rows*cols, sizeof(mpz_t) );
+
+ citer high = firstCell(dim);
+ citer low = firstCell(dim-1);
+
+
+ unsigned int i = 0;
+ while(i < rows*cols){
+ while(low != lastCell(dim-1)){
+ mpz_set_si(elems[i], kappa(*high, *low));
+ //printf(" %d %d, ",i, kappa(*high, *low));
+ i++;
+ low++;
+ }
+ low = firstCell(dim-1);
+ high++;
+ }
+ HMatrix.push_back(create_gmp_matrix(rows, cols,(const mpz_t *) elems));
+
+ //gmp_matrix_printf(_HMatrix[dim]);
+ }
+ return HMatrix;
+}
void CellComplex::getDomainVertices(std::set<MVertex*, Less_MVertex>& domainVertices, bool subdomain){
@@ -521,4 +585,21 @@ void CellComplex::printComplex(int dim){
}
}
+gmp_matrix* ChainComplex::ker(gmp_matrix* HMatrix){
+
+ // H = USV -> WH = US, W = inv(V)
+ gmp_normal_form* W_USform = create_gmp_Smith_normal_form(HMatrix, NOT_INVERTED, INVERTED);
+
+
+
+ return W_USform->canonical;
+
+}
+
+
+
+
+
#endif
+
+
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index 6670458a88..b7f4dcc619 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -3,7 +3,7 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
//
-// Contributed by Matti Pellikka, 16.3.2009.
+// Contributed by Matti Pellikka.
#ifndef _CELLCOMPLEX_H_
#define _CELLCOMPLEX_H_
@@ -155,6 +155,7 @@ class ZeroSimplex : public Simplex
virtual ~ZeroSimplex(){}
+ virtual int getDim() const { return 0; }
virtual int getNumVertices() const { return 1; }
virtual int getVertex(int vertex) const {return _v; }
virtual bool hasVertex(int vertex) const {return (_v == vertex); }
@@ -198,6 +199,7 @@ class OneSimplex : public Simplex
virtual ~OneSimplex(){}
+ virtual int getDim() const { return 1; }
virtual int getNumVertices() const { return 2; }
virtual int getVertex(int vertex) const {return _v[vertex]; }
virtual bool hasVertex(int vertex) const {return (_v[0] == vertex || _v[1] == vertex); }
@@ -236,6 +238,7 @@ class TwoSimplex : public Simplex
virtual ~TwoSimplex(){}
+ virtual int getDim() const { return 2; }
virtual int getNumVertices() const { return 3; }
virtual int getVertex(int vertex) const {return _v[vertex]; }
virtual bool hasVertex(int vertex) const {return
@@ -276,6 +279,7 @@ class ThreeSimplex : public Simplex
virtual ~ThreeSimplex(){}
+ virtual int getDim() const { return 3; }
virtual int getNumVertices() const { return 4; }
virtual int getVertex(int vertex) const {return _v[vertex]; }
virtual bool hasVertex(int vertex) const {return
@@ -294,18 +298,14 @@ class Less_Cell{
return (c1->getNumVertices() < c2->getNumVertices());
}
for(int i=0; i < c1->getNumVertices();i++){
- if(c1->getVertex(i) < c2->getVertex(i)){
- return true;
- }
- else if (c1->getVertex(i) > c2->getVertex(i)){
- return false;
- }
+ if(c1->getVertex(i) < c2->getVertex(i)) return true;
+ else if (c1->getVertex(i) > c2->getVertex(i)) return false;
}
return false;
}
};
-// Ordering for the finite element mesh vertices
+// Ordering for the finite element mesh vertices.
class Less_MVertex{
public:
bool operator()(const MVertex* v1, const MVertex* v2) const {
@@ -324,14 +324,13 @@ class CellComplex
// used in relative homology computation, may be empty
std::vector<GEntity*> _subdomain;
- // sorted container of unique cells in this cell complex
- // mapped according to their dimension
- std::map< int, std::set<Cell*, Less_Cell> > _cells;
-
- // boundary operator matrices for this cell complex
+ // sorted containers of unique cells in this cell complex
+ // one for each dimension
+ std::set<Cell*, Less_Cell> _cells[4];
+
+ // boundary operator matrices for this chain complex
// h_k: C_k -> C_(k-1)
- std::map<int, gmp_matrix*> _HMatrix;
-
+ //gmp_matrix* _HMatrix[4];
public:
// iterator for the cells of same dimension
@@ -355,6 +354,11 @@ class CellComplex
virtual void insertCells(bool subdomain);
public:
+ CellComplex( std::set<Cell*, Less_Cell>* cells ) {
+ for(int i = 0; i < 4; i++){
+ _cells[i] = cells[i];
+ }
+ }
CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> subdomain );
virtual ~CellComplex(){}
@@ -404,20 +408,53 @@ class CellComplex
// slower, but produces cleaner result
virtual int coreductionMrozek(Cell* generator);
-
- // construct boundary operator matrix dim->dim-1
- virtual void constructHMatrix(int dim);
- // get the boundary operator matrix dim->dim-1
- virtual gmp_matrix* getHMatrix(int dim) { return _HMatrix[dim]; }
-
// print the vertices of cells of certain dimension
virtual void printComplex(int dim);
// write this cell complex in legacy .msh format
virtual int writeComplexMSH(const std::string &name);
+ // construct boundary operator matrix dim->dim-1
+ //virtual void constructHMatrix(int dim);
+ // get the boundary operator matrix dim->dim-1
+ //virtual gmp_matrix* getHMatrix(int dim) { return _HMatrix[dim]; }
+
+ virtual std::vector<gmp_matrix*> constructHMatrices();
+
+};
+
+class ChainComplex{
+ private:
+ // boundary operator matrices for this chain complex
+ // h_k: C_k -> C_(k-1)
+ gmp_matrix* _HMatrix[4];
+
+ public:
+
+ ChainComplex( std::vector<gmp_matrix*> HMatrix ){
+ for(int i = 0; i < HMatrix.size(); i++){
+ _HMatrix[i] = HMatrix.at(i);
+ }
+ }
+ ChainComplex(){
+ for(int i = 0; i < 4; i++){
+ _HMatrix[i] = create_gmp_matrix_zero(1,1);
+ }
+ }
+ virtual ~ChainComplex(){}
+
+ // get the boundary operator matrix dim->dim-1
+ virtual gmp_matrix* getHMatrix(int dim) { return _HMatrix[dim]; }
+
+ virtual gmp_matrix* ker(gmp_matrix* HMatrix);
+
+ 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); }
+
};
+
#endif
#endif
diff --git a/utils/misc/celldriver.cpp b/utils/misc/celldriver.cpp
index c8b91e03a9..278eff7a47 100755
--- a/utils/misc/celldriver.cpp
+++ b/utils/misc/celldriver.cpp
@@ -1,10 +1,13 @@
// configure, compile and install Gmsh as a library with
//
-// ./configure --disable-gui --disable-netgen --disable-chaco
-// --disable-metis --disable-tetgen --prefix=/usr/local/
+// ./configure --disable-gui --disable-netgen --disable-chaco --prefix=/usr/local/
// make install-lib
//
-// Then compile this driver with "g++ kaka.cpp -lGmsh -llapack -lblas -lgmp"
+// 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.
+
#include <stdio.h>
#include <gmsh/Gmsh.h>
diff --git a/utils/misc/transformer.geo b/utils/misc/transformer.geo
index 9021d10d66..70212b799c 100755
--- a/utils/misc/transformer.geo
+++ b/utils/misc/transformer.geo
@@ -89,11 +89,11 @@ Line Loop(57) = {34, -39, -38, -44};
Ruled Surface(58) = {57};
Line Loop(59) = {33, 44, -37, -42};
Ruled Surface(60) = {59};
-Line Loop(61) = {28, -1, 25, 16, 29, -11, -32, 8};
-Line Loop(62) = {33, 34, 35, 36};
-Ruled Surface(63) = {61, 62};
-Line Loop(64) = {3, -26, -14, -30, 9, 31, -6, -27};
-Line Loop(65) = {37, 38, -40, 41};
-Ruled Surface(66) = {64, 65};
-Surface Loop(67) = {63, 46, 66, 18, 48, 24, 50, 22, 52, 20, 54, 60, 58, 56};
+Line Loop(61) = {27, 6, -31, -9, 30, 14, 26, -3};
+Line Loop(62) = {37, 38, -40, 41};
+Plane Surface(63) = {61, 62};
+Line Loop(64) = {25, 16, 29, -11, -32, 8, 28, -1};
+Line Loop(65) = {34, 35, 36, 33};
+Plane Surface(66) = {64, 65};
+Surface Loop(67) = {66, 48, 24, 63, 46, 20, 52, 22, 50, 18, 54, 60, 58, 56};
Volume(68) = {67};
--
GitLab