Merge pull request #106 from laetitia-m/feature/ls-cuttet

Level-set discretization function `PMMG_ls-cuttet`. Note that multi-material is not supported yet by this function.

CMakeLists.txt

@@ -252,7 +252,7 @@ IF ( DOWNLOAD_MMG )
 
   EXTERNALPROJECT_ADD ( Mmg
     GIT_REPOSITORY https://github.com/MmgTools/mmg.git
-    GIT_TAG 3a5e3abc38f7f100a481bcebe8aa23e105d7c381
+    GIT_TAG 5788cf5b872d6f2e627b4279260ac39706d2a7f6
     INSTALL_COMMAND ${CMAKE_MAKE_PROGRAM} install
     CMAKE_ARGS ${MMG_ARGS} -DUSE_ELAS=OFF ${COMPILER_CFG} ${FLAGS_CFG}
     ${SCOTCH_CFG} ${VTK_CFG} -DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE}

src/analys_pmmg.c

@@ -1602,7 +1602,7 @@ int PMMG_loopr(PMMG_pParMesh parmesh,PMMG_hn_loopvar *var ) {
   doublevalues = parmesh->int_node_comm->doublevalues;
 
   /* Loop on near-parallel edges */
-  for( k = 1; k <= var->hash->max; k++ ) {
+  for( k = 0; k <= var->hash->max; k++ ) {
     ph = &var->hash->geom[k];
     if( !ph->a ) continue;
 

src/communicators_pmmg.c

@@ -871,24 +871,7 @@ void PMMG_tria2elmFace_coords( PMMG_pParMesh parmesh ) {
     ifac   = ptt->cc%4;
 
     /* Get triangle node with highest coordinates */
-    iploc  = 0;
-    ppt = &mesh->point[ptt->v[0]];
-    cmax[0] = ppt->c[0];
-    cmax[1] = ppt->c[1];
-    cmax[2] = ppt->c[2];
-    for( iloc=1; iloc<3; iloc++ ) {
-      ppt = &mesh->point[ptt->v[iloc]];
-      for( idim=0; idim<3; idim++ ) {
-        if( ppt->c[idim] - cmax[idim] < -MMG5_EPSOK*20 ) break;
-        if( ppt->c[idim] - cmax[idim] >  MMG5_EPSOK*20 ) {
-          cmax[0] = ppt->c[0];
-          cmax[1] = ppt->c[1];
-          cmax[2] = ppt->c[2];
-          iploc  = iloc;
-          break;
-        }
-      }
-    }
+    iploc = PMMG_tria_highestcoord(mesh,ptt->v);
 
     /* Store ie-ifac-iploc in index1 */
     grp->face2int_face_comm_index1[i] = 12*ie+3*ifac+iploc;
@@ -900,6 +883,42 @@ void PMMG_tria2elmFace_coords( PMMG_pParMesh parmesh ) {
   }
 }
 
+/**
+ * \param mesh  pointer toward the mesh structure
+ * \param ptt_v indices of a triangle vertices
+ * \return iploc (0, 1 or 2) local node index
+ *
+ * Get triangle node with highest coordinates
+ *
+ */
+int PMMG_tria_highestcoord( MMG5_pMesh mesh, MMG5_int *ptt_v) {
+  MMG5_pPoint  ppt;
+  int          idim,iloc,iploc;
+  double       cmax[3];
+
+  /* Get triangle node with highest coordinates */
+  iploc  = 0;
+  ppt = &mesh->point[ptt_v[0]];
+  cmax[0] = ppt->c[0];
+  cmax[1] = ppt->c[1];
+  cmax[2] = ppt->c[2];
+  for( iloc=1; iloc<3; iloc++ ) {
+    ppt = &mesh->point[ptt_v[iloc]];
+    for( idim=0; idim<3; idim++ ) {
+      if( ppt->c[idim] - cmax[idim] < -MMG5_EPSOK*20 ) break;
+      if( ppt->c[idim] - cmax[idim] >  MMG5_EPSOK*20 ) {
+        cmax[0] = ppt->c[0];
+        cmax[1] = ppt->c[1];
+        cmax[2] = ppt->c[2];
+        iploc  = iloc;
+        break;
+      }
+    }
+  }
+
+  return iploc;
+}
+
 /**
  * \param parmesh pointer toward a parmesh structure
  * \return 0 if fail, 1 if success

src/hash_pmmg.c

@@ -28,6 +28,7 @@
  * \author Algiane Froehly (Inria/UBordeaux)
  * \author Nikolas Pattakos (Inria)
  * \author Luca Cirrottola (Inria)
+ * \author Laetitia Mottet (UBordeaux)
  * \version
  * \copyright
  *
@@ -273,3 +274,87 @@ int PMMG_bdryUpdate( MMG5_pMesh mesh )
 
   return PMMG_SUCCESS;
 }
+
+/**
+ * \param mesh pointer toward the mesh structure.
+ * \param hash pointer toward the hash table of edges.
+ * \param a index of the first extremity of the edge.
+ * \param b index of the second extremity of the edge.
+ * \param k index of new point along the edge [a,b].
+ * \param s If ls mode in ParMmg: index of new point in internal edge communicator;
+ *          otherwise, the value stored in variable s.
+ * \return PMMG_SUCCESS if success, PMMG_FAILURE if fail (edge is not found).
+ *
+ * Update the index of the new point stored along the edge \f$[a;b]\f$ (similar to MMG5_hashUpdate in mmg).
+ * If ls mode in ParMmg: update the index of the new point in internal edge communicator in variable s;
+ * otherwise, update the value stored in variable s.
+ *
+ */
+int PMMG_hashUpdate_all(MMG5_Hash *hash, MMG5_int a,MMG5_int b,MMG5_int k,MMG5_int s) {
+  MMG5_hedge  *ph;
+  MMG5_int     key;
+  MMG5_int    ia,ib;
+
+  ia  = MG_MIN(a,b);
+  ib  = MG_MAX(a,b);
+  key = (MMG5_KA*(int64_t)ia + MMG5_KB*(int64_t)ib) % hash->siz;
+  ph  = &hash->item[key];
+
+  while ( ph->a ) {
+    if ( ph->a == ia && ph->b == ib ) {
+      ph->k = k;
+      ph->s = s;
+      return PMMG_SUCCESS;
+    }
+
+    if ( !ph->nxt ) return PMMG_FAILURE;
+
+    ph = &hash->item[ph->nxt];
+
+  }
+
+  return PMMG_FAILURE;
+}
+
+/**
+ * \param hash pointer toward the hash table of edges.
+ * \param a index of the first extremity of the edge.
+ * \param b index of the second extremity of the edge.
+ * \param k index of new point along the edge [a,b].
+ * \param s If ls mode in ParMmg: index of new point in internal edge communicator;
+ *          otherwise, the value stored in variable s.
+ * \return PMMG_SUCCESS if success, PMMG_FAILURE if fail (edge is not found).
+ *
+ * Find the index of the new point stored along the edge \f$[a;b]\f$ (similar to MMG5_hashGet in mmg).
+ * If ls mode in ParMmg: find the index of the new point in internal edge communicator;
+ * otherwise, find the value stored in variable s.
+ *
+ */
+MMG5_int PMMG_hashGet_all(MMG5_Hash *hash,MMG5_int a,MMG5_int b,MMG5_int *k,MMG5_int *s) {
+  MMG5_hedge  *ph;
+  MMG5_int    key;
+  MMG5_int    ia,ib;
+
+  if ( !hash->item ) return 0;
+
+  ia  = MG_MIN(a,b);
+  ib  = MG_MAX(a,b);
+  key = (MMG5_KA*(int64_t)ia + MMG5_KB*(int64_t)ib) % hash->siz;
+  ph  = &hash->item[key];
+
+  if ( !ph->a )  return PMMG_FAILURE;
+  if ( ph->a == ia && ph->b == ib )  {
+    *k = ph->k;
+    *s = ph->s;
+    return PMMG_SUCCESS;
+  }
+  while ( ph->nxt ) {
+    ph = &hash->item[ph->nxt];
+    if ( ph->a == ia && ph->b == ib ) {
+      *k = ph->k;
+      *s = ph->s;
+      return PMMG_SUCCESS;
+    }
+  }
+  return PMMG_FAILURE;
+}

src/libparmmg.c

@@ -324,7 +324,7 @@ int PMMG_preprocessMesh_distributed( PMMG_pParMesh parmesh )
     mesh->info.sethmax = 1;
   }
 
-  /* Note :: Needed before ls discretization to include tetras with poor qualities
+  /* Note: Needed before ls discretization to include tetras with poor qualities
       inside one or the other part of the level-set */
   if ( !MMG3D_tetraQual( mesh, met, 0 ) ) {
     return PMMG_STRONGFAILURE;
@@ -356,7 +356,7 @@ int PMMG_preprocessMesh_distributed( PMMG_pParMesh parmesh )
 
       /* Convert tria index into iel face index (it needs a valid cc field in
        * each tria), and tag xtetra face as PARBDY before the tag is transmitted
-       * to edges and nodes - Why we do that? */
+       * to edges and nodes */
       if ( parmesh->myrank < parmesh->info.npartin ) {
         PMMG_tria2elmFace_coords( parmesh );
       }

src/libparmmg1.c

@@ -764,7 +764,7 @@ int PMMG_parmmglib1( PMMG_pParMesh parmesh )
         }
 
         if ( parmesh->iter < parmesh->niter-1 && (!parmesh->info.inputMet) ) {
-          /* Delete the metrec computed by Mmg except at last iter */
+          /* Delete the metric computed by Mmg except at last iter */
           PMMG_DEL_MEM(mesh,met->m,double,"internal metric");
         }