diff --git a/Plugin/CMakeLists.txt b/Plugin/CMakeLists.txt
index 671de2ce94ad8e26857449310a39b8be33241567..a557fd8eca265d6170fb72010d6298b5778a0a36 100644
--- a/Plugin/CMakeLists.txt
+++ b/Plugin/CMakeLists.txt
@@ -22,7 +22,7 @@ set(SRC
Annotate.cpp Remove.cpp
Probe.cpp MinMax.cpp
HarmonicToTime.cpp ModulusPhase.cpp
- HomologyComputation.cpp
+ HomologyComputation.cpp HomologyPostProcessing.cpp
Distance.cpp ExtractEdges.cpp NearestNeighbor.cpp
AnalyseCurvedMesh.cpp
FieldFromAmplitudePhase.cpp
diff --git a/Plugin/HomologyComputation.cpp b/Plugin/HomologyComputation.cpp
index 343337e6359b5f2f6f93cb431a20b376b97de50c..5c38a6e77f23a472709de87185de8d4cd9823412 100644
--- a/Plugin/HomologyComputation.cpp
+++ b/Plugin/HomologyComputation.cpp
@@ -6,6 +6,9 @@
// Contributed by Matti Pellikka <matti.pellikka@tut.fi>.
#include <stdlib.h>
+#include <string>
+#include <iostream>
+#include <sstream>
#include "Gmsh.h"
#include "GmshConfig.h"
#include "GModel.h"
@@ -15,15 +18,14 @@
#if defined(HAVE_KBIPACK)
StringXNumber HomologyComputationOptions_Number[] = {
- {GMSH_FULLRC, "PhysicalGroupForDomain1", NULL, 0.},
- {GMSH_FULLRC, "PhysicalGroupForDomain2", NULL, 0.},
- {GMSH_FULLRC, "PhysicalGroupForSubdomain1", NULL, 0.},
- {GMSH_FULLRC, "PhysicalGroupForSubdomain2", NULL, 0.},
- {GMSH_FULLRC, "CompututeHomology", NULL, 1.},
- {GMSH_FULLRC, "ComputeCohomology", NULL, 0.},
+ {GMSH_FULLRC, "ComputeHomology", NULL, 1.},
+ {GMSH_FULLRC, "ComputeCohomology", NULL, 0.}
};
StringXString HomologyComputationOptions_String[] = {
+ {GMSH_FULLRC, "DomainPhysicalGroups", NULL, ""},
+ {GMSH_FULLRC, "SubdomainPhysicalGroups", NULL, ""},
+ {GMSH_FULLRC, "DimensionOfChainsToSave", NULL, "1, 2"},
{GMSH_FULLRC, "Filename", NULL, "homology.msh"}
};
@@ -37,16 +39,16 @@ extern "C"
std::string GMSH_HomologyComputationPlugin::getHelp() const
{
- return "Plugin(Homology) computes ranks and basis elements "
- "of (relative) homology and cohomology spaces.\n\n"
+ return "Plugin(HomologyComputation) computes representative chains "
+ "of basis elements of (relative) homology and cohomology spaces.\n\n"
"Define physical groups in order to specify the computation "
"domain and the relative subdomain. Otherwise the whole mesh "
"is the domain and the relative subdomain is empty. \n\n"
- "Plugin(Homology) creates new views, one for each basis element. "
- "The resulting basis chains together with the mesh are saved to "
- "the file given.";
+ "Plugin(HomologyComputation) creates new views, one for each "
+ "basis element. The resulting basis chains of desired dimension "
+ "together with the mesh are saved to the given file.";
}
int GMSH_HomologyComputationPlugin::getNbOptions() const
@@ -69,27 +71,42 @@ StringXString *GMSH_HomologyComputationPlugin::getOptionStr(int iopt)
return &HomologyComputationOptions_String[iopt];
}
+bool GMSH_HomologyComputationPlugin::parseStringOpt
+(int stringOpt, std::vector<int>& intList)
+{
+ std::string list = HomologyComputationOptions_String[stringOpt].def;
+ intList.clear();
+
+ int n;
+ char a;
+ std::istringstream ss(list);
+ while(ss >> n){
+ intList.push_back(n);
+ if(ss >> a){
+ if(a != ',') {
+ Msg::Error("Unexpected character \'%c\' while parsing \'%s\'", a,
+ HomologyComputationOptions_String[stringOpt].str);
+ return false;
+ }
+ }
+ }
+ return true;
+}
+
PView *GMSH_HomologyComputationPlugin::execute(PView *v)
{
- std::string fileName = HomologyComputationOptions_String[0].def;
+ std::string fileName = HomologyComputationOptions_String[3].def;
+ int hom = (int)HomologyComputationOptions_Number[0].def;
+ int coh = (int)HomologyComputationOptions_Number[1].def;
std::vector<int> domain;
std::vector<int> subdomain;
+ std::vector<int> dimsave;
+ if(!parseStringOpt(0, domain)) return 0;
+ if(!parseStringOpt(1, subdomain)) return 0;
+ if(!parseStringOpt(2, dimsave)) return 0;
- int d1 = (int)HomologyComputationOptions_Number[0].def;
- int d2 = (int)HomologyComputationOptions_Number[1].def;
- if(d1 > 0) domain.push_back(d1);
- if(d2 > 0) domain.push_back(d2);
- d1 = (int)HomologyComputationOptions_Number[2].def;
- d2 = (int)HomologyComputationOptions_Number[3].def;
- if(d1 > 0) subdomain.push_back(d1);
- if(d2 > 0) subdomain.push_back(d2);
-
-
- int hom = (int)HomologyComputationOptions_Number[4].def;
- int coh = (int)HomologyComputationOptions_Number[5].def;
-
- GModel *m = GModel::current();
+ GModel* m = GModel::current();
Homology* homology = new Homology(m, domain, subdomain);
homology->setFileName(fileName);
@@ -101,6 +118,10 @@ PView *GMSH_HomologyComputationPlugin::execute(PView *v)
cc->restoreComplex();
homology->findCohomologyBasis(cc);
}
+ for(unsigned int i = 0; i < dimsave.size(); i++) {
+ int dim = dimsave.at(i);
+ if(dim > -1 && dim < 4) homology->addChainsToModel(dim);
+ }
delete cc;
delete homology;
diff --git a/Plugin/HomologyComputation.h b/Plugin/HomologyComputation.h
index 2cd18476f414eec2782f171d60fd7612e792ee92..ce9668ea77c6e6a0a998ba2bd52425dd59ba50b4 100644
--- a/Plugin/HomologyComputation.h
+++ b/Plugin/HomologyComputation.h
@@ -22,7 +22,7 @@ class GMSH_HomologyComputationPlugin : public GMSH_PostPlugin
{
public:
GMSH_HomologyComputationPlugin(){}
- std::string getName() const { return "Homology"; }
+ std::string getName() const { return "HomologyComputation"; }
std::string getShortHelp() const
{
return "Compute relative (co)homology spaces";
@@ -34,6 +34,7 @@ class GMSH_HomologyComputationPlugin : public GMSH_PostPlugin
int getNbOptionsStr() const;
StringXString *getOptionStr(int iopt);
PView *execute(PView *);
+ bool parseStringOpt(int stringOpt, std::vector<int>& intList);
};
#endif
diff --git a/Plugin/HomologyPostProcessing.cpp b/Plugin/HomologyPostProcessing.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3511f97881a41c948ce30ae744cd0eb082f28fc4
--- /dev/null
+++ b/Plugin/HomologyPostProcessing.cpp
@@ -0,0 +1,361 @@
+// Gmsh - Copyright (C) 1997-2012 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 <stdlib.h>
+#include <string>
+#include <iostream>
+#include <sstream>
+#include "Gmsh.h"
+#include "GmshConfig.h"
+#include "GModel.h"
+#include "Chain.h"
+#include "HomologyPostProcessing.h"
+
+#if defined(HAVE_KBIPACK)
+
+StringXNumber HomologyPostProcessingOptions_Number[] = {
+ {GMSH_FULLRC, "DimensionOfChains", NULL, 1},
+ {GMSH_FULLRC, "ApplyBoundaryOperatorToResults", NULL, 0}
+};
+
+StringXString HomologyPostProcessingOptions_String[] = {
+ {GMSH_FULLRC, "TransformationMatrix", NULL, "1, 0; 0, 1"},
+ {GMSH_FULLRC, "PhysicalGroupsOfOperatedChains", NULL, "1, 2"},
+ {GMSH_FULLRC, "PhysicalGroupsOfOperatedChains2", NULL, ""},
+ {GMSH_FULLRC, "PhysicalGroupsToTraceResults", NULL, ""},
+ {GMSH_FULLRC, "NameForResultChains", NULL, "c"},
+};
+
+extern "C"
+{
+ GMSH_Plugin *GMSH_RegisterHomologyPostProcessingPlugin()
+ {
+ return new GMSH_HomologyPostProcessingPlugin();
+ }
+}
+
+std::string GMSH_HomologyPostProcessingPlugin::getHelp() const
+{
+ return "Plugin(HomologyPostProcessing) operates on representative "
+ "basis chains of homology and cohomology spaces.\n\n"
+
+ "Define the physical groups of chains and their dimension "
+ "to be transformed by the "
+ "given integer transformation matrix. Number of matrix columns must "
+ "equal to the number of chains given. The resulting "
+ "chains are thus integer combinations of the given chains. ";
+}
+
+int GMSH_HomologyPostProcessingPlugin::getNbOptions() const
+{
+ return sizeof(HomologyPostProcessingOptions_Number) / sizeof(StringXNumber);
+}
+
+StringXNumber *GMSH_HomologyPostProcessingPlugin::getOption(int iopt)
+{
+ return &HomologyPostProcessingOptions_Number[iopt];
+}
+
+int GMSH_HomologyPostProcessingPlugin::getNbOptionsStr() const
+{
+ return sizeof(HomologyPostProcessingOptions_String) / sizeof(StringXString);
+}
+
+StringXString *GMSH_HomologyPostProcessingPlugin::getOptionStr(int iopt)
+{
+ return &HomologyPostProcessingOptions_String[iopt];
+}
+
+bool GMSH_HomologyPostProcessingPlugin::parseStringOpt
+(int stringOpt, std::vector<int>& intList)
+{
+ std::string list = HomologyPostProcessingOptions_String[stringOpt].def;
+ intList.clear();
+
+ int n;
+ char a;
+ std::istringstream ss(list);
+ while(ss >> n){
+ intList.push_back(n);
+ if(ss >> a){
+ if(a != ',') {
+ Msg::Error("Unexpected character \'%c\' while parsing \'%s\'", a,
+ HomologyPostProcessingOptions_String[stringOpt].str);
+ return false;
+ }
+ }
+ }
+ return true;
+}
+
+void GMSH_HomologyPostProcessingPlugin::createChains
+(int dim, GModel* m,
+ const std::vector<GEntity*>& chainEntities,
+ const std::vector<std::string>& chainNames,
+ std::vector<Chain<int> >& chains)
+{
+ for(unsigned int i = 0; i < chainEntities.size(); i++) {
+ GEntity* e = chainEntities.at(i);
+ Chain<int> chain;
+ for(unsigned int j = 0; j < e->getNumMeshElements(); j++) {
+ chain.addMeshElement(e->getMeshElement(j));
+ }
+ chains.push_back(chain);
+ chains.at(i).setName(chainNames.at(i));
+ }
+}
+
+int GMSH_HomologyPostProcessingPlugin::detIntegerMatrix
+(std::vector<int>& matrix)
+{
+ int n = sqrt(matrix.size());
+ fullMatrix<double> m(n,n);
+ for(int i = 0; i < n; i++)
+ for(int j = 0; j < n; j++)
+ m(i,j) = matrix.at(i*n+j);
+
+ return m.determinant();
+}
+
+bool GMSH_HomologyPostProcessingPlugin::invertIntegerMatrix
+(std::vector<int>& matrix)
+{
+ int n = sqrt(matrix.size());
+ fullMatrix<double> m(n,n);
+ for(int i = 0; i < n; i++)
+ for(int j = 0; j < n; j++)
+ m(i,j) = matrix.at(i*n+j);
+
+ if(!m.invertInPlace())
+ Msg::Error("Matrix is not unimodular");
+
+ for(int i = 0; i < n; i++)
+ for(int j = 0; j < n; j++)
+ matrix.at(i*n+j) = m(i,j);
+}
+
+PView *GMSH_HomologyPostProcessingPlugin::execute(PView *v)
+{
+ std::string matrixString = HomologyPostProcessingOptions_String[0].def;
+ std::string opString1 = HomologyPostProcessingOptions_String[1].def;
+ std::string opString2 = HomologyPostProcessingOptions_String[2].def;
+ std::string cname = HomologyPostProcessingOptions_String[4].def;
+ std::string traceString = HomologyPostProcessingOptions_String[3].def;
+ int dim = (int)HomologyPostProcessingOptions_Number[0].def;
+ int bd = (int)HomologyPostProcessingOptions_Number[1].def;
+
+ if(dim < 0 || dim > 3) {
+ Msg::Error("Invalid cell dimension %d", dim);
+ return 0;
+ }
+
+ GModel* m = GModel::current();
+
+ int n;
+ char a;
+ int cols = 0;
+ int col = 0;
+ std::vector<int> matrix;
+ if(matrixString != "I") {
+ std::istringstream ss(matrixString);
+ while(ss >> n){
+ matrix.push_back(n);
+ col++;
+ if(ss >> a){
+ if(a != ',' && a != ';') {
+ Msg::Error("Unexpected character \'%c\' while parsing \'%s\'", a,
+ HomologyPostProcessingOptions_String[0].str);
+ return 0;
+ }
+ if(a == ';') {
+ if(cols != 0 && cols != col) {
+ Msg::Error("Number of columns must match (%d != %d)", cols, col);
+ return 0;
+ }
+ cols = col;
+ col = 0;
+ }
+ }
+ }
+ if(cols != 0 && cols != col && col != 0) {
+ Msg::Error("Number of columns must match (%d != %d)", cols, col);
+ return 0;
+ }
+ if(cols == 0) cols = col;
+ }
+
+ int rows = 0;
+ if(!matrix.empty()) {
+ rows = matrix.size()/cols;
+ if((int)matrix.size() % cols != 0) {
+ Msg::Error("Number of matrix rows and columns aren't compatible (residual: %d)",
+ (int)matrix.size() % cols);
+ return 0;
+ }
+ }
+
+ std::vector<int> basisPhysicals;
+ if(!parseStringOpt(1, basisPhysicals)) return 0;
+ std::vector<int> basisPhysicals2;
+ if(!parseStringOpt(2, basisPhysicals2)) return 0;
+
+ if(matrixString != "I" &&
+ basisPhysicals.size() != cols &&
+ basisPhysicals2.empty()) {
+ Msg::Error("Number of matrix columns and operated chains must match (%d != %d)", cols, basisPhysicals.size());
+ return 0;
+ }
+ else if(matrixString == "I") {
+ cols = basisPhysicals.size();
+ rows = cols;
+ matrix = std::vector<int>(cols*cols, 0);
+ for(int i = 0; i < cols; i++)
+ matrix.at(i*cols+i) = 1;
+ }
+
+ if(!basisPhysicals2.empty() &&
+ basisPhysicals.size() != basisPhysicals2.size()) {
+ Msg::Error("Number of operated chains must match (%d != %d)",
+ basisPhysicals.size(), basisPhysicals2.size());
+ return 0;
+ }
+
+ std::vector<int> tracePhysicals;
+ if(!parseStringOpt(3, tracePhysicals)) return 0;
+
+ std::vector<GEntity*> basisEntities;
+ std::vector<GEntity*> basisEntities2;
+ std::map<int, std::vector<GEntity*> > groups[4];
+ m->getPhysicalGroups(groups);
+ std::map<int, std::vector<GEntity*> >::iterator it;
+
+ std::vector<std::string> curBasisNames;
+ for(unsigned int i = 0; i < basisPhysicals.size(); i++){
+ curBasisNames.push_back(m->getPhysicalName(dim, basisPhysicals.at(i)));
+ it = groups[dim].find(basisPhysicals.at(i));
+ if(it != groups[dim].end()){
+ std::vector<GEntity*> physicalGroup = it->second;
+ for(unsigned int k = 0; k < physicalGroup.size(); k++){
+ basisEntities.push_back(physicalGroup.at(k));
+ break;
+ }
+ }
+ else {
+ Msg::Error("%d-dimensional physical group %d does not exist",
+ dim, basisPhysicals.at(i));
+ return 0;
+ }
+ }
+ std::vector<std::string> curBasisNames2;
+ for(unsigned int i = 0; i < basisPhysicals2.size(); i++){
+ curBasisNames2.push_back(m->getPhysicalName(dim, basisPhysicals2.at(i)));
+ it = groups[dim].find(basisPhysicals2.at(i));
+ if(it != groups[dim].end()){
+ std::vector<GEntity*> physicalGroup = it->second;
+ for(unsigned int k = 0; k < physicalGroup.size(); k++){
+ basisEntities2.push_back(physicalGroup.at(k));
+ break;
+ }
+ }
+ else {
+ Msg::Error("%d-dimensional physical group %d does not exist",
+ dim, basisPhysicals2.at(i));
+ return 0;
+ }
+ }
+
+ std::vector<Chain<int> > curBasis;
+ createChains(dim, m, basisEntities, curBasisNames, curBasis);
+ std::vector<Chain<int> > curBasis2;
+ createChains(dim, m, basisEntities2, curBasisNames2, curBasis2);
+
+
+ if(!curBasis2.empty()) {
+ rows = curBasis2.size();
+ cols = curBasis.size();
+ matrix = std::vector<int>(rows*cols, 0);
+ for(int i = 0; i < rows; i++)
+ for(int j = 0; j < cols; j++)
+ matrix.at(i*cols+j) = incidence(curBasis2.at(i), curBasis.at(j));
+ }
+
+ if(!curBasis2.empty())
+ Msg::Debug("Incidence matrix: ");
+ else
+ Msg::Debug("Transformation matrix: ");
+ for(int i = 0; i < rows; i++)
+ for(int j = 0; j < cols; j++)
+ Msg::Debug("(%d, %d) = %d", i, j, matrix.at(i*cols+j));
+
+ std::vector<Chain<int> > newBasis(rows, Chain<int>());
+ if(!curBasis2.empty()) {
+ Msg::Info("Computing new basis %d-chains such that the incidence matrix of %d-chain bases %s and %s is the indentity matrix",
+ dim, dim, opString1.c_str(), opString2.c_str());
+ if(!invertIntegerMatrix(matrix)) return 0;
+ for(int i = 0; i < rows; i++)
+ for(int j = 0; j < cols; j++)
+ newBasis.at(i) += matrix.at(i*cols+j)*curBasis2.at(j); }
+ else {
+ Msg::Info("Applying transformation matrix [%s] to %d-chains %s",
+ matrixString.c_str(), dim, opString1.c_str());
+ if(rows == cols) {
+ int det = detIntegerMatrix(matrix);
+ if(det != 1 && det != -1)
+ Msg::Warning("Transformation matrix is not unimodular (det = %d)", det);
+ }
+ for(int i = 0; i < rows; i++)
+ for(int j = 0; j < cols; j++)
+ newBasis.at(i) += matrix.at(i*cols+j)*curBasis.at(j);
+
+ }
+
+ if(bd) {
+ Msg::Info("Applying boundary operator to the result %d-chains", dim);
+ for(unsigned int i = 0; i < newBasis.size(); i++)
+ newBasis.at(i) = newBasis.at(i).getBoundary();
+ }
+
+ if(!tracePhysicals.empty()) {
+ Msg::Info("Taking trace of result %d-chains to domain %s",
+ dim, traceString.c_str());
+ std::vector<GEntity*> traceEntities;
+ for(unsigned int i = 0; i < tracePhysicals.size(); i++){
+ bool found = false;
+ for(int j = 0; j < 4; j++){
+ it = groups[j].find(tracePhysicals.at(i));
+ if(it != groups[j].end()){
+ found = true;
+ std::vector<GEntity*> physicalGroup = it->second;
+ for(unsigned int k = 0; k < physicalGroup.size(); k++){
+ traceEntities.push_back(physicalGroup.at(k));
+ }
+ }
+ }
+ if(!found) {
+ Msg::Error("Physical group %d does not exist",
+ tracePhysicals.at(i));
+ return 0;
+ }
+ }
+ for(unsigned int i = 0; i < newBasis.size(); i++)
+ newBasis.at(i) = newBasis.at(i).getTrace(traceEntities);
+ ElemChain::clearVertexCache();
+ }
+
+ std::string dims = "";
+ std::string nums = "";
+ for(unsigned int i = 0; i < newBasis.size(); i++) {
+ convert(newBasis.at(i).getDim(), dims);
+ convert(i+1, nums);
+ newBasis.at(i).setName("C" + dims + " " + cname + nums);
+ newBasis.at(i).addToModel(m);
+ }
+
+ return 0;
+}
+
+#endif
diff --git a/Plugin/HomologyPostProcessing.h b/Plugin/HomologyPostProcessing.h
new file mode 100644
index 0000000000000000000000000000000000000000..bcaa5e109ef815b4b8f66f9f9e2e89c476ff61f4
--- /dev/null
+++ b/Plugin/HomologyPostProcessing.h
@@ -0,0 +1,49 @@
+// Gmsh - Copyright (C) 1997-2012 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>.
+
+#ifndef _HOMOLOGY_POST_PROCESSING_H_
+#define _HOMOLOGY_POST_PROCESSING_H_
+
+#include <string>
+#include "Plugin.h"
+#include "Chain.h"
+
+#if defined(HAVE_KBIPACK)
+
+extern "C"
+{
+ GMSH_Plugin *GMSH_RegisterHomologyPostProcessingPlugin();
+}
+
+class GMSH_HomologyPostProcessingPlugin : public GMSH_PostPlugin
+{
+ public:
+ GMSH_HomologyPostProcessingPlugin(){}
+ std::string getName() const { return "HomologyPostProcessing"; }
+ std::string getShortHelp() const
+ {
+ return "Post-process (co)homology space bases";
+ }
+ std::string getHelp() const;
+ std::string getAuthor() const { return "M. Pellikka"; }
+ int getNbOptions() const;
+ StringXNumber *getOption(int iopt);
+ int getNbOptionsStr() const;
+ StringXString *getOptionStr(int iopt);
+ PView *execute(PView *);
+ bool parseStringOpt(int stringOpt, std::vector<int>& intList);
+ void createChains(int dim, GModel* m,
+ const std::vector<GEntity*>& chainEntities,
+ const std::vector<std::string>& chainNames,
+ std::vector<Chain<int> >& chains);
+ bool invertIntegerMatrix(std::vector<int>& matrix);
+ int detIntegerMatrix(std::vector<int>& matrix);
+};
+
+#endif
+
+#endif
diff --git a/Plugin/PluginManager.cpp b/Plugin/PluginManager.cpp
index dbb85ddf591754f31e8596a9632461808c8185bb..7692045db9dffc88029828059e4603e7ae6c83af 100644
--- a/Plugin/PluginManager.cpp
+++ b/Plugin/PluginManager.cpp
@@ -48,6 +48,7 @@
#include "Probe.h"
#include "GSHHS.h"
#include "HomologyComputation.h"
+#include "HomologyPostProcessing.h"
#include "ExtractEdges.h"
#include "FieldFromAmplitudePhase.h"
#include "Bubbles.h"
@@ -245,6 +246,8 @@ void PluginManager::registerDefaultPlugins()
#if defined(HAVE_KBIPACK)
allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
("Homology", GMSH_RegisterHomologyComputationPlugin()));
+ allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
+ ("HomologyPost", GMSH_RegisterHomologyPostProcessingPlugin()));
#endif
#if defined(HAVE_SOLVER)
allPlugins.insert(std::pair<std::string, GMSH_Plugin*>