diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 3d59c7fac1631144aec101d897d60767d7cb5149..7ea6d4a226e8254f76003e76d4e79ce2a702764b 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -11,6 +11,8 @@
CellComplex::CellComplex(GModel* model,
std::vector<MElement*>& domainElements,
std::vector<MElement*>& subdomainElements,
+ std::vector<MElement*>& nondomainElements,
+ std::vector<MElement*>& nonsubdomainElements,
std::vector<MElement*>& immuneElements,
bool saveOriginalComplex) :
_model(model), _dim(0), _simplicial(true), _saveorig(saveOriginalComplex)
@@ -18,6 +20,8 @@ CellComplex::CellComplex(GModel* model,
_deleteCount = 0;
_insertCells(subdomainElements, 1);
_insertCells(domainElements, 0);
+ _removeCells(nonsubdomainElements, 1);
+ _removeCells(nondomainElements, 0);
_immunizeCells(immuneElements);
int num = 0;
for(int dim = 0; dim < 4; dim++){
@@ -78,6 +82,53 @@ bool CellComplex::_insertCells(std::vector<MElement*>& elements,
return true;
}
+bool CellComplex::_removeCells(std::vector<MElement*>& elements,
+ int domain)
+{
+
+ std::set<Cell*, Less_Cell> removed[4];
+
+ for(unsigned int i=0; i < elements.size(); i++){
+ MElement* element = elements.at(i);
+ int type = element->getType();
+ if(type == TYPE_PYR || type == TYPE_PRI ||
+ type == TYPE_POLYG || type == TYPE_POLYH) {
+ Msg::Error("Mesh element type %d not implemented in homology solver", type);
+ return false;
+ }
+ Cell* cell = new Cell(element, domain);
+ int dim = cell->getDim();
+ citer cit = _cells[dim].find(cell);
+ if(cit != lastCell(dim)) {
+ removeCell(*cit);
+ removed[dim].insert(cell);
+ }
+ else delete cell;
+ }
+
+ for (int dim = 3; dim > 0; dim--){
+ for(citer cit = removed[dim].begin(); cit != removed[dim].end(); cit++){
+ Cell* cell = *cit;
+ for(int i = 0; i < cell->getNumBdElements(); i++){
+ Cell* newCell = new Cell(cell, i);
+
+ citer cit2 = _cells[dim-1].find(newCell);
+ if(cit2 != lastCell(dim-1)) {
+ removeCell(*cit2);
+ removed[dim-1].insert(newCell);
+ }
+ else delete newCell;
+ }
+ }
+ }
+ for (int dim = 3; dim > 0; dim--){
+ for(citer cit = removed[dim].begin(); cit != removed[dim].end(); cit++){
+ delete *cit;
+ }
+ }
+ return true;
+}
+
bool CellComplex::_immunizeCells(std::vector<MElement*>& elements)
{
for(unsigned int i=0; i < elements.size(); i++){
@@ -291,7 +342,7 @@ int CellComplex::reduceComplex(bool docombine, bool omit)
std::vector<Cell*> empty;
for(int i = 3; i > 0; i--) count = count + reduction(i, false, empty);
- if(omit && !count){
+ if(omit){
removeSubdomain();
std::vector<Cell*> newCells;
@@ -362,7 +413,8 @@ int CellComplex::coreduceComplex(bool docombine, bool omit)
}
}
- if(omit && !count){
+ if(omit){
+
std::vector<Cell*> newCells;
while (getSize(0) != 0){
citer cit = firstCell(0);
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index c6d4a76bf48d3ba516e8e55a18526ee636634458..326521e7bf47241528b9e1ac9cdcf3b3534c26bd 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -48,6 +48,7 @@ class CellComplex
// for constructor
bool _insertCells(std::vector<MElement*>& elements, int domain);
+ bool _removeCells(std::vector<MElement*>& elements, int domain);
bool _immunizeCells(std::vector<MElement*>& elements);
@@ -67,6 +68,8 @@ class CellComplex
CellComplex(GModel* model,
std::vector<MElement*>& domainElements,
std::vector<MElement*>& subdomainElements,
+ std::vector<MElement*>& nondomainElements,
+ std::vector<MElement*>& nonsubdomainElements,
std::vector<MElement*>& immuneElements,
bool saveOriginalComplex=true);
~CellComplex();
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index 9f9b8ccfd8f063959ebc283943de3bbd1cc4774b..e163d8aad3c5d31835b1ea624aad04faf7f502a6 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -36,41 +36,43 @@ Homology::Homology(GModel* model, std::vector<int> physicalDomain,
}
}
+ _getEntities(_domain, _domainEntities);
+ _getEntities(_subdomain, _subdomainEntities);
+ _getEntities(_nondomain, _nondomainEntities);
+ _getEntities(_nonsubdomain, _nonsubdomainEntities);
+ _getEntities(_imdomain, _immuneEntities);
+
+}
+
+void Homology::_getEntities(const std::vector<int>& physicalGroups,
+ std::vector<GEntity*>& entities)
+{
+ entities.clear();
std::map<int, std::vector<GEntity*> > groups[4];
- model->getPhysicalGroups(groups);
+ _model->getPhysicalGroups(groups);
std::map<int, std::vector<GEntity*> >::iterator it;
- for(unsigned int i = 0; i < _domain.size(); i++){
- for(int j = 0; j < 4; j++){
- it = groups[j].find(_domain.at(i));
- if(it != groups[j].end()){
- std::vector<GEntity*> physicalGroup = (*it).second;
- for(unsigned int k = 0; k < physicalGroup.size(); k++){
- _domainEntities.push_back(physicalGroup.at(k));
- }
- }
- }
- }
- for(unsigned int i = 0; i < _subdomain.size(); i++){
+ for(unsigned int i = 0; i < physicalGroups.size(); i++){
for(int j = 0; j < 4; j++){
- it = groups[j].find(_subdomain.at(i));
+ it = groups[j].find(physicalGroups.at(i));
if(it != groups[j].end()){
std::vector<GEntity*> physicalGroup = (*it).second;
for(unsigned int k = 0; k < physicalGroup.size(); k++){
- _subdomainEntities.push_back(physicalGroup.at(k));
+ entities.push_back(physicalGroup.at(k));
}
}
}
}
- for(unsigned int i = 0; i < _imdomain.size(); i++){
- for(int j = 0; j < 4; j++){
- it = groups[j].find(_imdomain.at(i));
- if(it != groups[j].end()){
- std::vector<GEntity*> physicalGroup = (*it).second;
- for(unsigned int k = 0; k < physicalGroup.size(); k++){
- _immuneEntities.push_back(physicalGroup.at(k));
- }
- }
+}
+
+void Homology::_getElements(const std::vector<GEntity*>& entities,
+ std::vector<MElement*>& elements)
+{
+ elements.clear();
+ for(unsigned int j=0; j < entities.size(); j++) {
+ for(unsigned int i=0; i < entities.at(j)->getNumMeshElements(); i++){
+ MElement* element = entities.at(j)->getMeshElement(i);
+ elements.push_back(element);
}
}
}
@@ -85,32 +87,22 @@ void Homology::_createCellComplex()
std::vector<MElement*> domainElements;
std::vector<MElement*> subdomainElements;
+ std::vector<MElement*> nondomainElements;
+ std::vector<MElement*> nonsubdomainElements;
std::vector<MElement*> immuneElements;
- for(unsigned int j=0; j < _domainEntities.size(); j++) {
- for(unsigned int i=0; i < _domainEntities.at(j)->getNumMeshElements(); i++){
- MElement* element = _domainEntities.at(j)->getMeshElement(i);
- domainElements.push_back(element);
- }
- }
- for(unsigned int j=0; j < _subdomainEntities.size(); j++) {
- for(unsigned int i=0; i < _subdomainEntities.at(j)->getNumMeshElements();
- i++){
- MElement* element = _subdomainEntities.at(j)->getMeshElement(i);
- subdomainElements.push_back(element);
- }
- }
- for(unsigned int j=0; j < _immuneEntities.size(); j++) {
- for(unsigned int i=0; i < _immuneEntities.at(j)->getNumMeshElements();
- i++){
- MElement* element = _immuneEntities.at(j)->getMeshElement(i);
- immuneElements.push_back(element);
- }
- }
+
+ _getElements(_domainEntities, domainElements);
+ _getElements(_subdomainEntities, subdomainElements);
+ _getElements(_nondomainEntities, nondomainElements);
+ _getElements(_nonsubdomainEntities, nonsubdomainElements);
+ _getElements(_immuneEntities, immuneElements);
if(_cellComplex != NULL) delete _cellComplex;
_cellComplex = new CellComplex(_model,
domainElements,
subdomainElements,
+ nondomainElements,
+ nonsubdomainElements,
immuneElements,
_saveOrig);
diff --git a/Geo/Homology.h b/Geo/Homology.h
index 16bf78964604f9368c937f8d5174cfd1c8142da3..a84124a951b49cc9a565d9c8e7b8965c4b61da02 100644
--- a/Geo/Homology.h
+++ b/Geo/Homology.h
@@ -30,10 +30,14 @@ class Homology
// physical group IDs
std::vector<int> _domain;
std::vector<int> _subdomain;
+ std::vector<int> _nondomain;
+ std::vector<int> _nonsubdomain;
std::vector<int> _imdomain;
// corresponding geometrical entities
std::vector<GEntity*> _domainEntities;
std::vector<GEntity*> _subdomainEntities;
+ std::vector<GEntity*> _nondomainEntities;
+ std::vector<GEntity*> _nonsubdomainEntities;
std::vector<GEntity*> _immuneEntities;
// save original cell complex
@@ -61,6 +65,11 @@ class Homology
typedef std::map<Cell*, int, Less_Cell>::iterator citer;
+ void _getEntities(const std::vector<int>& physicalGroups,
+ std::vector<GEntity*>& entities);
+ void _getElements(const std::vector<GEntity*>& entities,
+ std::vector<MElement*>& elements);
+
// create a string describing the generator
std::string _getDomainString(const std::vector<int>& domain,
const std::vector<int>& subdomain);
diff --git a/doc/texinfo/gmsh.texi b/doc/texinfo/gmsh.texi
index d938c021a78d6a50831091e1516556ca1f8b799e..86018f0c08c13e2fd69167d5ad68665137356a56 100644
--- a/doc/texinfo/gmsh.texi
+++ b/doc/texinfo/gmsh.texi
@@ -2525,6 +2525,7 @@ Show the entities listed in @var{expression-list}, if
Show the entity @var{char-expression}, if @code{General.VisibilityMode} is
set to @code{0} or @code{1} (@var{char-expression} can for example be
@code{"*"}).
+
@end ftable
@c -------------------------------------------------------------------------
@@ -3023,6 +3024,20 @@ Shows the mesh of the entity @var{char-expression}, if
Sets number of elliptic smoothing steps for the surfaces listed in
@var{expression-list} (smoothing only applies to transfinite meshes at
the moment).
+
+@item Homology @{ @{ @var{expression-list} @}, @{ @var{expression-list} @} @};
+Compute a basis representation for homology spaces after a mesh has
+been generated. The first @var{expression-list} is a list of physical
+groups that constitute the computation domain. The second
+@var{expression-list} is a list of physical groups that constitute the
+relative subdomain of relative homology computation; if empty,
+absolute homology is computed. Resulting basis representation chains
+are stored as physical groups in the mesh.
+
+@item Cohomology @{ @{ @var{expression-list} @}, @{ @var{expression-list} @} @};
+Similar to command @code{Homology}, but computes a basis representation
+for cohomology spaces instead.
+
@end ftable
@c -------------------------------------------------------------------------
diff --git a/tutorial/t14.geo b/tutorial/t14.geo
index c945f69ae5f1d5a0b40c1226956174ba4b784fe8..8af18545f5cbf1d7352bed0c487d035dfc82f378 100644
--- a/tutorial/t14.geo
+++ b/tutorial/t14.geo
@@ -1,19 +1,15 @@
-/*********************************************************************
+/*********************************************************************
*
* Gmsh tutorial 14
*
- * Homology computation
+ * Homology and cohomology computation
*
*********************************************************************/
-
+
// Homology computation in Gmsh finds representative chains of
-// (relative) homology spaces using a mesh of a model. Those
-// representatives generate the (relative) homology spaces of the
-// model.
-
-// The generator chains are stored in a given .msh-file as physical
-// groups, whose mesh elements are oriented such that their
-// coefficients are 1 in the generator chain.
+// (relative) (co)homology space bases using a mesh of a model.
+// The representative basis chains are stored in the mesh as
+// physical groups of Gmsh, one for each chain.
// Create an example geometry
@@ -40,7 +36,7 @@ Plane Surface(15) = {13, 14};
Extrude {0, 0, h}{ Surface{15}; }
// Create physical groups, which are used to define the domain of the
-// homology computation and the subdomain of the relative homology
+// (co)homology computation and the subdomain of the relative (co)homology
// computation.
// Whole domain
@@ -60,17 +56,20 @@ Physical Surface(80) = bnd[];
bnd[] -= {36, 44, 52, 60};
Physical Surface(75) = bnd[];
-// Find generators of relative homology spaces of the domain modulo
-// the four terminals.
+// Find bases for relative homology spaces of
+// the domain modulo the four terminals.
Homology {{1}, {70, 71, 72, 73}};
-// Find the corresponding thick cuts
-Cohomology {{1}, {70, 71, 72, 73}};
-
-// Find the corresponding thin cuts, generators of relative homology
-// spaces modulo the non-terminal domain surface.
+// Find homology space bases isomorphic to the previous bases:
+// homology spaces modulo the non-terminal domain surface,
+// a.k.a the thin cuts.
Homology {{1}, {75}};
+// Find cohomology space bases isomorphic to the previous bases:
+// cohomology spaces of the domain modulo the four terminals,
+// a.k.a the thick cuts.
+Cohomology {{1}, {70, 71, 72, 73}};
+
// More examples:
// Homology {1};
// Homology;