diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index 2d8fb36c989e8583c1c0a0452614ee542867e9bc..fabbf3b36f9a901345eb639c99ebdcf585c5b89f 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -668,7 +668,7 @@ void ChainComplex::smoothenChain(std::map<Cell*, int, Less_Cell>& cells)
}
eraseNullCells(cells);
int size = cells.size();
- Msg::Info("Smoothened a %d-chain from %d cells to %d cells",
+ Msg::Info("Simplified a %d-chain from %d cells to %d cells",
dim, start, size);
}
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index ed620e77db2729e6701566bcd4bfd8ebb2f8ae54..28231946bfffca78d85056541b91d12982d953a2 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -3220,19 +3220,25 @@ void GModel::computeHomology()
}
std::string dims = ss.str();
- if(type == "Homology") {
+ if(type == "Homology" && !homology->isHomologyComputed(dim)) {
+
homology->findHomologyBasis(dim);
+
Msg::Info("Homology space basis chains to save: %s.", dims.c_str());
for(unsigned int i = 0; i < dim.size(); i++)
if(dim.at(i) >= 0 && dim.at(i) <= getDim())
homology->addChainsToModel(dim.at(i));
+
}
- else if(type == "Cohomology") {
+ else if(type == "Cohomology" && !homology->isCohomologyComputed(dim)) {
+
homology->findCohomologyBasis(dim);
+
Msg::Info("Cohomology space basis cochains to save: %s.", dims.c_str());
for(unsigned int i = 0; i < dim.size(); i++)
if(dim.at(i) >= 0 && dim.at(i) <= getDim())
homology->addCochainsToModel(dim.at(i));
+
}
else
Msg::Error("Unknown type of homology computation: %s", type.c_str());
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index a9db40d532232c597332285637a8e4d78e1435c9..69a14f6d1aa0ea7caae57fe3e07b46466df537f9 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -137,8 +137,8 @@ void Homology::_deleteChains(std::vector<int> dim)
if(d < 0 || d > 3) continue;
for(unsigned int i = 0; i < _chains[d].size(); i++) {
delete _chains[d].at(i);
- _chains[d].clear();
}
+ _chains[d].clear();
_homologyComputed[d] = false;
}
}
@@ -157,8 +157,8 @@ void Homology::_deleteCochains(std::vector<int> dim)
if(d < 0 || d > 3) continue;
for(unsigned int i = 0; i < _cochains[d].size(); i++) {
delete _cochains[d].at(i);
- _cochains[d].clear();
}
+ _cochains[d].clear();
_cohomologyComputed[d] = false;
}
}
@@ -179,7 +179,7 @@ void Homology::findHomologyBasis(std::vector<int> dim)
double t1 = Cpu();
int omitted = _cellComplex->reduceComplex(_combine, _omit);
- if(!dim.empty() && _combine > 1) {
+ if(!dim.empty() && _combine > 1 && !_smoothen) {
for(int i = 1; i <= 3; i++) {
if(!std::binary_search(dim.begin(), dim.end(), i)) {
_cellComplex->cocombine(i-1);
@@ -297,7 +297,7 @@ void Homology::findCohomologyBasis(std::vector<int> dim)
std::string name = "H^" + dimension + domain + generator;
std::map<Cell*, int, Less_Cell> chain;
- chainComplex.getBasisChain(chain, i, j, 3, _smoothen);
+ chainComplex.getBasisChain(chain, i, j, 3, false);
int torsion = chainComplex.getTorsion(j,i);
if(!chain.empty()) {
_createChain(chain, name, true);
@@ -331,6 +331,102 @@ void Homology::findCohomologyBasis(std::vector<int> dim)
}
}
+bool Homology::isHomologyComputed(std::vector<int> dim)
+{
+ if(dim.empty()) return (_homologyComputed[0] &&
+ _homologyComputed[1] &&
+ _homologyComputed[2] &&
+ _homologyComputed[3]);
+
+ bool computed = true;
+ for(unsigned int i = 0; i < dim.size(); i++) {
+ int d = dim.at(i);
+ if(d < 0 || d > 3) continue;
+ computed = computed && _homologyComputed[d];
+ }
+ return computed;
+}
+
+bool Homology::isCohomologyComputed(std::vector<int> dim)
+{
+ if(dim.empty()) return (_cohomologyComputed[0] &&
+ _cohomologyComputed[1] &&
+ _cohomologyComputed[2] &&
+ _cohomologyComputed[3]);
+
+ bool computed = true;
+ for(unsigned int i = 0; i < dim.size(); i++) {
+ int d = dim.at(i);
+ if(d < 0 || d > 3) continue;
+ computed = computed && _cohomologyComputed[d];
+ }
+ return computed;
+}
+
+void Homology::findCompatibleBasisPair(int master, std::vector<int> dim)
+{
+ if(!this->isHomologyComputed(dim)) this->findHomologyBasis(dim);
+ if(!this->isCohomologyComputed(dim)) this->findCohomologyBasis(dim);
+ for(unsigned int idim = 0 ; idim < dim.size(); idim++) {
+ int d = dim.at(idim);
+ if(d < 1 || d > 2) continue;
+ int n = this->betti(d);
+ if(n < 2) continue;
+ if((int)_chains[d].size() != n || (int)_cochains[d].size() != n) {
+ Msg::Warning("Cannot produce compatible %d-(co)homology bases.", d);
+ Msg::Debug("%d basis %d-chains and %d basis %d-cochains.",
+ (int)_chains[d].size(), d, (int)_cochains[d].size(), d);
+ continue;
+ }
+ fullMatrix<double> m(n,n);
+ for(int i = 0; i < n; i++) {
+ for(int j = 0; j < n; j++) {
+ if(master==0) m(i,j) = incidence(*_cochains[d].at(i), *_chains[d].at(j));
+ else m(i,j) = incidence(*_chains[d].at(i), *_cochains[d].at(j));
+ }
+ }
+
+ int det = m.determinant();
+ if(abs(det) != 1 || !m.invertInPlace()) {
+ Msg::Warning("Cannot produce compatible %d-(co)homology bases.", d);
+ Msg::Debug("Incidence matrix: ");
+ for(int i = 0; i < n; i++)
+ for(int j = 0; j < n; j++)
+ Msg::Debug("(%d, %d) = %d", i, j, m(i,j));
+ continue;
+ }
+
+ std::vector<Chain<int>*> newBasis(n, NULL);
+
+ if(master==0) {
+ for(int i = 0; i < n; i++) {
+ newBasis.at(i) = new Chain<int>();
+ for(int j = 0; j < n; j++) {
+ *newBasis.at(i) += (int)m(i,j)*(*_cochains[d].at(j));
+ }
+ }
+ for(int i = 0; i < n; i++) {
+ newBasis.at(i)->setName(_cochains[d].at(i)->getName());
+ delete _cochains[d].at(i);
+ _cochains[d].at(i) = newBasis.at(i);
+ }
+ }
+ else {
+ for(int i = 0; i < n; i++) {
+ newBasis.at(i) = new Chain<int>();
+ for(int j = 0; j < n; j++) {
+ *newBasis.at(i) += (int)m(i,j)*(*_chains[d].at(j));
+ }
+ }
+ for(int i = 0; i < n; i++) {
+ newBasis.at(i)->setName(_chains[d].at(i)->getName());
+ delete _chains[d].at(i);
+ _chains[d].at(i) = newBasis.at(i);
+ }
+ }
+ }
+}
+
void Homology::addChainsToModel(int dim, bool post, int physicalNumRequest)
{
int pgnum = -1;
diff --git a/Geo/Homology.h b/Geo/Homology.h
index 0af68cbd80ef5d7619a7f91230f6515d89a2dcee..386e6a7e7f122d448a4d6588b1b462079de87a62 100644
--- a/Geo/Homology.h
+++ b/Geo/Homology.h
@@ -99,12 +99,28 @@ class Homology
GModel* getModel() const { return _model; }
void setFileName(std::string fileName) { _fileName = fileName; }
+ void getDomain(std::vector<int>& domain) { domain = _domain; }
+ void getSubdomain(std::vector<int>& subdomain) { subdomain = _subdomain; }
+
// find the bases of (co)homology spaces
// if dim is empty, find 0-,1-,2-,3-(co)homology spaces bases
// otherwise only find those indicated in dim
void findHomologyBasis(std::vector<int> dim=std::vector<int>());
void findCohomologyBasis(std::vector<int> dim=std::vector<int>());
+ // find a homology and cohomology basis pair such that
+ // the incidence matrix of the bases is an identity matrix
+ // if master==0, homology basis determines the cohomology basis
+ // if dim is empty, find 0-,1-,2-,3-(co)homology spaces bases
+ // otherwise only find those indicated in dim
+ void findCompatibleBasisPair(int master=0,
+ std::vector<int> dim=std::vector<int>());
+
+ // is the (co)homology in given dimensions already compited
+ // if dim is empty, return true only if computed in all dimensions
+ bool isHomologyComputed(std::vector<int> dim=std::vector<int>());
+ bool isCohomologyComputed(std::vector<int> dim=std::vector<int>());
+
// add chains to Gmsh model
// dim: only add dim-chains if dim != -1
// post: create a post-processing view
diff --git a/Plugin/HomologyComputation.cpp b/Plugin/HomologyComputation.cpp
index 9d52c8714a0800318f9ed08980208849a713dd47..e3794c78d4521336fdb3f6daf6368d0b306d2e27 100644
--- a/Plugin/HomologyComputation.cpp
+++ b/Plugin/HomologyComputation.cpp
@@ -25,7 +25,7 @@ StringXNumber HomologyComputationOptions_Number[] = {
{GMSH_FULLRC, "CreatePostProcessingViews", NULL, 1.},
{GMSH_FULLRC, "ReductionOmit", NULL, 1.},
{GMSH_FULLRC, "ReductionCombine", NULL, 2.},
- {GMSH_FULLRC, "PostProcessSmoothen", NULL, 1.}
+ {GMSH_FULLRC, "PostProcessSimplify", NULL, 1.}
};
StringXString HomologyComputationOptions_String[] = {