diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp index 90c2c9f7fc8a67b2f1119514c93d9fd1c4c194fb..53d0dce5e9e8f5552505669f1b02442270dc9097 100644 --- a/Geo/ChainComplex.cpp +++ b/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,28 +568,25 @@ 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); } } return false; } -void Chain::smoothenChain(){ - +void Chain::smoothenChain(){ if(!_cellComplex->simplicial()) return; int start = getSize(); @@ -656,9 +664,16 @@ void Chain::createPView(){ int coeff = (*cit).second; std::vector<MVertex*> v = cell->getVertexVector(); if(cell->getDim() > 0 && coeff < 0){ // flip orientation - MVertex* temp = v[0]; - v[0] = v[1]; - v[1] = temp; + 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); diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp index 5b81afce1e155806ecd471715a27849b858ec2b7..9d83783a416f046b7f54c306bdfd594edc5bb2b9 100644 --- a/Geo/Homology.cpp +++ b/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]); diff --git a/Geo/Homology.h b/Geo/Homology.h index b7212549cb582335ce99d3e41f69d2b9ff921f82..99af3629f14a3a693deafdf7428c94b7182995f0 100644 --- a/Geo/Homology.h +++ b/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; } diff --git a/tutorial/t10.geo b/tutorial/t10.geo index bd92364cc00540b73522d3059918daa5d1d18c70..464d688417c380680ea017f3403f3950a28e3c41 100644 --- a/tutorial/t10.geo +++ b/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}};