From 1a5bd0372eb7d0b5c1da71541ad7e867d201d3c9 Mon Sep 17 00:00:00 2001
From: Matti Pellika <matti.pellikka@tut.fi>
Date: Thu, 5 Nov 2009 08:53:30 +0000
Subject: [PATCH] Faster 2-chain smoothening. Updated api demo.
---
Geo/CellComplex.cpp | 2 +-
Geo/ChainComplex.cpp | 129 ++++++++++++++------------
Geo/ChainComplex.h | 2 +-
Geo/Homology.cpp | 2 +-
Geo/Homology.h | 12 ++-
utils/api_demos/mainHomology.cpp | 150 ++++++-------------------------
6 files changed, 111 insertions(+), 186 deletions(-)
diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 09876bbe70..01164cc1e4 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.
+// Contributed by Matti Pellikka <matti.pellikka@tut.fi>.
#include "CellComplex.h"
#include "OS.h"
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index ed1fa49b20..f25ec4c13d 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -783,8 +783,24 @@ void Chain::smoothenChain(){
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(!straightenChain(*cit) && getDim() == 2) bendChain(*cit);
+ removeBoundary(*cit);
+ }
+ deImmuneCells();
+ eraseNullCells();
+ if (size >= getSize()) useless++;
+ else useless = 0;
+ if (useless > 5) break;
+ }
+
+
+ /*
const int MAXROUNDS = 2;
-
bool smoothened = true;
int round = 0;
while(smoothened){
@@ -798,67 +814,70 @@ void Chain::smoothenChain(){
}
eraseNullCells();
- if(getDim() == 2){
- Chain* bestChain = new Chain(this);
- int before = getSize();
- deImmuneCells();
- srand ( time(NULL) );
- int n = 0;
- while( n < 10){
- double t = 1;
- while(t>0){
- double tt = 1;
-
- for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
- int random = rand() % 100 + 1;
- double r = random*t;
- //printf("random: %d, t: %d \n", random, t);
- if(r > 80) {
- bendChain(*cit);
- //printf("random: %d, t: %g, random*t: %g.\n", random, t, r);
- //cit = _cells.begin();
- //tt = tt - 0.1;
- //if(tt <= 0) tt = 0;
- }
- }
- eraseNullCells();
-
- smoothened = true;
- round = 0;
- while(smoothened){
- round++;
- smoothened = false;
+ if(getDim() == 2){
+ Chain* bestChain = new Chain(this);
+ int before = getSize();
+ deImmuneCells();
+ srand ( time(NULL) );
+ int n = 0;
+ while( n < 10){
+ double t = 1;
+ while(t>0){
+ double tt = 1;
+
for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
- if(straightenChain(*cit)){smoothened = true;}
- if(removeBoundary(*cit)) { smoothened = true; }
+ int random = rand() % 100 + 1;
+ double r = random*t;
+ //printf("random: %d, t: %d \n", random, t);
+ if(r > 80) {
+ bendChain(*cit);
+ //printf("random: %d, t: %g, random*t: %g.\n", random, t, r);
+ //cit = _cells.begin();
+ //tt = tt - 0.1;
+ //if(tt <= 0) tt = 0;
+ }
}
- if(round > MAXROUNDS) smoothened = false;
+ eraseNullCells();
+
+ smoothened = true;
+ round = 0;
+ while(smoothened){
+ round++;
+ smoothened = false;
+ for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
+ if(straightenChain(*cit)){smoothened = true;}
+ if(removeBoundary(*cit)) { smoothened = true; }
+ }
+ if(round > MAXROUNDS) smoothened = false;
+ }
+ eraseNullCells();
+ if(this->getSize() < bestChain->getSize()) bestChain = this;
+ t = t - 0.1;
}
- eraseNullCells();
- if(this->getSize() < bestChain->getSize()) bestChain = this;
- t = t - 0.1;
+ deImmuneCells();
+ n=n+1;
}
+
+
deImmuneCells();
- n=n+1;
- }
-
-
- deImmuneCells();
- smoothened = true;
- round = 0;
- while(smoothened){
- round++;
- smoothened = false;
- for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
- if(straightenChain(*cit)){smoothened = true;}
- if(removeBoundary(*cit)) { smoothened = true; }
+ smoothened = true;
+ round = 0;
+ while(smoothened){
+ round++;
+ smoothened = false;
+ for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
+ if(straightenChain(*cit)){smoothened = true;}
+ if(removeBoundary(*cit)) { smoothened = true; }
+ }
+ if(round > MAXROUNDS) smoothened = false;
}
- if(round > MAXROUNDS) smoothened = false;
- }
- eraseNullCells();
- if(this->getSize() < bestChain->getSize()) bestChain = this;
- printf("%d-chain simulated annealing removed %d cells. \n", getDim(), before - getSize());
+
+
+ eraseNullCells();
+ if(this->getSize() < bestChain->getSize()) bestChain = this;
+ printf("%d-chain simulated annealing removed %d cells. \n", getDim(), before - getSize());
}
+ */
double t2 = Cpu();
printf("Smoothened a %d-chain from %d cells to %d cells (%g s).\n", getDim(), start, getSize(), t2-t1);
return;
diff --git a/Geo/ChainComplex.h b/Geo/ChainComplex.h
index e18aec776e..6a7dc87091 100644
--- a/Geo/ChainComplex.h
+++ b/Geo/ChainComplex.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.
+// Contributed by Matti Pellikka <matti.pellikka@tut.fi>.
#ifndef _CHAINCOMPLEX_H_
#define _CHAINCOMPLEX_H_
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index f2ca589ee6..198f0566d2 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.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.
+// Contributed by Matti Pellikka <matti.pellikka@tut.fi>.
#include "Homology.h"
diff --git a/Geo/Homology.h b/Geo/Homology.h
index 5308ce105a..2c8e69a5bc 100644
--- a/Geo/Homology.h
+++ b/Geo/Homology.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.
+// Contributed by Matti Pellikka <matti.pellikka@tut.fi>.
#ifndef _HOMOLOGY_H_
#define _HOMOLOGY_H_
@@ -36,15 +36,19 @@ class Homology
Homology(GModel* model, std::vector<int> physicalDomain, std::vector<int> physicalSubdomain);
~Homology(){ delete _cellComplex; }
+ // Find the generators/duals of homology spaces, or just compute the ranks of homology spaces
void findGenerators(std::string fileName);
void findDualGenerators(std::string fileName);
-
void computeBettiNumbers();
+
void swapSubdomain() { _cellComplex->swapSubdomain(); }
-
+
+ // Restore the cell complex to its original state before cell reductions
+ void restoreHomology() { _cellComplex->restoreComplex(); }
+
+ // Create a string describing the generator
std::string getDomainString() {
-
std::string domainString = "({";
for(unsigned int i = 0; i < _domain.size(); i++){
std::string temp = "";
diff --git a/utils/api_demos/mainHomology.cpp b/utils/api_demos/mainHomology.cpp
index 4fe82d739b..acf2812edc 100644
--- a/utils/api_demos/mainHomology.cpp
+++ b/utils/api_demos/mainHomology.cpp
@@ -1,14 +1,10 @@
-// configure, compile and install Gmsh as a library with
-//
-// ./configure --disable-gui --disable-netgen --disable-chaco --prefix=/usr/local/
-// make install-lib
-//
-// Then compile this driver with "g++ celldriver.cpp -lGmsh -llapack -lblas -lgmp"
-//
-//
-// This program creates .msh file chains.msh which represents the generators and
-// the cuts of an two port transformer model's homology groups.
-
+// 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>.
+//
#include <stdio.h>
#include <sstream>
@@ -17,130 +13,36 @@
#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();
-}
-
+#include <gmsh/Homology.h>
int main(int argc, char **argv)
{
GmshInitialize(argc, argv);
GModel *m = new GModel();
- m->readGEO("transformer.geo");
- m->mesh(3);
-
- printf("This model has: %d GRegions, %d GFaces, %d GEdges and %d GVertices. \n" , m->getNumRegions(), m->getNumFaces(), m->getNumEdges(), m->getNumVertices());
-
- std::vector<GEntity*> domain;
- std::vector<GEntity*> subdomain;
-
- // whole model the domain
- for(GModel::riter rit = m->firstRegion(); rit != m->lastRegion(); rit++){
- GEntity *r = *rit;
- domain.push_back(r);
- }
-
- // the ports as subdomain
- GModel::fiter sub = m->firstFace();
- GEntity *s = *sub;
- subdomain.push_back(s);
- s = *(++sub);
- subdomain.push_back(s);
- s =*(++sub);
- subdomain.push_back(s);
- s= *(++sub);
- subdomain.push_back(s);
-
- // create a cell complex
- CellComplex* complex = new CellComplex(domain, subdomain);
-
- // write the cell complex to a .msh file
- complex->writeComplexMSH("chains.msh");
-
- 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));
- // reduce the complex in order to ease homology computation
- complex->reduceComplex(true);
-
- // perform series of cell combinatios in order to further ease homology computation
- complex->combine(3);
- complex->combine(2);
- complex->combine(1);
-
- // 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();
-
- // 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;
-
- }
- }
-
- delete chains;
- delete complex;
-
- // create new complex to compute the cuts
- CellComplex* complex2 = new CellComplex(domain, subdomain);
-
- // reduce the complex by coreduction, this removes all vertices, hence 0 as an argument
- complex2->coreduceComplex(true);
-
- // perform series of "dual complex" cell combinations in order to ease homology computation
- complex2->cocombine(0);
- complex2->cocombine(1);
- complex2->cocombine(2);
+ m->readGEO("model.geo");
+ m->mesh(3);
+ // OR
+ // m->readMSH("model.msh");
- // create a chain complex of a cell complex (construct boundary operator matrices)
- ChainComplex* dualChains = new ChainComplex(complex2);
+ // List of physical regions as domain for homology computation (relative to subdomain).
+ std::vector<int> domain;
+ std::vector<int> subdomain;
- // transpose the boundary operator matrices,
- // these are the boundary operator matrices for the dual chain complex
- dualChains->transposeHMatrices();
+ Homology* homology = new Homology(m, domain, subdomain);
+ std::string fileName = "homology_generators.msh";
+ homology->findGenerators(fileName);
+ homology->restoreHomology();
- // compute the homology of the dual complex
- dualChains->computeHomology(true);
+ Homology* homology = new Homology(m, domain, subdomain);
+ fileName = "homology_dualgenerators.msh";
+ homology->findDualGenerators(fileName);
+ homology->restoreHomology();
- // Append the duals of the cuts of the homology generators to the .msh file.
- for(int j = 3; j > -1; j--){
- for(int i = 1; i <= dualChains->getBasisSize(j); i++){
-
- std::string generator;
- std::string dimension;
- convert(i, generator);
- convert(3-j, dimension);
-
- std::string name = dimension + "D Dual cut " + generator;
- Chain* chain = new Chain(complex2->getCells(j), dualChains->getCoeffVector(j,i), complex2, name);
- chain->writeChainMSH("chains.msh");
- delete chain;
- }
- }
-
- delete dualChains;
-
+ Homology* homology = new Homology(m, domain, subdomain);
+ homology->computeBettiNumbers();
+ delete homology;
- delete complex2;
delete m;
GmshFinalize();
--
GitLab