Merge pull request MmgTools#102 from MmgTools/feature/fix-edgeLength

Feature/fix edge length
laetitia-m · Dec 8, 2023 · eaebcbe · eaebcbe
2 parents 4f0a95d + 5ef4be1
commit eaebcbe
Show file tree

Hide file tree

Showing 5 changed files with 65 additions and 25 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -252,11 +252,12 @@ IF ( DOWNLOAD_MMG )
 
   EXTERNALPROJECT_ADD ( Mmg
     GIT_REPOSITORY https://github.com/MmgTools/mmg.git
-    GIT_TAG 41b2e777aea1dfee594f02fa7fcdbe8630d1001d
+    GIT_TAG 3a5e3abc38f7f100a481bcebe8aa23e105d7c381
     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}
-    -DBUILD=MMG3D -DBUILD_SHARED_LIBS=${LIBPARMMG_SHARED}  -DUSE_VTK=${USE_VTK}
+    -DBUILD=MMG3D -DBUILD_SHARED_LIBS=${LIBPARMMG_SHARED}
+    -DUSE_VTK=${USE_VTK} -DPMMG_CALL=1
     -DMMG_INSTALL_PRIVATE_HEADERS=ON
     -DUSE_POINTMAP=${USE_POINTMAP} -DCMAKE_INSTALL_PREFIX=../Mmg-install)
 

diff --git a/cmake/modules/optional-dependencies.cmake b/cmake/modules/optional-dependencies.cmake
@@ -108,7 +108,7 @@ IF ( NOT DEFINED USE_VTK OR USE_VTK STREQUAL "" OR USE_VTK MATCHES " +" OR USE_V
   # Before v9
   FIND_PACKAGE(VTK QUIET)
   IF ( VTK_FOUND )
-    message (STATUS "VTK_VERSION: ${VTK_VERSION}")
+    message (STATUS "ParMmg - VTK_VERSION: ${VTK_VERSION}")
     IF (VTK_VERSION VERSION_LESS "9.0.0")
       find_package(VTK  COMPONENTS
         vtkCommonCore
@@ -137,15 +137,15 @@ IF ( NOT DEFINED USE_VTK OR USE_VTK STREQUAL "" OR USE_VTK MATCHES " +" OR USE_V
   ELSEIF ( USE_VTK )
     # USE_VTK is not empty so user explicitely ask for VTK...
     # but it is not found: raise an error
-    MESSAGE(FATAL_ERROR "VTK library not found.")
+    MESSAGE(FATAL_ERROR "ParMmg - VTK library not found.")
   ENDIF()
 ENDIF()
 
 IF ( VTK_FOUND AND NOT USE_VTK MATCHES OFF)
 
   add_definitions(-DUSE_VTK)
 
-  MESSAGE ( STATUS "Compilation with VTK: add vtk, vtp and vtu I/O." )
+  MESSAGE ( STATUS "ParMmg - Compilation with VTK: add vtk, vtp and vtu I/O." )
 
   IF( "${VTK_MAJOR_VERSION}.${VTK_MINOR_VERSION}" LESS 8.90 )
     INCLUDE ( ${VTK_USE_FILE} )

diff --git a/src/libparmmg.c b/src/libparmmg.c
@@ -130,6 +130,8 @@ int PMMG_preprocessMesh( PMMG_pParMesh parmesh )
   char       stim[32];
   mytime     ctim[TIMEMAX];
 
+  /* Chrono initialization */
+  tminit(ctim,TIMEMAX);
 
   mesh = parmesh->listgrp[0].mesh;
   met  = parmesh->listgrp[0].met;
@@ -243,6 +245,9 @@ int PMMG_preprocessMesh_distributed( PMMG_pParMesh parmesh )
   mytime     ctim[TIMEMAX];
   int ier = PMMG_SUCCESS;
 
+  /* Chrono initialization */
+  tminit(ctim,TIMEMAX);
+
   mesh = parmesh->listgrp[0].mesh;
   met  = parmesh->listgrp[0].met;
   ls   = parmesh->listgrp[0].ls;

diff --git a/src/moveinterfaces_pmmg.c b/src/moveinterfaces_pmmg.c
@@ -1281,6 +1281,10 @@ int PMMG_set_ifcDirection( PMMG_pParMesh parmesh,int **displsgrp,int **mapgrp )
   for( k=0; k<nproc; k++ )
     ngrps[k] = (*displsgrp)[k+1]-(*displsgrp)[k];
 
+  // Remark: with open-mpi (v4.1.2) the next allgatherv raises a valgrind error
+  // on rank 0 due to the use of memcpy with overlapping locations. It seems to
+  // be an error in the open-mpi implementation and to have no consequences on
+  // the results
   MPI_CHECK(
     MPI_Allgatherv(MPI_IN_PLACE,0,MPI_DATATYPE_NULL,
                    &(*mapgrp)[0],ngrps,&(*displsgrp)[0],MPI_INT,

diff --git a/src/quality_pmmg.c b/src/quality_pmmg.c
@@ -24,6 +24,7 @@
 #include "parmmg.h"
 #include <stddef.h>
 #include "inlined_functions_3d_private.h"
+#include "mmgexterns_private.h"
 
 typedef struct {
   double min;
@@ -371,7 +372,7 @@ int PMMG_qualhisto( PMMG_pParMesh parmesh, int opt, int isCentral, MPI_Comm comm
  * read_comm communicator (i.e. the communicator used to provide the inputs).
  * For all ather calls, comm has to be the communicator to use for computations.
  *
- * \return 0 if fail, 1 otherwise.
+ * \return 2 without metric, 0 if fail, 1 otherwise.
  *
  * Compute the required information to print the length histogram
  *
@@ -405,6 +406,12 @@ int PMMG_computePrilen( PMMG_pParMesh parmesh,MMG5_pMesh mesh, MMG5_pSol met, do
   *amin = *amax = *bmin = *bmax = 0;
   *nullEdge = 0;
 
+  if ( (!met) || (!met->m) ) {
+    /* the functions that computes the edge length cannot be called without an
+     * allocated metric */
+    return 2;
+  }
+
   /* Hash parallel edges in the mesh */
   if ( PMMG_hashPar(mesh,&hpar) != PMMG_SUCCESS ) return 0;
 
@@ -469,11 +476,18 @@ int PMMG_computePrilen( PMMG_pParMesh parmesh,MMG5_pMesh mesh, MMG5_pSol met, do
       /* Remove edge from hash ; ier = 1 if edge has been found */
       ier = MMG5_hashPop(&hash,np,nq);
       if( ier ) {
-        if ( (!metRidTyp) && met->size==6 && met->m ) {
+        assert ( met->m );
+        if ( (!metRidTyp) && met->size==6 ) {
+          assert ( met->m );
+          /* We pass here if metric is aniso without metRidTyp */
           len = MMG5_lenSurfEdg33_ani(mesh,met,np,nq,(tag & MG_GEO));
         }
-        else
-          len = MMG5_lenSurfEdg_iso(mesh,met,np,nq,0);
+        else {
+          /* We pass here if metric is aniso with metRidTyp or iso with
+           * allocated metric. Note that the lenSurfEdg function segfault if called
+           * with met==NULL or met->m==NULL */
+          len = MMG5_lenSurfEdg(mesh,met,np,nq,0);
+        }
 
 
         if ( !len ) {
@@ -534,6 +548,7 @@ int PMMG_computePrilen( PMMG_pParMesh parmesh,MMG5_pMesh mesh, MMG5_pSol met, do
       /* Remove edge from hash ; ier = 1 if edge has been found */
       ier = MMG5_hashPop(&hash,np,nq);
       if( ier ) {
+        assert ( met->m );
         if ( (!metRidTyp) && met->size==6 && met->m ) {
           len = MMG5_lenedg33_ani(mesh,met,ia,pt);
         }
@@ -656,21 +671,35 @@ int PMMG_prilen( PMMG_pParMesh parmesh, int8_t metRidTyp, int isCentral, MPI_Com
   if ( parmesh->ngrp==1 ) {
     mesh = parmesh->listgrp[0].mesh;
     met = parmesh->listgrp[0].met;
-    if ( met && met->m ) {
-      if( isCentral )
-        ier = MMG3D_computePrilen( mesh, met,
-                                   &lenStats.avlen, &lenStats.lmin,
-                                   &lenStats.lmax, &lenStats.ned, &lenStats.amin,
-                                   &lenStats.bmin, &lenStats.amax, &lenStats.bmax,
-                                   &lenStats.nullEdge, metRidTyp, &bd,
-                                   lenStats.hl );
-      else
-        ier = PMMG_computePrilen( parmesh, mesh, met,
-                                  &lenStats.avlen, &lenStats.lmin,
-                                  &lenStats.lmax, &lenStats.ned, &lenStats.amin,
-                                  &lenStats.bmin, &lenStats.amax, &lenStats.bmax,
-                                  &lenStats.nullEdge, metRidTyp, &bd,
-                                  lenStats.hl,comm );
+    if( isCentral ) {
+      /* If metric is not allocated or if hash table alloc fails, the next
+       * function returns 0, which allows to detect that we cannot print the
+       * edge length histo. */
+      ier = MMG3D_computePrilen( mesh, met,
+                                 &lenStats.avlen, &lenStats.lmin,
+                                 &lenStats.lmax, &lenStats.ned, &lenStats.amin,
+                                 &lenStats.bmin, &lenStats.amax, &lenStats.bmax,
+                                 &lenStats.nullEdge, metRidTyp, &bd,
+                                 lenStats.hl );
+    }
+    else {
+      /* The next function returns 0 if the hash table alloc fails and 2 if
+       * called without metric (in this case we are not able to compute the edge
+       * lengths). It allows to detect:
+       *   - if we can't print the histo due to an alloc error (ier = 0 on 1 MPI
+       *     process at least)
+       *   - if we can't print the histo because the metric is not allocated
+       *     (ier=2 on all the MPI process)
+       *   - if we can print the histo (metric is allocated on at least 1 MPI
+       *     process and no MPI process fail, thus ier is at least 1 on all the
+       *     MPI proc but may be 2 on some of them)
+       */
+      ier = PMMG_computePrilen( parmesh, mesh, met,
+                                &lenStats.avlen, &lenStats.lmin,
+                                &lenStats.lmax, &lenStats.ned, &lenStats.amin,
+                                &lenStats.bmin, &lenStats.amax, &lenStats.bmax,
+                                &lenStats.nullEdge, metRidTyp, &bd,
+                                lenStats.hl,comm );
     }
   }
 
@@ -679,7 +708,8 @@ int PMMG_prilen( PMMG_pParMesh parmesh, int8_t metRidTyp, int isCentral, MPI_Com
   else
     MPI_Reduce( &ier, &ieresult,1, MPI_INT, MPI_MIN, parmesh->info.root, comm );
 
-  if ( !ieresult ) {
+  if ( (ieresult==0) || ieresult==2 ) {
+    /* We are not able to print the histogram */
     MPI_Type_free( &mpi_lenStats_t );
     MPI_Op_free( &mpi_lenStats_op );
     return 0;