diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index 93f4e3237f7ace378b7cafe837846ca80d3abc09..a77267522e7aff16bd465c3d1d75967f75817624 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -441,12 +441,8 @@ gmp_matrix* ChainComplex::getBasis(int dim, int basis)
void ChainComplex::getBasisChain(std::map<Cell*, int, Less_Cell>& chain,
int num, int dim, int basis, bool deform)
{
- gmp_matrix* basisMatrix;
- if(basis == 0) basisMatrix = getBasis(dim, 0);
- else if(basis == 1) basisMatrix = getBasis(dim, 1);
- else if(basis == 2) basisMatrix = getBasis(dim, 2);
- else if(basis == 3) basisMatrix = getBasis(dim, 3);
- else return;
+ if(basis < 0 || basis > 3) return;
+ gmp_matrix* basisMatrix = getBasis(dim, basis);
chain.clear();
if(dim < 0 || dim > 4) return;
@@ -473,8 +469,8 @@ void ChainComplex::getBasisChain(std::map<Cell*, int, Less_Cell>& chain,
std::map<Cell*, int, Less_Cell > subCells;
cell->getCells(subCells);
for(Cell::citer it = subCells.begin(); it != subCells.end(); it++){
- Cell* subCell = (*it).first;
- int coeff = (*it).second;
+ Cell* subCell = it->first;
+ int coeff = it->second;
chain[subCell] = coeff*elemi*torsion;
}
}
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index f17915de01562d9412bbe0f4e5081935998fed21..e21bdffaa58128dab1a09d69490676fe2bb1d4a6 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -984,9 +984,15 @@ void GModel::storeChain(int dim,
std::map<int, std::vector<MElement*> > &entityMap,
std::map<int, std::map<int, std::string> > &physicalMap)
{
- // create new discrete entities that have no associated MVertices
_storeElementsInEntities(entityMap);
_storePhysicalTagsInEntities(dim, physicalMap);
+ std::map<int, std::vector<MElement*> >::iterator it;
+ for(it = entityMap.begin(); it != entityMap.end(); it++) {
+ if(dim == 0) _chainVertices.insert(getVertexByTag(it->first));
+ else if(dim == 1) _chainEdges.insert(getEdgeByTag(it->first));
+ else if(dim == 2) _chainFaces.insert(getFaceByTag(it->first));
+ else if(dim == 3) _chainRegions.insert(getRegionByTag(it->first));
+ }
}
template<class T>
@@ -1057,11 +1063,11 @@ void GModel::_storeElementsInEntities(std::map< int, std::vector<MElement* > >&
}
template<class T>
-static void _associateEntityWithElementVertices(GEntity *ge, std::vector<T*> &elements)
+static void _associateEntityWithElementVertices(GEntity *ge, std::vector<T*> &elements, bool force=false)
{
for(unsigned int i = 0; i < elements.size(); i++){
for(int j = 0; j < elements[i]->getNumVertices(); j++){
- if (!elements[i]->getVertex(j)->onWhat() ||
+ if (force || !elements[i]->getVertex(j)->onWhat() ||
elements[i]->getVertex(j)->onWhat()->dim() > ge->dim())
elements[i]->getVertex(j)->setEntity(ge);
}
@@ -1137,6 +1143,25 @@ void GModel::_pruneMeshVertexAssociations()
entities[i]->mesh_vertices.clear();
}
_associateEntityWithMeshVertices();
+ // associate mesh vertices primarily with chain entities
+ for(riter it = _chainRegions.begin(); it != _chainRegions.end(); ++it){
+ _associateEntityWithElementVertices(*it, (*it)->tetrahedra, true);
+ _associateEntityWithElementVertices(*it, (*it)->hexahedra, true);
+ _associateEntityWithElementVertices(*it, (*it)->prisms, true);
+ _associateEntityWithElementVertices(*it, (*it)->pyramids, true);
+ _associateEntityWithElementVertices(*it, (*it)->polyhedra, true);
+ }
+ for(fiter it = _chainFaces.begin(); it != _chainFaces.end(); ++it){
+ _associateEntityWithElementVertices(*it, (*it)->triangles, true);
+ _associateEntityWithElementVertices(*it, (*it)->quadrangles, true);
+ _associateEntityWithElementVertices(*it, (*it)->polygons, true);
+ }
+ for(eiter it = _chainEdges.begin(); it != _chainEdges.end(); ++it){
+ _associateEntityWithElementVertices(*it, (*it)->lines, true);
+ }
+ for(viter it = _chainVertices.begin(); it != _chainVertices.end(); ++it){
+ _associateEntityWithElementVertices(*it, (*it)->points, true);
+ }
_storeVerticesInEntities(vertices);
}
@@ -1739,7 +1764,7 @@ void GModel::createTopologyFromFaces(std::vector<discreteFace*> &discFaces, int
for (std::vector<discreteFace*>::iterator it = discFaces.begin();
it != discFaces.end(); it++)
(*it)->findEdges(map_edges);
-
+
// return if no boundary edges (torus, sphere, ...)
if (map_edges.empty()) return;
diff --git a/Geo/GModel.h b/Geo/GModel.h
index f19bc97bd12e719e3cf10d4e52e714b252633617..7ab64098818e05519563901e365ba9781af65294 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -38,6 +38,10 @@ class GModel
private:
friend class OCCFactory;
std::multimap<std::pair<std::vector<int>, std::vector<int> >, std::string> _homologyRequests;
+ std::set<GRegion*, GEntityLessThan> _chainRegions;
+ std::set<GFace*, GEntityLessThan> _chainFaces;
+ std::set<GEdge*, GEntityLessThan> _chainEdges;
+ std::set<GVertex*, GEntityLessThan> _chainVertices;
protected:
// the name of the model
@@ -376,12 +380,12 @@ class GModel
// mesh the model
int mesh(int dimension);
- // adapt the mesh anisotropically using metrics that are computed from a set of functions f(x,y,z).
+ // adapt the mesh anisotropically using metrics that are computed from a set of functions f(x,y,z).
// For all algorithms
// parameters[1] = lcmin (default : in global gmsh options CTX::instance()->mesh.lcMin)
// parameters[2] = lcmax (default : in global gmsh options CTX::instance()->mesh.lcMax)
- // niter is the maximum number of iterations
- // Available algorithms:
+ // niter is the maximum number of iterations
+ // Available algorithms:
// 1) Assume that the function is a levelset -> adapt using Coupez technique (technique = 1)
// parameters[0] = thickness of the interface (mandatory)
// 2) Assume that the function is a physical quantity -> adapt using the Hessian (technique = 2)
@@ -397,7 +401,7 @@ class GModel
// parameters[4] = thickness of the interface in the negative ls direction (=parameters[0] if not specified)
// parameters[3] = the required minimum number of elements to represent a circle - used for curvature-based metric (default: = 15)
// 5) Same as 4, except that the transition in band E uses linear interpolation of h, instead of linear interpolation of metric
- // The algorithm first generate a mesh if no one is available
+ // The algorithm first generate a mesh if no one is available
// In this first attempt, only the highest dimensional mesh is adapted, which is ok if
// we assume that boundaries are already adapted.
@@ -487,15 +491,7 @@ class GModel
void storeChain(int dim,
std::map<int, std::vector<MElement*> > &entityMap,
std::map<int, std::map<int, std::string> > &physicalMap);
- void storeChain(std::vector<MVertex*> &vertices, int dim,
- std::map<int, std::vector<MElement*> > &entityMap,
- std::map<int, std::map<int, std::string> > &physicalMap)
- {
- _storeVerticesInEntities(vertices);
- _storeElementsInEntities(entityMap);
- _storePhysicalTagsInEntities(dim, physicalMap);
- _associateEntityWithMeshVertices();
- }
+
// request homology computation
void addHomologyRequest(const std::string &type, std::vector<int> &domain,
std::vector<int> &subdomain);
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index 9596a1715d9b90d05194100b59259ce72c1473f5..fe390b4e07cc08d9b8f6b79bbf84ef696694098b 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -130,7 +130,6 @@ void Homology::findHomologyBasis(CellComplex* cellComplex)
_subdomainEntities);
ownComplex = true;
}
- std::string domainString = _getDomainString(_domain, _subdomain);
Msg::StatusBar(2, true, "Reducing cell complex...");
@@ -151,7 +150,7 @@ void Homology::findHomologyBasis(CellComplex* cellComplex)
Msg::StatusBar(2, true, "Done computing homology space basis (%g s)", t2 - t1);
int dim = cellComplex->getDim();
-
+ std::string domain = _getDomainString(_domain, _subdomain);
int HRank[4];
for(int j = 0; j < 4; j++){
HRank[j] = 0;
@@ -162,7 +161,7 @@ void Homology::findHomologyBasis(CellComplex* cellComplex)
std::string generator = "";
convert(i, generator);
- std::string name = "H" + dimension + domainString + generator;
+ std::string name = "H" + dimension + domain + generator;
std::map<Cell*, int, Less_Cell> chain;
chainComplex.getBasisChain(chain, i, j, 3, _smoothen);
int torsion = chainComplex.getTorsion(j,i);
@@ -181,7 +180,7 @@ void Homology::findHomologyBasis(CellComplex* cellComplex)
if(ownComplex) delete cellComplex;
- Msg::Info("Ranks of homology spaces:");
+ Msg::Info("Ranks of domain %shomology spaces:", domain.c_str());
Msg::Info("H0 = %d", HRank[0]);
Msg::Info("H1 = %d", HRank[1]);
Msg::Info("H2 = %d", HRank[2]);
@@ -221,7 +220,7 @@ void Homology::findCohomologyBasis(CellComplex* cellComplex)
Msg::StatusBar(2, true, "Done computing cohomology space basis (%g s)", t2- t1);
int dim = cellComplex->getDim();
-
+ std::string domain = _getDomainString(_domain, _subdomain);
int HRank[4];
for(int i = 0; i < 4; i++) HRank[i] = 0;
for(int j = 3; j > -1; j--){
@@ -234,7 +233,7 @@ void Homology::findCohomologyBasis(CellComplex* cellComplex)
convert(i, generator);
std::string name = "H" + dimension + "*" +
- _getDomainString(_domain, _subdomain) + generator;
+ domain + generator;
std::map<Cell*, int, Less_Cell> chain;
chainComplex.getBasisChain(chain, i, j, 3, _smoothen);
int torsion = chainComplex.getTorsion(j,i);
@@ -252,7 +251,7 @@ void Homology::findCohomologyBasis(CellComplex* cellComplex)
if(ownComplex) delete cellComplex;
- Msg::Info("Ranks of cohomology spaces:");
+ Msg::Info("Ranks of domain %scohomology spaces:", domain.c_str());
Msg::Info("H0* = %d", HRank[0]);
Msg::Info("H1* = %d", HRank[1]);
Msg::Info("H2* = %d", HRank[2]);
@@ -313,9 +312,6 @@ void Homology::_createPhysicalGroup(std::map<Cell*, int, Less_Cell>& chain, std:
MElementFactory factory;
int dim = 0;
- std::map<Cell*, int, Less_Cell> chain2;
- std::string name2 = name + "i";
-
typedef std::map<Cell*, int, Less_Cell>::iterator citer;
for(citer cit = chain.begin(); cit != chain.end(); cit++){
Cell* cell = cit->first;
@@ -338,15 +334,8 @@ void Homology::_createPhysicalGroup(std::map<Cell*, int, Less_Cell>& chain, std:
elements.push_back(ecopy);
}
- //std::vector<double> coeffs (1,1);
std::vector<double> coeffs (1,abs(coeff));
data[e->getNum()] = coeffs;
-
- /*if(abs(coeff) > 1) {
- if(coeff < 0) chain2[cell] = coeff+1;
- else chain2[cell] = coeff-1;
- }*/
-
}
GModel* m = this->getModel();
@@ -383,7 +372,6 @@ void Homology::_createPhysicalGroup(std::map<Cell*, int, Less_Cell>& chain, std:
view->setOptions(opt);
#endif
}
- //if(!chain2.empty()) _createPhysicalGroup(chain2, name2);
}
bool Homology::writeBasisMSH(bool binary)
diff --git a/Geo/Homology.h b/Geo/Homology.h
index a781b714a5c0bebf0804f3c58b884c0548a75150..c5dbdc9e51c663939edb175b91271f7c37906411 100644
--- a/Geo/Homology.h
+++ b/Geo/Homology.h
@@ -52,7 +52,7 @@ class Homology
// create a string describing the generator
std::string _getDomainString(const std::vector<int>& domain,
- const std::vector<int>& subdomain);
+ const std::vector<int>& subdomain);
// remove cells with coefficient 0
int _eraseNullCells(std::map<Cell*, int, Less_Cell>& chain);