diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 726882c9f41d9103f9a06a66da04a0ca45dfcec1..b6f025039847ac954258a87a7b8c2644cfa6cb7e 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -7,6 +7,9 @@
#include "CellComplex.h"
#include "MElement.h"
+#include "OS.h"
+
+double CellComplex::_patience = 10;
CellComplex::CellComplex(GModel* model,
std::vector<MElement*>& domainElements,
@@ -72,7 +75,10 @@ bool CellComplex::_insertCells(std::vector<MElement*>& elements,
}
_dim = 0;
+ double t1 = Cpu();
+
for(unsigned int i=0; i < elements.size(); i++){
+
MElement* element = elements.at(i);
int dim = element->getDim();
int type = element->getType();
@@ -112,6 +118,15 @@ bool CellComplex::_insertCells(std::vector<MElement*>& elements,
_biggestCell = biggestElement[_dim];
for (int dim = 3; dim > 0; dim--){
+
+ double t2 = Cpu();
+ if(t2-t1 > CellComplex::_patience && dim > 1) {
+ if(domain == 0)
+ Msg::Info(" ... creating domain %d-cells", dim);
+ else if(domain == 1)
+ Msg::Info(" ... creating subdomain %d-cells", dim);
+ }
+
for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
Cell* cell = *cit;
for(int i = 0; i < cell->getNumBdElements(); i++){
@@ -470,14 +485,13 @@ int CellComplex::reduceComplex(int combine, bool omit, bool homseq)
getSize(3), getSize(2), getSize(1), getSize(0));
if(!getSize(0)) return 0;
+
+ double t1 = Cpu();
int count = 0;
if(relative() && !homseq) removeSubdomain();
std::vector<Cell*> empty;
for(int i = 3; i > 0; i--) count = count + reduction(i, -1, empty);
- Msg::Debug(" %d volumes, %d faces, %d edges, and %d vertices",
- getSize(3), getSize(2), getSize(1), getSize(0));
-
if(omit && !homseq){
std::vector<Cell*> newCells;
@@ -495,15 +509,26 @@ int CellComplex::reduceComplex(int combine, bool omit, bool homseq)
}
}
- Msg::Debug(" %d volumes, %d faces, %d edges, and %d vertices",
- getSize(3), getSize(2), getSize(1), getSize(0));
+ double t2 = Cpu();
+ if(t2-t1 > CellComplex::_patience) {
+ Msg::Info(" .. %d volumes, %d faces, %d edges, and %d vertices",
+ getSize(3), getSize(2), getSize(1), getSize(0));
+ }
if(combine > 0) this->combine(3);
- reduction(2, -1, empty);
+
+ if(combine > 2) for(int i = 3; i > 0; i--) reduction(i, -1, empty);
+ else reduction(2, -1, empty);
+
if(combine > 0) this->combine(2);
- reduction(1, -1, empty);
+
+ if(combine > 2) for(int i = 3; i > 0; i--) reduction(i, -1, empty);
+ else reduction(1, -1, empty);
+
if(combine > 0) this->combine(1);
+ if(combine > 2) for(int i = 3; i > 0; i--) reduction(i, -1, empty);
+
Msg::Debug(" %d volumes, %d faces, %d edges, and %d vertices",
getSize(3), getSize(2), getSize(1), getSize(0));
@@ -542,6 +567,9 @@ int CellComplex::coreduceComplex(int combine, bool omit, int heuristic)
getSize(3), getSize(2), getSize(1), getSize(0));
if(!getSize(0)) return 0;
+
+ double t1 = Cpu();
+
int count = 0;
if(relative()) removeSubdomain();
std::vector<Cell*> empty;
@@ -557,9 +585,6 @@ int CellComplex::coreduceComplex(int combine, bool omit, int heuristic)
for(int j = 1; j <= getDim(); j++)
count += coreduction(j, -1, empty);
- Msg::Debug(" %d volumes, %d faces, %d edges, and %d vertices",
- getSize(3), getSize(2), getSize(1), getSize(0));
-
if(omit){
std::vector<Cell*> newCells;
@@ -589,15 +614,26 @@ int CellComplex::coreduceComplex(int combine, bool omit, int heuristic)
}
- Msg::Debug(" %d volumes, %d faces, %d edges, and %d vertices",
- getSize(3), getSize(2), getSize(1), getSize(0));
+ double t2 = Cpu();
+ if(t2-t1 > CellComplex::_patience) {
+ Msg::Info(" .. %d volumes, %d faces, %d edges, and %d vertices",
+ getSize(3), getSize(2), getSize(1), getSize(0));
+ }
if(combine > 0) this->cocombine(0);
- coreduction(1, -1, empty);
+
+ if(combine > 2) for(int i = 1; i < 4; i++) coreduction(i, -1, empty);
+ else coreduction(1, -1, empty);
+
if(combine > 0) this->cocombine(1);
- coreduction(2, -1, empty);
+
+ if(combine > 2) for(int i = 1; i < 4; i++) coreduction(i, -1, empty);
+ else coreduction(2, -1, empty);
+
if(combine > 0) this->cocombine(2);
- coreduction(3, -1, empty);
+
+ if(combine > 2) for(int i = 1; i < 4; i++) coreduction(i, -1, empty);
+ else coreduction(3, -1, empty);
coherent();
Msg::Debug(" %d volumes, %d faces, %d edges, and %d vertices",
@@ -669,15 +705,26 @@ int CellComplex::combine(int dim)
// getSize(3), getSize(2), getSize(1), getSize(0));
if(dim < 1 || dim > 3) return 0;
+ double t1 = Cpu();
+
std::queue<Cell*> Q;
std::set<Cell*, Less_Cell> Qset;
std::map<Cell*, short int, Less_Cell> bd_c;
int count = 0;
for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
+
+ double t2 = Cpu();
+ if(t2-t1 > CellComplex::_patience) {
+ t1 = Cpu();
+ Msg::Info(" ... %d volumes, %d faces, %d edges, and %d vertices",
+ getSize(3), getSize(2), getSize(1), getSize(0));
+ }
+
Cell* cell = *cit;
cell->getBoundary(bd_c);
enqueueCells(bd_c, Q, Qset);
+
while(Q.size() != 0){
Cell* s = Q.front();
Q.pop();
@@ -741,12 +788,22 @@ int CellComplex::cocombine(int dim)
if(dim < 0 || dim > 2) return 0;
+ double t1 = Cpu();
+
std::queue<Cell*> Q;
std::set<Cell*, Less_Cell> Qset;
std::map<Cell*, short int, Less_Cell> cbd_c;
int count = 0;
for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
+
+ double t2 = Cpu();
+ if(t2-t1 > CellComplex::_patience) {
+ t1 = Cpu();
+ Msg::Info(" ... %d volumes, %d faces, %d edges, and %d vertices",
+ getSize(3), getSize(2), getSize(1), getSize(0));
+ }
+
Cell* cell = *cit;
cell->getCoboundary(cbd_c);
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index b782873549bc2f8ff58969b649789ae052f47f25..e80238ce42844496fd00ad11815e27d18edc8056 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -83,6 +83,8 @@ class CellComplex
int coreduction(Cell* startCell, int omit,
std::vector<Cell*>& omittedCells);
+ static double _patience;
+
public:
CellComplex(GModel* model,
std::vector<MElement*>& domainElements,
diff --git a/Geo/Chain.cpp b/Geo/Chain.cpp
index 29e06b67d0d3d66852c9120aad22280fefa7138f..142315c9bdba782ffda9c1025ca88f5cafc15c85 100644
--- a/Geo/Chain.cpp
+++ b/Geo/Chain.cpp
@@ -48,6 +48,33 @@ inline void ElemChain::_sortVertexIndices()
_si.push_back(it->second);
}
+void findEntitiesInPhysicalGroups
+(GModel* m, const std::vector<int>& physicalGroups,
+std::vector<GEntity*>& entities)
+{
+ std::map<int, std::vector<GEntity*> > groups[4];
+ m->getPhysicalGroups(groups);
+ std::map<int, std::vector<GEntity*> >::iterator it;
+ for(unsigned int i = 0; i < physicalGroups.size(); i++){
+ bool found = false;
+ for(int j = 0; j < 4; j++){
+ it = groups[j].find(physicalGroups.at(i));
+ if(it != groups[j].end()){
+ found = true;
+ std::vector<GEntity*> physicalGroup = it->second;
+ for(unsigned int k = 0; k < physicalGroup.size(); k++){
+ entities.push_back(physicalGroup.at(k));
+ }
+ }
+ }
+ if(!found) {
+ Msg::Error("Physical group %d does not exist",
+ physicalGroups.at(i));
+ }
+ }
+}
+
+
bool ElemChain::_equalVertices(const std::vector<MVertex*>& v2) const {
if(_v.size() != v2.size()) return false;
for(unsigned int i = 0; i < _v.size(); i++)
diff --git a/Geo/Chain.h b/Geo/Chain.h
index 221ee3dabadc2c75bd3b9b06a4bcac9866a7aa93..721a429df04692c1ec4389f2409d83fd871b68e4 100644
--- a/Geo/Chain.h
+++ b/Geo/Chain.h
@@ -164,6 +164,10 @@ class ElemChain : public PosetCat<ElemChain>
};
+void findEntitiesInPhysicalGroups
+(GModel* m, const std::vector<int>& physicalGroups,
+std::vector<GEntity*>& entities);
+
// Class to represent a chain, formal sum of elementary chains
template <class C>
class Chain : public VectorSpaceCat<Chain<C>, C>
@@ -176,6 +180,9 @@ private:
// A name for the chain
std::string _name;
+ Chain<C> _getTraceOrProject(const std::vector<GEntity*>& entities,
+ bool trace) const;
+
public:
// Elementary chain iterators
typedef typename std::map<ElemChain, C>::iterator ecit;
@@ -222,11 +229,22 @@ public:
Chain<C> getBoundary() const;
// Get a chain which contains elementary chains that are
- // fully in the given physical group or elementary entities
+ // in the given physical group or elementary entities
Chain<C> getTrace(GModel* m, int physicalGroup) const;
Chain<C> getTrace(GModel* m, const std::vector<int>& physicalGroups) const;
Chain<C> getTrace(const std::vector<GEntity*>& entities) const;
+ // Get a chain which contains elementary chains that are *not*
+ // in the given physical group or elementary entities
+ Chain<C> getProject(GModel* m, int physicalGroup) const;
+ Chain<C> getProject(GModel* m, const std::vector<int>& physicalGroups) const;
+ Chain<C> getProject(const std::vector<GEntity*>& entities) const;
+
+ // The above two methods decompose a chain c so that
+ // (c - c.getTrace(...) - c.getProject(...)).isZero() == true
+ // holds
+
+
// Add chain to Gmsh model as a physical group,
// elementary chains are turned into mesh elements with
// orientation and multiplicity given by elementary chain coefficient
@@ -239,26 +257,13 @@ public:
template <class C>
Chain<C>::Chain(GModel* m, int physicalGroup)
{
+ std::vector<int> groups(1, physicalGroup);
std::vector<GEntity*> entities;
- std::map<int, std::vector<GEntity*> > groups[4];
- m->getPhysicalGroups(groups);
- std::map<int, std::vector<GEntity*> >::iterator it;
-
- for(int j = 0; j < 4; j++){
- it = groups[j].find(physicalGroup);
- if(it != groups[j].end()){
- _dim = j;
- std::vector<GEntity*> physicalGroup = it->second;
- for(unsigned int k = 0; k < physicalGroup.size(); k++){
- entities.push_back(physicalGroup.at(k));
- }
- }
- }
- if(entities.empty()) {
- Msg::Error("Physical group %d does not exist", physicalGroup);
- }
+ findEntitiesInPhysicalGroups(m, groups, entities);
+
for(unsigned int i = 0; i < entities.size(); i++) {
GEntity* e = entities.at(i);
+ _dim = e->dim();
for(unsigned int j = 0; j < e->getNumMeshElements(); j++) {
this->addMeshElement(e->getMeshElement(j));
}
@@ -356,37 +361,36 @@ Chain<C> Chain<C>::getTrace(GModel* m, int physicalGroup) const
return this->getTrace(m, groups);
}
+template <class C>
+Chain<C> Chain<C>::getProject(GModel* m, int physicalGroup) const
+{
+ std::vector<int> groups(1, physicalGroup);
+ return this->getProject(m, groups);
+}
+
template <class C>
Chain<C> Chain<C>::getTrace(GModel* m,
const std::vector<int>& physicalGroups) const
{
- std::map<int, std::vector<GEntity*> > groups[4];
- m->getPhysicalGroups(groups);
- std::map<int, std::vector<GEntity*> >::iterator it;
std::vector<GEntity*> entities;
- for(unsigned int i = 0; i < physicalGroups.size(); i++){
- bool found = false;
- for(int j = 0; j < 4; j++){
- it = groups[j].find(physicalGroups.at(i));
- if(it != groups[j].end()){
- found = true;
- std::vector<GEntity*> physicalGroup = it->second;
- for(unsigned int k = 0; k < physicalGroup.size(); k++){
- entities.push_back(physicalGroup.at(k));
- }
- }
- }
- if(!found) {
- Msg::Error("Physical group %d does not exist",
- physicalGroups.at(i));
- }
- }
+ findEntitiesInPhysicalGroups(m, physicalGroups, entities);
if(entities.empty()) return Chain<C>();
- return getTrace(entities);
+ return this->_getTraceOrProject(entities, true);
}
template <class C>
-Chain<C> Chain<C>::getTrace(const std::vector<GEntity*>& entities) const
+Chain<C> Chain<C>::getProject(GModel* m,
+ const std::vector<int>& physicalGroups) const
+{
+ std::vector<GEntity*> entities;
+ findEntitiesInPhysicalGroups(m, physicalGroups, entities);
+ if(entities.empty()) return Chain<C>();
+ return this->_getTraceOrProject(entities, false);
+}
+
+template <class C>
+Chain<C> Chain<C>::_getTraceOrProject
+(const std::vector<GEntity*>& entities, bool trace) const
{
Chain<C> result;
for(cecit it = _elemChains.begin(); it != _elemChains.end(); it++) {
@@ -397,11 +401,24 @@ Chain<C> Chain<C>::getTrace(const std::vector<GEntity*>& entities) const
break;
}
}
- if(inDomain) result.addElemChain(it->first, it->second);
+ if(inDomain && trace) result.addElemChain(it->first, it->second);
+ if(!inDomain && !trace) result.addElemChain(it->first, it->second);
}
return result;
}
+template <class C>
+Chain<C> Chain<C>::getTrace(const std::vector<GEntity*>& entities) const
+{
+ return this->_getTraceOrProject(entities, true);
+}
+
+template <class C>
+Chain<C> Chain<C>::getProject(const std::vector<GEntity*>& entities) const
+{
+ return this->_getTraceOrProject(entities, false);
+}
+
template <class C>
void Chain<C>::addMeshElement(MElement* e, C coeff)
{
diff --git a/Geo/Homology.h b/Geo/Homology.h
index 7489cb5743981af90cc05da6f5a1728ace07b0cf..1a997a72c45ad02a21e1739c2166ceb72192b7cb 100644
--- a/Geo/Homology.h
+++ b/Geo/Homology.h
@@ -100,7 +100,7 @@ class Homology
std::vector<int> physicalSubdomain,
std::vector<int> physicalIm,
bool saveOrig=true,
- int combine=2, bool omit=true, bool smoothen=true,
+ int combine=3, bool omit=true, bool smoothen=true,
int heuristic=1);
~Homology();
diff --git a/Plugin/HomologyComputation.cpp b/Plugin/HomologyComputation.cpp
index 5cde0ceb8883d8b7723abdc4e3a386825676dfe8..78d0ba427c2e273c5b50e4583bfcc84654d6ae4b 100644
--- a/Plugin/HomologyComputation.cpp
+++ b/Plugin/HomologyComputation.cpp
@@ -24,7 +24,7 @@ StringXNumber HomologyComputationOptions_Number[] = {
{GMSH_FULLRC, "CohomologyPhysicalGroupsBegin", NULL, -1.},
{GMSH_FULLRC, "CreatePostProcessingViews", NULL, 1.},
{GMSH_FULLRC, "ReductionOmit", NULL, 1.},
- {GMSH_FULLRC, "ReductionCombine", NULL, 2.},
+ {GMSH_FULLRC, "ReductionCombine", NULL, 3.},
{GMSH_FULLRC, "PostProcessSimplify", NULL, 1.},
{GMSH_FULLRC, "ReductionHeuristic", NULL, 1.}
};
diff --git a/Plugin/HomologyPostProcessing.cpp b/Plugin/HomologyPostProcessing.cpp
index fd0d912896799a655f938bcdb11959441b025d01..aa9f95bf1ddd2e5366fb6179ba8d6781c17d7908 100644
--- a/Plugin/HomologyPostProcessing.cpp
+++ b/Plugin/HomologyPostProcessing.cpp
@@ -26,6 +26,7 @@ StringXString HomologyPostProcessingOptions_String[] = {
{GMSH_FULLRC, "PhysicalGroupsOfOperatedChains", NULL, "1, 2"},
{GMSH_FULLRC, "PhysicalGroupsOfOperatedChains2", NULL, ""},
{GMSH_FULLRC, "PhysicalGroupsToTraceResults", NULL, ""},
+ {GMSH_FULLRC, "PhysicalGroupsToProjectResults", NULL, ""},
{GMSH_FULLRC, "NameForResultChains", NULL, "c"},
};
@@ -54,7 +55,8 @@ std::string GMSH_HomologyPostProcessingPlugin::getHelp() const
"Results a new basis for the homology space such that the incidence matrix of the new basis and the basis of the cohomology space is the identity matrix.\n\n"
"Options:\n"
- "'PhysicalGroupsToTraceResults': Trace the resulting cochains to the given physical groups.\n"
+ "'PhysicalGroupsToTraceResults': Trace the resulting (co)chains to the given physical groups.\n"
+ "'PhysicalGroupsToProjectResults': Project the resulting (co)chains to the complement of the given physical groups.\n"
"'NameForResultChains': Post-processing view name prefix for the results.\n"
"'ApplyBoundaryOperatorToResults': Apply boundary operator to the resulting chains.\n";
@@ -139,8 +141,9 @@ 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 cname = HomologyPostProcessingOptions_String[5].def;
std::string traceString = HomologyPostProcessingOptions_String[3].def;
+ std::string projectString = HomologyPostProcessingOptions_String[4].def;
int bd = (int)HomologyPostProcessingOptions_Number[0].def;
GModel* m = GModel::current();
@@ -216,6 +219,8 @@ PView *GMSH_HomologyPostProcessingPlugin::execute(PView *v)
std::vector<int> tracePhysicals;
if(!parseStringOpt(3, tracePhysicals)) return 0;
+ std::vector<int> projectPhysicals;
+ if(!parseStringOpt(4, projectPhysicals)) return 0;
std::vector<Chain<int> > curBasis;
for(unsigned int i = 0; i < basisPhysicals.size(); i++){
@@ -285,8 +290,15 @@ PView *GMSH_HomologyPostProcessingPlugin::execute(PView *v)
dim, traceString.c_str());
for(unsigned int i = 0; i < newBasis.size(); i++)
newBasis.at(i) = newBasis.at(i).getTrace(m, tracePhysicals);
- ElemChain::clearVertexCache();
}
+ if(!projectPhysicals.empty()) {
+ Msg::Info("Taking projection of result %d-chains to the complement of the domain %s", dim, projectString.c_str());
+ for(unsigned int i = 0; i < newBasis.size(); i++)
+ newBasis.at(i) = newBasis.at(i).getProject(m, projectPhysicals);
+ }
+ if(!tracePhysicals.empty() || !projectPhysicals.empty())
+ ElemChain::clearVertexCache();
+
for(unsigned int i = 0; i < newBasis.size(); i++) {
std::string dims = convertInt(newBasis.at(i).getDim());