Merge branch 'develop' of https://github.com/MmgTools/ParMmg into dev…

…elop
laetitia-m · Sep 4, 2023 · 6dccbff · 6dccbff
2 parents df44b15 + f981ff8
commit 6dccbff
Show file tree

Hide file tree

Showing 11 changed files with 391 additions and 119 deletions.
diff --git a/.github/workflows/long-tests.yml b/.github/workflows/long-tests.yml
@@ -5,7 +5,6 @@ on:
   push:
   # run tests on PR events
   pull_request:
-    types: [synchronize]
 
   # run tests manually on a given branch (default is master)
   workflow_dispatch:

diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -254,7 +254,7 @@ IF ( DOWNLOAD_MMG )
 
   EXTERNALPROJECT_ADD ( Mmg
     GIT_REPOSITORY https://github.com/MmgTools/mmg.git
-    GIT_TAG 6de02d0f
+    GIT_TAG 31769770b
     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}

diff --git a/cmake/testing/pmmg_tests.cmake b/cmake/testing/pmmg_tests.cmake
@@ -229,9 +229,8 @@ IF( BUILD_TESTING )
 
     set_property(TEST PvtuOut-RenameOut-2
       PROPERTY PASS_REGULAR_EXPRESSION
-      "${OutputVtkRenameFilename}.*${OutputVtkRenameWarning}.*${OutputVtkErr};
-       ${OutputVtkRenameWarning}.*${OutputVtkErr}.*${OutputVtkRenameFilename};
-       ${OutputVtkRenameFilename}.*${OutputVtkErr}.*${OutputVtkRenameWarning}")
+      "${OutputVtkRenameFilename}.*${OutputVtkRenameWarning};
+       ${OutputVtkRenameWarning}.*${OutputVtkRenameFilename}")
 
     ###############################################################################
     #####
@@ -419,8 +418,10 @@ IF( BUILD_TESTING )
       -sol ${CI_DIR}/LevelSet/centralized/3D-cube-ls.sol
       -out ${CI_DIR_RESULTS}/3D-cube-ls-CenIn-DisOut-${NP}-out.pvtu)
 
-    set_property(TEST ls-CenIn-DisOut-${NP}
-      PROPERTY PASS_REGULAR_EXPRESSION "${OutputVtkErr}")
+    IF ( (NOT VTK_FOUND) OR USE_VTK MATCHES OFF )
+      set_property(TEST ls-CenIn-DisOut-${NP}
+        PROPERTY PASS_REGULAR_EXPRESSION "${OutputVtkErr}")
+    ENDIF ( )
 
   endforeach()
 
@@ -524,7 +525,7 @@ IF( BUILD_TESTING )
   # and  to write distributed output fields in VTK   format
   add_test( NAME fields-DisIn-DisOutVTK-2
     COMMAND ${MPIEXEC} ${MPI_ARGS} ${MPIEXEC_NUMPROC_FLAG} 2 $<TARGET_FILE:${PROJECT_NAME}>
-    ${CI_DIR}/LevelSet/distributed/3D-cube.mesh
+    ${CI_DIR}/LevelSet/distributed/3D-cube.mesh -v 10
     -field ${CI_DIR}/LevelSet/distributed/3D-cube-fields.sol
     -out ${CI_DIR_RESULTS}/3D-cube-fields-DisIn-DisOutVTK-2-out.pvtu)
 
@@ -533,9 +534,8 @@ IF( BUILD_TESTING )
 
   set_property(TEST fields-DisIn-DisOutVTK-2
     PROPERTY PASS_REGULAR_EXPRESSION
-    "${InputDistributedFields}.*${OutputVtkFields}.*${OutputVtkErr};
-     ${OutputVtkFields}.*${OutputVtkErr}.*${InputDistributedFields};
-     ${InputDistributedFields}.*${OutputVtkErr}.*${OutputVtkFields}")
+    "${InputDistributedFields}.*${OutputVtkFields};
+     ${OutputVtkFields}.*${InputDistributedFields}")
 
   # Test to write distributed output fields and metric in Medit format
   add_test( NAME fields-DisIn-DisOutMesh-2

diff --git a/src/API_functionsf_pmmg.c b/src/API_functionsf_pmmg.c
@@ -1062,6 +1062,16 @@ FORTRAN_NAME(PMMG_PARMMGLIB_DISTRIBUTED,pmmg_parmmglib_distributed,
   return;
 }
 
+/**
+ * See \ref PMMG_parmmgls_distributed function in \ref libparmmg.h file.
+ */
+FORTRAN_NAME(PMMG_PARMMGLS_DISTRIBUTED,pmmg_parmmgls_distributed,
+             (PMMG_pParMesh *parmesh,int* retval),
+             (parmesh,retval)) {
+  *retval = PMMG_parmmgls_distributed(*parmesh);
+  return;
+}
+
 /**
  * See \ref PMMG_parmmglib_centralized function in \ref libparmmg.h file.
  */

diff --git a/src/analys_pmmg.c b/src/analys_pmmg.c
@@ -1510,6 +1510,11 @@ int PMMG_update_nmgeom(PMMG_pParMesh parmesh,MMG5_pMesh mesh){
   return 1;
 }
 
+static inline
+int MMG5_skip_nonOldParBdy ( int8_t tag ) {
+  return !(tag & MG_OLDPARBDY);
+}
+
 /**
  * \param mesh pointer toward the mesh  structure.
  * \return 0 if point is singular, 1 otherwise.
@@ -1520,88 +1525,9 @@ int PMMG_update_nmgeom(PMMG_pParMesh parmesh,MMG5_pMesh mesh){
  */
 static inline
 int PMMG_update_singul(PMMG_pParMesh parmesh,MMG5_pMesh mesh) {
-  MMG5_pTetra         ptet;
-  MMG5_pPoint         ppt;
-  MMG5_Hash           hash;
-  int                 k,i;
-  int                 nc, nre, ng, nrp, nm, ier;
-
-  /* Second: seek the non-required non-manifold points and try to analyse
-   * whether they are corner or required. */
-
-  /* Hash table used by boulernm to store the special edges passing through
-   * a given point */
-  if ( ! MMG5_hashNew(mesh,&hash,mesh->np,(int)(3.71*mesh->np)) ) return 0;
-
-  nc = nre = 0;
-  ++mesh->base;
-  for (k=1; k<=mesh->ne; ++k) {
-    ptet = &mesh->tetra[k];
-    if ( !MG_EOK(ptet) ) continue;
-
-    for ( i=0; i<4; ++i ) {
-      ppt = &mesh->point[ptet->v[i]];
-
-      /* Skip non-previously-parallel points */
-      if ( !(ppt->tag & MG_OLDPARBDY) ) continue;
-
-      if ( (!MG_VOK(ppt)) || (ppt->flag==mesh->base)  ) continue;
-      ppt->flag = mesh->base;
-
-      if ( (!(ppt->tag & MG_NOM)) || (ppt->tag & MG_REQ) ) continue;
 
-      ier = MMG5_boulernm(mesh,&hash, k, i, &ng, &nrp, &nm);
-      if ( ier < 0 ) return 0;
-      else if ( !ier ) continue;
+  return MMG5_setVertexNmTag(mesh,MMG5_skip_nonOldParBdy);
 
-      if ( (ng+nrp+nm) > 2 ) {
-        /* More than 2 feature edges are passing through the point: point is
-         * marked as corner */
-        ppt->tag |= MG_CRN + MG_REQ;
-        ppt->tag &= ~MG_NOSURF;
-        nre++;
-        nc++;
-      }
-      else if ( (ng == 2) || (nrp == 2) || (nm == 2) ) {
-        /* Exactly 2 edges of same type are passing through the point: do
-         * nothing */
-        continue;
-      }
-      else if ( (ng+nrp+nm) == 2 ) {
-        /* 2 edges of different type are passing through the point: point is
-         * marked as required */
-        ppt->tag |= MG_REQ;
-        ppt->tag &= ~MG_NOSURF;
-        nre++;
-      }
-      else if ( ng == 1 && !nrp ){
-        ppt->tag |= MG_CRN + MG_REQ;
-        ppt->tag &= ~MG_NOSURF;
-        nre++;
-        nc++;
-      }
-      else if ( (ng+nrp+nm) == 1 ){
-        /* Only 1 feature edge is passing through the point: point is
-         * marked as corner */
-        assert ( (ng == 1) || (nrp==1) || (nm==1) );
-        ppt->tag |= MG_CRN + MG_REQ;
-        ppt->tag &= ~MG_NOSURF;
-        nre++;
-        nc++;
-      }
-      else {
-        assert ( 0 && "unexpected case");
-      }
-    }
-  }
-
-  /* Free the edge hash table */
-  MMG5_DEL_MEM(mesh,hash.item);
-
-  if ( mesh->info.ddebug || abs(mesh->info.imprim) > 3 )
-    fprintf(stdout,"     %d corner and %d required vertices added\n",nc,nre);
-
-  return 1;
 }
 
 /**
@@ -2399,9 +2325,10 @@ int PMMG_setdhd(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MMG5_HGeom *pHash,MPI_Comm
  * Check all boundary triangles.
  */
 int PMMG_analys_tria(PMMG_pParMesh parmesh,MMG5_pMesh mesh) {
+  int       ier;
 
   /**--- stage 1: data structures for surface */
-  if ( abs(mesh->info.imprim) > 3 )
+  if ( parmesh->info.imprim > PMMG_VERB_VERSION )
     fprintf(stdout,"\n  ** SURFACE ANALYSIS\n");
 
   /* create tetra adjacency */
@@ -2410,6 +2337,16 @@ int PMMG_analys_tria(PMMG_pParMesh parmesh,MMG5_pMesh mesh) {
     return 0;
   }
 
+  /* Update the xtetra data after the ls discretization */
+  if ( mesh->info.iso && mesh->info.opnbdy ) {
+    ier = MMG3D_update_xtetra ( mesh );
+    if ( !ier ) {
+      fprintf(stderr,"\n  ## Problem when updating the xtetra data after ls discretization."
+              " Exit program.\n");
+      return 0;
+    }
+  }
+
   /* create prism adjacency */
   if ( !MMG3D_hashPrism(mesh) ) {
     fprintf(stderr,"\n  ## Prism hashing problem. Exit program.\n");
@@ -2665,7 +2602,7 @@ int PMMG_analys(PMMG_pParMesh parmesh,MMG5_pMesh mesh,MPI_Comm comm) {
   }
 
   /**--- stage 2: surface analysis */
-  if ( abs(mesh->info.imprim) > 5  || mesh->info.ddebug )
+  if ( parmesh->info.imprim > PMMG_VERB_VERSION )
     fprintf(stdout,"  ** SETTING TOPOLOGY\n");
 
   /* identify connexity */

diff --git a/src/libparmmg.c b/src/libparmmg.c
@@ -70,9 +70,6 @@ int PMMG_check_inputData(PMMG_pParMesh parmesh)
       fprintf(stderr,
               "  ## Error: lagrangian mode unavailable (MMG3D_IPARAM_lag):\n");
       return 0;
-    } else if ( mesh->info.iso ) {
-      fprintf(stderr,"\n\n  ## WARNING: level-set discretisation under construction. \n\n");
-      // return 0;
     } else if ( mesh->info.optimLES && met->size==6 ) {
       fprintf(stdout,"  ## Error: strong mesh optimization for LES methods"
               " unavailable (MMG3D_IPARAM_optimLES) with an anisotropic metric.\n");
@@ -240,25 +237,28 @@ int PMMG_preprocessMesh( PMMG_pParMesh parmesh )
 int PMMG_preprocessMesh_distributed( PMMG_pParMesh parmesh )
 {
   MMG5_pMesh mesh;
-  MMG5_pSol  met;
+  MMG5_pSol  met,ls;
+  int8_t     tim;
+  char       stim[32];
+  mytime     ctim[TIMEMAX];
   int ier = PMMG_SUCCESS;
 
   mesh = parmesh->listgrp[0].mesh;
   met  = parmesh->listgrp[0].met;
+  ls   = parmesh->listgrp[0].ls;
 
   assert ( ( mesh != NULL ) && ( met != NULL ) && "Preprocessing empty args");
 
-  if (mesh->info.iso) {
-// Just print a warning saying that the feature is not implemented and
-// deallocate the isovalue structure (as the ls is not interpolated during the
-// remeshing step, the continuous integration tests will fail otherwise)
-    if ( parmesh->myrank == parmesh->info.root) {
-      fprintf(stdout,"Isovalue discretization is not yet implemented:"
-        " Deallocation of the level-set structure.\n");
-    }
-    PMMG_DEL_MEM(mesh,parmesh->listgrp[0].ls->m,double,"ls structure");
-    parmesh->listgrp[0].ls->np = 0;
-  }
+  // if (mesh->info.iso) {
+  //   // Just print a warning saying that the feature is not implemented and
+  //   // deallocate the isovalue structure (as the ls is not interpolated during the
+  //   // remeshing step, the continuous integration tests will fail otherwise)
+  //   if ( parmesh->myrank == parmesh->info.root) {
+  //     fprintf(stdout,"Isovalue discretization is under development.\n");
+  //   }
+  //   PMMG_DEL_MEM(mesh,parmesh->listgrp[0].ls->m,double,"ls structure");
+  //   parmesh->listgrp[0].ls->np = 0;
+  // }
 
 
   /** Check distributed API mode. Interface faces OR nodes need to be set by the
@@ -287,15 +287,15 @@ int PMMG_preprocessMesh_distributed( PMMG_pParMesh parmesh )
   MMG3D_Set_commonFunc();
 
   /** Mesh scaling and quality histogram */
-  if ( !MMG5_scaleMesh(mesh,met,NULL) ) {
+  if ( !MMG5_scaleMesh(mesh,met,ls) ) {
     return PMMG_LOWFAILURE;
   }
 
   /* Set mmg3d  function pointers here to assign dosol */
   MMG3D_setfunc(mesh,met);
   PMMG_setfunc(parmesh);
 
-  /** specific meshing */
+  /** Specific meshing */
   if ( mesh->info.optim && !met->np ) {
     // Warning: doSol would need a clean // implementation along interfaces
     if ( !MMG3D_doSol(mesh,met) ) {
@@ -319,6 +319,8 @@ int PMMG_preprocessMesh_distributed( PMMG_pParMesh parmesh )
     mesh->info.sethmax = 1;
   }
 
+  /* 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;
   }
@@ -329,11 +331,12 @@ int PMMG_preprocessMesh_distributed( PMMG_pParMesh parmesh )
   }
 
   /** Mesh unscaling */
-  if ( !MMG5_unscaleMesh(mesh,met,NULL) ) {
+  if ( !MMG5_unscaleMesh(mesh,met,ls) ) {
     return PMMG_STRONGFAILURE;
   }
 
-  /** Mesh analysis I: check triangles, create xtetras */
+  /** Mesh analysis I: Needed to create communicators
+   *  Check triangles, create xtetras */
   if ( parmesh->myrank < parmesh->info.npartin ) {
     if ( !PMMG_analys_tria(parmesh,mesh) ) {
       return PMMG_STRONGFAILURE;
@@ -348,7 +351,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 */
+       * to edges and nodes - Why we do that? */
       if ( parmesh->myrank < parmesh->info.npartin ) {
         PMMG_tria2elmFace_coords( parmesh );
       }
@@ -376,6 +379,38 @@ int PMMG_preprocessMesh_distributed( PMMG_pParMesh parmesh )
       break;
   }
 
+  /** Discretization of the isovalue  */
+  if (mesh->info.iso) {
+
+    /* Destroy adja and adjat */
+    MMG5_DEL_MEM(mesh,mesh->adja);
+    MMG5_DEL_MEM(mesh,mesh->adjt);
+
+    tim = 1;
+    chrono(ON,&(ctim[tim]));
+    if ( parmesh->info.imprim > PMMG_VERB_VERSION ) {
+      fprintf(stdout,"\n  -- PHASE 1a: ISOVALUE DISCRETIZATION     \n");
+      fprintf(stdout,"  --    under development     \n");
+    }
+    if ( !PMMG_ls(parmesh,mesh,ls,met) ) {
+      return PMMG_STRONGFAILURE;
+    }
+    chrono(OFF,&(ctim[tim]));
+    printim(ctim[tim].gdif,stim);
+    if ( parmesh->info.imprim > PMMG_VERB_VERSION ) {
+      fprintf(stdout,"\n  -- PHASE 1a COMPLETED     %s\n",stim);
+    }
+
+    /** Mesh analysis Ib : After LS discretization
+     * Check triangles, create xtetras */
+    if ( parmesh->myrank < parmesh->info.npartin ) {
+      if ( !PMMG_analys_tria(parmesh,mesh) ) {
+        return PMMG_STRONGFAILURE;
+      }
+    }
+  }
+
+  /** Mesh analysis II: Perform surface analysis */
   if ( parmesh->myrank < parmesh->info.npartin ) {
     if ( !PMMG_analys(parmesh,mesh,parmesh->info.read_comm) ) {
       return PMMG_STRONGFAILURE;
@@ -1727,8 +1762,16 @@ int PMMG_parmmg_centralized(PMMG_pParMesh parmesh) {
 }
 
 int PMMG_parmmglib_distributed(PMMG_pParMesh parmesh) {
+    return PMMG_parmmg_distributed(parmesh);
+}
+
+int PMMG_parmmgls_distributed(PMMG_pParMesh parmesh) {
+  return PMMG_parmmg_distributed(parmesh);
+}
+
+int PMMG_parmmg_distributed(PMMG_pParMesh parmesh) {
   MMG5_pMesh       mesh;
-  MMG5_pSol        met;
+  MMG5_pSol        met,ls;
   int              ier,iresult,ierlib;
   mytime           ctim[TIMEMAX];
   int8_t           tim;
@@ -1790,7 +1833,8 @@ int PMMG_parmmglib_distributed(PMMG_pParMesh parmesh) {
     ier  = PMMG_preprocessMesh_distributed( parmesh );
     mesh = parmesh->listgrp[0].mesh;
     met  = parmesh->listgrp[0].met;
-    if ( (ier==PMMG_STRONGFAILURE) && (parmesh->nprocs == parmesh->info.npartin) && MMG5_unscaleMesh( mesh, met, NULL ) ) {
+    ls   = parmesh->listgrp[0].ls;
+    if ( (ier==PMMG_STRONGFAILURE) && (parmesh->nprocs == parmesh->info.npartin) && MMG5_unscaleMesh( mesh, met, ls ) ) {
       ier = PMMG_LOWFAILURE;
     }
   }
@@ -1834,7 +1878,7 @@ int PMMG_parmmglib_distributed(PMMG_pParMesh parmesh) {
   chrono(ON,&(ctim[tim]));
   if ( parmesh->info.imprim > PMMG_VERB_VERSION ) {
     fprintf( stdout,"\n  -- PHASE 2 : %s MESHING\n",
-             parmesh->listgrp[0].met->size < 6 ? "ISOTROPIC" : "ANISOTROPIC" );
+            parmesh->listgrp[0].met->size < 6 ? "ISOTROPIC" : "ANISOTROPIC" );
   }
 
   ier = PMMG_parmmglib1(parmesh);