Correction of a bug: hashTria was bad at eliminating duplicated faces in NOM situation; added function remdup

Author: dapogny

src/common/hash.c

@@ -44,7 +44,7 @@
  *
  * Create surface adjacency
  *
- * \remark the ph->s field computation is useless in mmgs.
+ * \remark the ph->s field computation (containing the number of triangles shared by a hashed edge)  is useless in mmgs.
  *
  *
  * \remark: as all triangles are mesh boundaries, we do not need to mark their

src/mmg3d/API_functions_3d.c

@@ -2153,6 +2153,9 @@ int MMG3D_Set_iparameter(MMG5_pMesh mesh, MMG5_pSol sol, int iparam,MMG5_int val
       mesh->info.dhd    = MMG5_ANGEDG;
     }
     break;
+  case MMG3D_IPARAM_kiso :
+    mesh->info.kiso = val;
+    break;
   case MMG3D_IPARAM_nofem :
     mesh->info.setfem = (val==1)? 0 : 1;
     break;

src/mmg3d/analys_3d.c

@@ -1368,6 +1368,14 @@ int MMG3D_analys(MMG5_pMesh mesh) {
 
   /* Set surface triangles to required in nosurf mode or for parallel boundaries */
   MMG3D_set_reqBoundaries(mesh);
+  
+  /* Remove duplicate triangles */
+  memset ( &hash, 0x0, sizeof(MMG5_Hash));
+  if ( !MMG5_remdup(mesh,&hash) ) {
+    MMG5_DEL_MEM(mesh,hash.item);
+    fprintf(stderr,"\n  ## Problem in removing duplicate faces. Exit program.\n");
+    return 0;
+  }
 
   /* create surface adjacency */
   memset ( &hash, 0x0, sizeof(MMG5_Hash));
@@ -1449,6 +1457,11 @@ int MMG3D_analys(MMG5_pMesh mesh) {
     MMG5_DEL_MEM(mesh,mesh->xpoint);
     return 0;
   }
+  
+  /* printf("On est bien ici maintenant??\n");
+  MMG3D_packMesh(mesh,NULL,NULL);
+  MMG3D_saveMesh(mesh,mesh->nameout);
+  exit(0);*/
 
   /* check subdomains connected by a vertex and mark these vertex as corner and
      required */

src/mmg3d/hash_3d.c

@@ -504,11 +504,10 @@ int MMG5_setEdgeNmTag(MMG5_pMesh mesh, MMG5_Hash *hash) {
         /* Set edge tag and point tags to MG_REQ if the non-manifold edge shared
          * separated domains */
         if ( ph->s > 3 ) {
-          start = ptt->cc/4;
+          start = ptt->cc/4; // ptt->cc contains information about a supporting tetrahedron
           assert(start);
           pt = &mesh->tetra[start];
 
-
           for (ia=0; ia<6; ++ia) {
             ipa = MMG5_iare[ia][0];
             ipb = MMG5_iare[ia][1];
@@ -517,7 +516,6 @@ int MMG5_setEdgeNmTag(MMG5_pMesh mesh, MMG5_Hash *hash) {
           }
           assert(ia<6);
 
-
           /* Travel throug the shell of the edge until reaching a tetra without adjacent
            * or until reaching the starting tetra */
           iface = ptt->cc%4;
@@ -802,6 +800,102 @@ int MMG5_setVertexNmTag(MMG5_pMesh mesh,uint16_t func(uint16_t) ) {
   return 1;
 }
 
+/**
+ * \param mesh pointer towar the mesh structure.
+ * \param hash  hash table.
+ * \return 1 if success, 0 if failed.
+ *
+ * Remove duplicate faces from the  triangle mesh
+ *
+ */
+int MMG5_remdup(MMG5_pMesh mesh,MMG5_Hash *hash) {
+  MMG5_pTria     pt;
+  MMG5_hedge     *ph;
+  MMG5_int       dup,k,kk,hmax,ip0,ip1,ip2,min,max,sum,key;
+
+  dup = 0;
+  
+  /* Allocate hash table params */
+  hash->siz  = mesh->np;
+  hmax       = MG_MAX(hash->siz,mesh->nt);
+  hash->max  = hmax + 1;
+  hash->nxt  = hash->siz;
+  MMG5_ADD_MEM(mesh,(hash->max+1)*sizeof(MMG5_hedge),"hash table",return 0);
+  MMG5_SAFE_CALLOC(hash->item,hash->max+1,MMG5_hedge,return 0);
+  
+  /* Hash triangles */
+  for(k=1; k<=mesh->nt; k++) {
+    pt = &mesh->tria[k];
+    if ( !MG_EOK(pt) )  continue;
+    
+    ip0  = pt->v[0];
+    ip1  = pt->v[1];
+    ip2  = pt->v[2];
+    
+    min = MG_MIN(ip0,MG_MIN(ip1,ip2));
+    max = MG_MAX(ip0,MG_MAX(ip1,ip2));
+    sum = ip0 + ip1 + ip2;
+    
+    key = (MMG5_KA*(int64_t)min + MMG5_KB*(int64_t)max + MMG5_KC*(int64_t)sum) % hash->siz;
+    ph  = &hash->item[key];
+    
+    /* Empty cell: store tria */
+    if ( ph->a == 0 ) {
+      ph->a   = min;
+      ph->b   = max;
+      ph->k   = sum;
+      ph->nxt = 0;
+      continue;
+    }
+
+    while ( ph->a ) {
+      /* Remove duplicate */
+      if ( ph->a == min && ph->b == max && ph->k == sum ) {
+        pt->v[0] = 0;
+        dup++;
+        break;
+      }
+      else if ( !ph->nxt ) {
+        ph->nxt = hash->nxt;
+        hash->nxt++;
+        ph = &hash->item[ph->nxt];
+        assert(ph);
+
+        if ( hash->nxt >= hash->max-1 ) {
+          if ( mesh->info.ddebug ) {
+            fprintf(stderr,"\n  ## Warning: %s: memory alloc problem (edge):"
+                    " %" MMG5_PRId "\n",__func__,hash->max);
+          }
+          MMG5_TAB_RECALLOC(mesh,hash->item,hash->max,MMG5_GAP,MMG5_hedge,
+                              "MMG5_edge",
+                              MMG5_DEL_MEM(mesh,hash->item);
+                              return 0);
+
+          ph = &hash->item[hash->nxt];
+        }
+        
+        ph->a     = min;
+        ph->b     = max;
+        ph->k     = sum;
+        ph->nxt   = 0;
+        break;
+      }
+      else {
+        ph = &hash->item[ph->nxt];
+      }
+    }
+  }
+    
+  /* Display information */
+  if ( abs(mesh->info.imprim) > 5 && dup > 0 ) {
+    fprintf(stdout,"  ## ");  fflush(stdout);
+    if ( dup > 0 )  fprintf(stdout," %" MMG5_PRId " duplicate removed",dup);
+    fprintf(stdout,"\n");
+  }
+    
+  return 1;
+}
+
 /**
  * \param mesh pointer towar the mesh structure.
  * \param hash edges hash table.

src/mmg3d/inout_3d.c

@@ -140,7 +140,6 @@ int MMG3D_loadMesh_opened(MMG5_pMesh mesh,FILE *inm,int bin) {
   mesh->np = mesh->nt = mesh->ne = 0;
   nref = 0;
 
-
   if (!bin) {
     strcpy(chaine,"D");
     while(fscanf(inm,"%127s",&chaine[0])!=EOF && strncmp(chaine,"End",strlen("End")) ) {
@@ -523,7 +522,6 @@ int MMG3D_loadMesh_opened(MMG5_pMesh mesh,FILE *inm,int bin) {
     }
   }
 
-
   /* get corners */
   if(ncor) {
     rewind(inm);
@@ -1056,6 +1054,7 @@ int MMG3D_loadMesh(MMG5_pMesh mesh,const char *filename) {
   if( ier < 1 ) return ier;
 
   fclose(inm);
+  
   return 1;
 }
 

src/mmg3d/libmmg3d.c

@@ -101,7 +101,7 @@ MMG5_int MMG3D_bdryBuild(MMG5_pMesh mesh) {
   MMG5_hgeom  *ph;
   MMG5_int    k,nr;
   int         i;
-
+  
   /* rebuild triangles*/
   if ( mesh->tria )
     MMG5_DEL_MEM(mesh,mesh->tria);
@@ -1296,9 +1296,11 @@ int MMG3D_mmg3dls(MMG5_pMesh mesh,MMG5_pSol sol,MMG5_pSol umet) {
   }
 
   /* clean old isosurface */
-  if ( !MMG3D_Clean_isoSurf(mesh) ) {
-    fprintf(stderr,"\n  ## Unable to clean old isosurface.\n");
-    _LIBMMG5_RETURN(mesh,met,sol,MMG5_STRONGFAILURE);
+  if ( ! mesh->info.kiso ) {
+    if ( !MMG3D_Clean_isoSurf(mesh) ) {
+      fprintf(stderr,"\n  ## Unable to clean old isosurface.\n");
+      _LIBMMG5_RETURN(mesh,met,sol,MMG5_STRONGFAILURE);
+    }
   }
 
   chrono(OFF,&(ctim[1]));

src/mmg3d/libmmg3d.h

@@ -167,6 +167,7 @@ enum MMG3D_Param {
   MMG3D_IPARAM_anisosize,                 /*!< [1/0], Turn on/off anisotropic metric creation when no metric is provided */
   MMG3D_IPARAM_octree,                    /*!< [n], Max number of vertices per PROctree cell (DELAUNAY) */
   MMG3D_IPARAM_nosizreq,                  /*!< [0/1], Allow/avoid overwriting of sizes at required vertices (advanced usage) */
+  MMG3D_IPARAM_kiso,                      /*!< [1/0], Keep previous isoline in mesh (in ls mode) */
   MMG3D_IPARAM_isoref,                    /*!< [0/n], Isosurface boundary material reference */
   MMG3D_DPARAM_angleDetection,            /*!< [val], Value for angle detection (degrees) */
   MMG3D_DPARAM_hmin,                      /*!< [val], Minimal edge length */

src/mmg3d/libmmg3d_private.h

@@ -278,6 +278,7 @@ int16_t MMG5_openCoquilTravel(MMG5_pMesh,MMG5_int,MMG5_int,MMG5_int*,MMG5_int*,i
 int  MMG3D_get_shellEdgeTag(MMG5_pMesh,MMG5_int,int8_t,uint16_t*,MMG5_int *);
 int  MMG5_settag(MMG5_pMesh,MMG5_int,int,uint16_t,int);
 int  MMG5_deltag(MMG5_pMesh,MMG5_int,int,uint16_t);
+int  MMG5_remdup(MMG5_pMesh mesh, MMG5_Hash *hash);
 int  MMG5_setNmTag(MMG5_pMesh mesh, MMG5_Hash *hash);
 int  MMG5_setVertexNmTag(MMG5_pMesh mesh,uint16_t func(uint16_t) );
 int  MMG5_chkcol_int(MMG5_pMesh,MMG5_pSol,MMG5_int,int8_t,int8_t,int64_t*,int,int8_t);

src/mmg3d/libmmg3d_tools.c

@@ -382,6 +382,12 @@ int MMG3D_storeknownar(int argc,char *argv[],MMG5_pMesh mesh,MMG5_pSol met,
           MMG_ARGV_APPEND(argv, mmgArgv, i, *mmgArgc,return 0);
         }
         break;
+      case 'k':
+        if ( !strcmp(argv[i],"-kiso") ) {
+          if ( !MMG3D_Set_iparameter(mesh,met,MMG3D_IPARAM_kiso,1) )
+            return 0;
+        }
+        break;
       case 'l':
         if ( !strcmp(argv[i],"-lag") ) {
           if ( ++i < argc && isdigit(argv[i][0]) ) {
@@ -719,7 +725,7 @@ int MMG3D_parsar(int argc,char *argv[],MMG5_pMesh mesh,MMG5_pSol met,MMG5_pSol s
   MMG5_SAFE_MALLOC( mmgArgv, argc, char*,return 0);
   MMG_ARGV_APPEND ( argv, mmgArgv, 0, mmgArgc,return 0);
 
-  int ier = MMG3D_storeknownar( argc,argv,mesh,met,sol,&mmgArgc,mmgArgv);
+  int ier = MMG3D_storeknownar(argc,argv,mesh,met,sol,&mmgArgc,mmgArgv);
 
   /* Third step: treat unknown args */
   if ( ier ) {