[Gmsh] possible bug in the Frontal Algo

[walter steffe]

Hello,

I have new information on this bug which seems to be REAL and not just
possible.

The problem arises inside of the function edgeSwapPass() called from
bowyerWatsonFrontal() (the second call of edgeSwapPass() just before the
call of transferDataStructure(gf, AllTris, DATA).

To analyze the problem I have added following line of codes (marked
below as w.steffe_begin ... w.steffe_end) into the file
meshGFaceOptimize.cpp:


int edgeSwapPass(GFace *gf, std::set<MTri3*, compareTri3Ptr> &allTris,
                 const swapCriterion &cr,bidimMeshData & data)
{

//w.steffe_begin

//Here I put into a set the mesh vertices which are also face vertices:
      std::list<GEdge*> edges = gf->edges();
      std::set<MVertex*, MVertexLessThanNum> all_geom_vertices;
      std::list<GEdge*>::iterator ite = edges.begin();
      while(ite != edges.end()){
        if(!(*ite)->isMeshDegenerated()){
          MVertex *v1 = (*ite)->lines[0]->getVertex(0);
          int i=(*ite)->lines.size()-1;
          MVertex *v2 = (*ite)->lines[i]->getVertex(1);
          all_geom_vertices.insert(v1);
          all_geom_vertices.insert(v2);
        }
        ++ite;
      }

//Here I put into a set all the triangle vertices (before edge swap):

      std::set<MVertex*, MVertexLessThanNum> all_tria_vertices;
      int NTbefore=0;
      for(std::set<MTri3*,compareTri3Ptr>::iterator ti =
allTris.begin(); ti != allTris.end(); ++ti)
      if(!(*ti)->isDeleted()) {
	NTbefore++;
        MTriangle *t = (*ti)->tri();
	for(int j = 0; j < 3; ++j){
	  MVertex *v = t->getVertex(j);
          all_tria_vertices.insert(v);
	}
      }
      int NVbefore= all_tria_vertices.size();
      std::set<MVertex*, MVertexLessThanNum>::iterator
itv=all_geom_vertices.begin();
      while(itv != all_geom_vertices.end()){
        MVertex *v = *itv;
//   face vertices must belong to the set of all_tria_vertices
//   The test is OK at this point
        assert(all_tria_vertices.find(v)!=all_tria_vertices.end());
       ++itv;
      }
//w.steffe_end

  typedef std::set<MTri3*, compareTri3Ptr> CONTAINER;

  int nbSwapTot = 0;
  std::set<swapquad> configs;
  for(int iter = 0; iter < 1200; iter++){
    int nbSwap = 0;
    std::vector<MTri3*> newTris;

//w.steffe_begin
//  deleted_tria_vertices is a global variable handled inside of
//  the function edgeSwap()
    deleted_tria_vertices.clear();
//w.steffe_end

    for(CONTAINER::iterator it = allTris.begin(); it != allTris.end();
++it){
      if(!(*it)->isDeleted()){
        for(int i = 0; i < 3; i++){
          if(edgeSwap(configs, *it, gf, i, newTris, cr, data)){
            nbSwap++;
            break;
          }
        }
      }
      else{
        delete *it;
        CONTAINER::iterator itb = it;
        ++it;
        allTris.erase(itb);
        if(it == allTris.end()) break;
      }
    }
    allTris.insert(newTris.begin(), newTris.end());
    nbSwapTot += nbSwap;
//w.steffe_begin
// new_tria_vertices is the set of vertices belonging to newTris:
      std::set<MVertex*, MVertexLessThanNum> new_tria_vertices;
      for(int i=0; i<newTris.size(); i++){
        MTriangle *t = newTris[i]->tri();
	for(int j = 0; j < 3; ++j){
	  MVertex *v = t->getVertex(j);
          new_tria_vertices.insert(v);
	}
      }
      int NTafter=0;
      all_tria_vertices.clear();
      for(std::set<MTri3*,compareTri3Ptr>::iterator ti =
allTris.begin(); ti != allTris.end(); ++ti)
      if(!(*ti)->isDeleted()) {
	NTafter++;
        MTriangle *t = (*ti)->tri();
	for(int j = 0; j < 3; ++j){
	  MVertex *v = t->getVertex(j);
          all_tria_vertices.insert(v);
	}
      }
      int NVafter= all_tria_vertices.size();

//   face vertices must belong to the set of all_tria_vertices
//   The test is OK at this point
      itv=all_geom_vertices.begin();
      while(itv != all_geom_vertices.end()){
        MVertex *v = *itv;
//    TEST FAILS AT THIS POINT:
        assert(all_tria_vertices.find(v)!=all_tria_vertices.end());
       ++itv;

//    Using gdb I have seen:
//    COMPARISON BETWEEB NVafter and NVbefore shows that one vertex was
//    lost
//    COMPARISON BETWEEB NTafter and NTbefore shows that triangle number
//    decreased by 1 which should not happen in Edge Swap
//
//    all deleted_tria_vertices (registered inside of the function
//    edgeSwap) were included in new_tria_vertices.
//
//    I DO NOT UNDERSTAND WHERE THE VERTEX/TRIANGLE LOSS MAY
//    HAPPEN !!!!!!!!!

      }
//w.steffe_end
    if(nbSwap == 0) break;
  }


  return nbSwapTot;
}


ANY SUGGESTION ?

Walter