Update for chain deformation orientation bug. Enabled some gmsh options by

default after homology computatation for better generator chain visualization. 

Geo/ChainComplex.cpp

@@ -8,6 +8,7 @@
 #include "ChainComplex.h"
 #include "OS.h"
 #include "PView.h"
+#include "Options.h"
 
 #if defined(HAVE_KBIPACK)
 
@@ -484,34 +485,44 @@ Chain::Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs, CellComp
 }
 
 bool Chain::deform(std::map<Cell*, int, Less_Cell> &cellsInChain, std::map<Cell*, int, Less_Cell> &cellsNotInChain){
-  
+  //printf("--- \n");
   std::vector<int> cc;
   std::vector<int> bc;
   
   for(citer cit = cellsInChain.begin(); cit != cellsInChain.end(); cit++){
     Cell* c = (*cit).first;
     c->setImmune(false);
+    //c->printCell();
     if(!c->inSubdomain()) {
       cc.push_back(getCoeff(c));
       bc.push_back((*cit).second);
     }
-    removeCell(c);
   }
   
-  // FIXME: orientations don't get right on 2-chains
-  int flip = 1;
-  if(cc[0] != bc[0]) flip = flip*-1;
+  // FIXME: orientations don't get right on 2-chains with bending or subdomain cells (disabled for now).
+  
+  bool flip = false;
+  if(cc[0] != bc[0]) flip = true; 
+  for(int i = 0; i < cc.size(); i++){
+    //printf("cc: %d, bc: %d \n", cc.at(i), bc.at(i));
+    if(flip && cc[i] == bc[i]) return false;
+    if(!flip && cc[i] != bc[i]) return false; 
+  }
+  
+  for(citer cit = cellsInChain.begin(); cit != cellsInChain.end(); cit++) removeCell((*cit).first);
   
   int n = 1;
   for(citer cit = cellsNotInChain.begin(); cit != cellsNotInChain.end(); cit++){
     Cell* c = (*cit).first;
     if(n == 2) c->setImmune(true);
     else c->setImmune(false);
-    int coeff = -1*(*cit).second*flip;
+    int coeff = -1*(*cit).second;
+    if(flip) coeff = coeff*-1;
     addCell(c, coeff);
+    //c->printCell();
     n++;
   }
-  
+  //printf("--- \n");
   
   return true;
 }
@@ -541,7 +552,7 @@ bool Chain::deformChain(std::pair<Cell*, int> cell, bool straighten, bool bend){
       Cell* c = (*cit2).first;
       int coeff = getCoeff(c);
       if(c->getImmune()) next = true;
-      if(bend && c->inSubdomain()) next = true;
+      if(/* FIXME: bend && */c->inSubdomain()) next = true;
       if(coeff > 1 || coeff < -1) next = true; 
     }
     
@@ -557,20 +568,18 @@ bool Chain::deformChain(std::pair<Cell*, int> cell, bool straighten, bool bend){
     if( (getDim() == 1 && cellsInChain.size() == 2 && cellsNotInChain.size() == 1 && straighten) ||
         (getDim() == 2 && cellsInChain.size() == 3 && cellsNotInChain.size() == 1 && straighten)){
       //printf("straighten \n");
-      deform(cellsInChain, cellsNotInChain);
-      return true;
+      return deform(cellsInChain, cellsNotInChain);
     }
     else if ( (getDim() == 1 && cellsInChain.size() == 1 && cellsNotInChain.size() == 2 && bend) ||
               (getDim() == 2 && cellsInChain.size() == 2 && cellsNotInChain.size() == 2 && bend)){
       //printf("bend \n");
-      deform(cellsInChain, cellsNotInChain);
-      return true;
+      //FIXME: return deform(cellsInChain, cellsNotInChain);
+      return false;
     }
     else if ((getDim() == 1 && cellsInChain.size() == 3 && cellsNotInChain.size() == 0) ||
              (getDim() == 2 && cellsInChain.size() == 4 && cellsNotInChain.size() == 0)){
       //printf("remove boundary \n");
-      deform(cellsInChain, cellsNotInChain);
-      return true;
+      return deform(cellsInChain, cellsNotInChain);
     }
   }
   
@@ -578,7 +587,6 @@ bool Chain::deformChain(std::pair<Cell*, int> cell, bool straighten, bool bend){
 }
 
 void Chain::smoothenChain(){
-  
   if(!_cellComplex->simplicial()) return;
   
   int start = getSize();
@@ -656,10 +664,17 @@ void Chain::createPView(){
     int coeff = (*cit).second;
     std::vector<MVertex*> v = cell->getVertexVector();
     if(cell->getDim() > 0 && coeff < 0){ // flip orientation
+      if(getDim() != 1){
+        MVertex* temp = v[1];
+        v[1] = v[v.size()-1];
+        v[v.size()-1] = temp;
+      }
+      else{
         MVertex* temp = v[0];
         v[0] = v[1];
         v[1] = temp;
       }
+    }
     MElement *e = factory.create(cell->getTypeForMSH(), v, cell->getNum(), cell->getPartition());
     for(int i = 0; i < abs(coeff); i++) elements.push_back(e);
     
@@ -683,6 +698,21 @@ void Chain::createPView(){
   physicalInfo[physicalNum]=getName();
   physicalMap[entityNum] = physicalInfo;
   
+  // hide mesh
+  opt_mesh_points(0, GMSH_SET, 0);
+  opt_mesh_lines(0, GMSH_SET, 0);
+  opt_mesh_triangles(0, GMSH_SET, 0);
+  opt_mesh_quadrangles(0, GMSH_SET, 0);
+  opt_mesh_tetrahedra(0, GMSH_SET, 0);
+  opt_mesh_hexahedra(0, GMSH_SET, 0);
+  opt_mesh_prisms(0, GMSH_SET, 0);
+  opt_mesh_pyramids(0, GMSH_SET, 0);
+
+  // show post-processing normals, tangents and element boundaries
+  opt_view_normals(0, GMSH_SET, 20);
+  opt_view_tangents(0, GMSH_SET, 20);
+  opt_view_show_element(0, GMSH_SET, 1);
+  
   if(!data.empty()){
     _model->storeChain(getDim(), entityMap, physicalMap);
     _model->setPhysicalName(getName(), getDim(), physicalNum);

Geo/Homology.cpp

@@ -147,7 +147,7 @@ void Homology::findGenerators(std::string fileName){
   }
   
   createPViews();
-  writeGeneratorsMSH(fileName);
+  if(fileName != "") writeGeneratorsMSH(fileName);
   
   Msg::Info("Ranks of homology groups for primal cell complex:");
   Msg::Info("H0 = %d", HRank[0]);
@@ -253,7 +253,7 @@ void Homology::findDualGenerators(std::string fileName){
   }
    
   createPViews();
-  writeGeneratorsMSH(fileName);
+  if(fileName != "") writeGeneratorsMSH(fileName);
   
   Msg::Info("Ranks of homology groups for dual cell complex:");
   Msg::Info("H0* = %d", HRank[0]);

Geo/Homology.h

@@ -117,8 +117,8 @@ class Homology
          if(!chain->writeChainMSH(fileName)) return false;
        }
      }*/
-     Msg::Info("Wrote results to %s.", fileName.c_str());
-     printf("Wrote results to %s. \n", fileName.c_str());
+     Msg::Info("Wrote homology computation results to %s.", fileName.c_str());
+     printf("Wrote homology computation results to %s. \n", fileName.c_str());
      
      return true;
    }

tutorial/t10.geo

@@ -130,22 +130,16 @@ Physical Surface(75) = {46, 63, 66, 52, 50, 48, 54, 60, 58, 56};
 Mesh 3;
 
 // Find generators of relative homology spaces of the domain modulo the four terminals.
-// Save the generator chains to t10_homgen.msh.
-HomGen("t10_homgen.msh") = {{69}, {70, 71, 72, 73}};
+// Save the generator chains to t10_hom.msh.
+HomGen("t10_hom.msh") = {{69}, {70, 71, 72, 73}};
 
 // Find the corresponding cuts.
-// Save the cut chains to t10_homcut.msh.
-HomCut("t10_homcut.msh") = {{69}, {70, 71, 72, 73}};
+// Save the cut chains to t10_hom.msh.
+HomCut("t10_hom.msh") = {{69}, {70, 71, 72, 73}};
 
 // Only find and print the ranks of the relative homology spaces (Betti numbers).
 HomRank {{69},{70, 71, 72, 73}};
 
-// Hide mesh elements for instant visualization
-Mesh.Points = 0;
-Mesh.Lines = 0;
-Mesh.Triangles = 0;
-Mesh.Tetrahedra = 0; 
-
 // More examples (uncomment):
 //  HomGen("t10_homgen.msh_1") = {{69}, {}}; 
 //  HomGen("t10_homgen.msh_2") = {{69}, {74}};