pruneMeshVertexAssociations, no not save 0- and highest dim-chains

Geo/GModel.cpp

@@ -969,14 +969,7 @@ void GModel::storeChain(int dim,
 {
   // create new discrete entities that have no associated MVertices
   _storeElementsInEntities(entityMap);
-   // give them physical tags
   _storePhysicalTagsInEntities(dim, physicalMap);
-  // For some weird reason, mesh element color by elementary/physical
-  // entity causes crashes when drawing mesh elements of the new discrete
-  // entities
-  // (somehow related to MVertex.onWhat() as MVertices are unaware of the
-  //  new discrete entities)
-  CTX::instance()->mesh.colorCarousel = 0; // color by element type
 }
 
 template<class T>
@@ -1113,6 +1106,23 @@ void GModel::_storeVerticesInEntities(std::vector<MVertex*> &vertices)
   }
 }
 
+void GModel::_pruneMeshVertexAssociations()
+{
+  std::vector<GEntity*> entities;
+  std::vector<MVertex*> vertices;
+  getEntities(entities);
+  for(unsigned int i = 0; i < entities.size(); i++) {
+    for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++) {
+      MVertex* v = entities[i]->mesh_vertices[j];
+      v->setEntity(0);
+      vertices.push_back(v);
+    }
+    entities[i]->mesh_vertices.clear();
+  }
+  _associateEntityWithMeshVertices();
+  _storeVerticesInEntities(vertices);
+}
+
 void GModel::checkMeshCoherence(double tolerance)
 {
   int numEle = getNumMeshElements();
@@ -2752,7 +2762,7 @@ void GModel::computeHomology()
               std::multimap<dpair, std::string>::iterator> itp =
       _homologyRequests.equal_range(*it);
     Homology* homology = new Homology(this, itp.first->first.first,
-                                      itp.first->first.second);
+                                      itp.first->first.second, false);
     CellComplex *cellcomplex = homology->createCellComplex();
     if(cellcomplex->getSize(0)){
       for(std::multimap<dpair, std::string>::iterator itt = itp.first;
@@ -2767,6 +2777,7 @@ void GModel::computeHomology()
         else
           Msg::Error("Unknown type of homology computation: %s", type.c_str());
       }
+      _pruneMeshVertexAssociations();
     }
     delete cellcomplex;
     delete homology;

Geo/GModel.h

@@ -102,6 +102,12 @@ class GModel
   void _storeVerticesInEntities(std::map<int, MVertex*> &vertices);
   void _storeVerticesInEntities(std::vector<MVertex*> &vertices);
 
+  // remove all mesh vertex associations to geometrical entities and
+  // remove vertices from geometrical entities, then
+  // _associateEntityWithMeshVertices and _storeVerticesInEntities
+  // are called to rebuild the associations
+  void _pruneMeshVertexAssociations();
+
   // store the physical tags in the geometrical entities
   void _storePhysicalTagsInEntities(int dim,
                                     std::map<int, std::map<int, std::string> > &map);

Geo/Homology.cpp

@@ -15,10 +15,10 @@
 #if defined(HAVE_KBIPACK)
 
 Homology::Homology(GModel* model, std::vector<int> physicalDomain,
-		   std::vector<int> physicalSubdomain,
+		   std::vector<int> physicalSubdomain, bool save0N,
 		   bool combine, bool omit, bool smoothen) :
   _model(model), _domain(physicalDomain), _subdomain(physicalSubdomain),
-  _combine(combine), _omit(omit), _smoothen(smoothen)
+  _save0N(save0N), _combine(combine), _omit(omit), _smoothen(smoothen)
 {
   _fileName = "";
 
@@ -137,6 +137,8 @@ void Homology::findGenerators(CellComplex* cellComplex)
   t2 = Cpu();
   Msg::StatusBar(2, true, "Done computing homology spaces (%g s)", t2 - t1);
 
+  int dim = cellComplex->getDim();
+
   int HRank[4];
   for(int j = 0; j < 4; j++){
     HRank[j] = 0;
@@ -165,7 +167,7 @@ void Homology::findGenerators(CellComplex* cellComplex)
       // FIXME: Cell* and CellComplex* pointers should not outlive
       // the objects, don't store Chains containing them for now
       //_basisChains[chain->createPGroup()] = chain;
-      chain->createPGroup();
+      if(!_save0N && (j != 0 && j != dim)) chain->createPGroup();
       delete chain;
     }
   }
@@ -246,7 +248,7 @@ void Homology::findDualGenerators(CellComplex* cellComplex)
       // FIXME: Cell* and CellComplex* pointers should not outlive
       // the objects, don't store Chains containing them for now
       //_basisChains[chain->createPGroup()] = chain;
-      chain->createPGroup();
+      if(!_save0N && (j != 0 && j != dim)) chain->createPGroup();
       delete chain;
     }
   }
@@ -266,114 +268,6 @@ void Homology::findDualGenerators(CellComplex* cellComplex)
 		 HRank[0], HRank[1], HRank[2], HRank[3]);
 }
 
-/*void Homology::findHomSequence(){
-  CellComplex* cellComplex = createCellComplex(_domainEntities,
-					       _subdomainEntities);
-  Msg::StatusBar(2, true, "Reducing cell complex...");
-
-  double t1 = Cpu();
-  cellComplex->reduceComplex();
-  double t2 = Cpu();
-
-  Msg::StatusBar(2, true, "Done reducing cell complex (%g s)", t2 - t1);
-  Msg::Info("%d volumes, %d faces, %d edges and %d vertices",
-            cellComplex->getSize(3), cellComplex->getSize(2),
-	    cellComplex->getSize(1), cellComplex->getSize(0));
-
-  Msg::StatusBar(2, true, "Computing homology spaces...");
-  t1 = Cpu();
-
-  ChainComplex* subcomplex = new ChainComplex(cellComplex, 2);
-  ChainComplex* complex = new ChainComplex(cellComplex, 1);
-  ChainComplex* relcomplex = new ChainComplex(cellComplex, 0);
-
-  subcomplex->computeHomology();
-  complex->computeHomology();
-  relcomplex->computeHomology();
-
-  t2 = Cpu();
-  Msg::StatusBar(2, true, "Done compuring homology spaces (%g s)", t2 - t1);
-
-  Msg::StatusBar(2, true, "Computing homology sequence...");
-  HomologySequence* seq = new HomologySequence(subcomplex,
-					       complex, relcomplex);
-  t1 = Cpu();
-  Msg::StatusBar(2, true, "Done computing homology sequence (%g s)", t1 - t2);
-
-  for(int task = 0; task < 3; task++){
-    ChainComplex* chains;
-
-    std::string domainString = "";
-    std::vector<int> empty;
-    if(task == 0) {
-      chains = subcomplex;
-      domainString = getDomainString(_subdomain, empty);
-    }
-    else if(task == 1) {
-      chains = complex;
-      domainString = getDomainString(_domain, empty);
-    }
-    else{
-      chains = relcomplex;
-      domainString = getDomainString(_domain, _subdomain);
-    }
-
-    int HRank[4];
-    for(int j = 0; j < 4; j++){
-      HRank[j] = 0;
-      std::string dimension = "";
-      convert(j, dimension);
-      for(int i = 1; i <= chains->getBasisSize(j, 3); i++){
-
-	std::string generator = "";
-	convert(i, generator);
-
-	std::string name = "H" + dimension + domainString + generator;
-	std::map<Cell*, int, Less_Cell> protoChain;
-	chains->getBasisChain(protoChain, i, j, 3, true);
-	Chain* chain = new Chain(protoChain, cellComplex, _model, name,
-				 chains->getTorsion(j,i));
-	if(chain->getSize() == 0) {
-	  delete chain;
-	  continue;
-	}
-	HRank[j] = HRank[j] + 1;
-	if(chain->getTorsion() != 1){
-	  Msg::Warning("H%d %d has torsion coefficient %d!",
-		       j, i, chain->getTorsion());
-	}
-	_generators.push_back(chain->createPGroup());
-	delete chain;
-      }
-
-    }
-    if(task == 0){
-      Msg::Info("Ranks of relative homology spaces:");
-    }
-    if(task == 1){
-      Msg::Info("Ranks of absolute homology spaces:");
-    }
-    if(task == 2){
-      Msg::Info("Ranks of absolute homology spaces of relative subdomain:");
-    }
-    Msg::Info("H0 = %d", HRank[0]);
-    Msg::Info("H1 = %d", HRank[1]);
-    Msg::Info("H2 = %d", HRank[2]);
-    Msg::Info("H3 = %d", HRank[3]);
-
-    Msg::StatusBar(2, false, "H0: %d, H1: %d, H2: %d, H3: %d",
-		   HRank[0], HRank[1], HRank[2], HRank[3]);
-  }
-
-  delete subcomplex;
-  delete complex;
-  delete relcomplex;
-
-  if(_fileName != "") writeGeneratorsMSH();
-  delete cellComplex;
-  delete seq;
-}
-*/
 std::string Homology::getDomainString(const std::vector<int>& domain,
 				      const std::vector<int>& subdomain)
 {

Geo/Homology.h

@@ -50,6 +50,9 @@ class Homology
   // use chain smoothning
   bool _smoothen;
 
+  // save chains of 0 and highest dimension
+  bool _save0N;
+
   int _maxdomain;
   int _maxnum;
 
@@ -61,7 +64,7 @@ class Homology
  public:
 
   Homology(GModel* model, std::vector<int> physicalDomain,
-	   std::vector<int> physicalSubdomain,
+	   std::vector<int> physicalSubdomain, bool save0N=false,
 	   bool combine=true, bool omit=true, bool smoothen=true);
   ~Homology();