Feature/ls skeleton centralized

cmake/testing/pmmg_tests.cmake

@@ -316,18 +316,58 @@ IF( BUILD_TESTING )
   #####        Test isovalue mode - ls discretization
   #####
   ###############################################################################
+  # The following tests were to test the implementation of the ls option
+  # Not relevant anymore
+  # foreach( NP 1 2 4 8 )
+  #   add_test( NAME ls-arg-option-${NP}
+  #     COMMAND ${MPIEXEC} ${MPI_ARGS} ${MPIEXEC_NUMPROC_FLAG} ${NP} $<TARGET_FILE:${PROJECT_NAME}>
+  #     ${CI_DIR}/LevelSet/3D-cube.mesh
+  #     -ls 0.01
+  #     -sol ${CI_DIR}/LevelSet/3D-cube-ls.sol
+  #     -out ${CI_DIR_RESULTS}/${MESH}-${NP}.o.mesh)
+  #   set(lsNotImplemented "## Error: level-set discretisation unavailable")
+  #   set_property(TEST ls-arg-option-${NP}
+  #     PROPERTY PASS_REGULAR_EXPRESSION "${lsNotImplemented}")
+  # endforeach()
+
+  #--------------------------------
+  #--- CENTRALIZED INPUT (CenIn)
+  #--------------------------------
+  # Tests of ls discretization for centralized mesh input
   foreach( NP 1 2 4 8 )
-    add_test( NAME ls-arg-option-${NP}
+    add_test( NAME ls-CenIn-${NP}
       COMMAND ${MPIEXEC} ${MPI_ARGS} ${MPIEXEC_NUMPROC_FLAG} ${NP} $<TARGET_FILE:${PROJECT_NAME}>
       ${CI_DIR}/LevelSet/3D-cube.mesh
-      -ls 0.01
+      -ls 0.0
       -sol ${CI_DIR}/LevelSet/3D-cube-ls.sol
-      -out ${CI_DIR_RESULTS}/${MESH}-${NP}.o.mesh)
-    set(lsNotImplemented "## Error: level-set discretisation unavailable")
-    set_property(TEST ls-arg-option-${NP}
-      PROPERTY PASS_REGULAR_EXPRESSION "${lsNotImplemented}")
+      -out ${CI_DIR_RESULTS}/3D-cube-ls-CenIn-${NP}.o.mesh)
   endforeach()
 
+  # Tests for ls + met for centralized mesh input
+  foreach( NP 1 2 4 8 )
+  add_test( NAME ls-CenIn-met-${NP}
+    COMMAND ${MPIEXEC} ${MPI_ARGS} ${MPIEXEC_NUMPROC_FLAG} ${NP} $<TARGET_FILE:${PROJECT_NAME}>
+    ${CI_DIR}/LevelSet/3D-cube.mesh
+    -ls 0.0
+    -sol ${CI_DIR}/LevelSet/3D-cube-ls.sol
+    -met ${CI_DIR}/LevelSet/3D-cube-metric.sol
+    -out ${CI_DIR_RESULTS}/3D-cube-ls-CenIn-met-${NP}.o.mesh)
+  endforeach()
+
+  # Tests of pvtu output when ls mode
+  foreach( NP 1 2 4 8 )
+    add_test( NAME ls-CenIn-DisOut-${NP}
+      COMMAND ${MPIEXEC} ${MPI_ARGS} ${MPIEXEC_NUMPROC_FLAG} ${NP} $<TARGET_FILE:${PROJECT_NAME}>
+      ${CI_DIR}/LevelSet/3D-cube.mesh
+      -ls 0.0
+      -sol ${CI_DIR}/LevelSet/3D-cube-ls.sol
+      -out ${CI_DIR_RESULTS}/3D-cube-ls-CenIn-DisOut-${NP}.o.pvtu)
+  endforeach()
+
+  #--------------------------------
+  #--- DISTRIBUTED INPUT (DisIn)
+  #--------------------------------
+
 
   ###############################################################################
   #####

src/API_functionsf_pmmg.c

@@ -1051,6 +1051,16 @@ FORTRAN_NAME(PMMG_PARMMGLIB_CENTRALIZED,pmmg_parmmglib_centralized,
   return;
 }
 
+/**
+ * See \ref PMMG_parmmgls_centralized function in \ref libparmmg.h file.
+ */
+FORTRAN_NAME(PMMG_PARMMGLS_CENTRALIZED,pmmg_parmmgls_centralized,
+             (PMMG_pParMesh *parmesh,int* retval),
+             (parmesh,retval)) {
+  *retval = PMMG_parmmgls_centralized(*parmesh);
+  return;
+}
+
 /**
  * See \ref PMMG_loadMesh function in \ref libparmmg.h file.
  */

src/libparmmg.c

@@ -55,7 +55,7 @@ extern int (*PMMG_interp2bar)(MMG5_pMesh mesh,MMG5_pSol met,MMG5_pSol oldMet,MMG
 int PMMG_check_inputData(PMMG_pParMesh parmesh)
 {
   MMG5_pMesh mesh;
-  MMG5_pSol  met;
+  MMG5_pSol  met,ls;
   int        k;
 
   if ( parmesh->info.imprim > PMMG_VERB_VERSION )
@@ -71,14 +71,14 @@ int PMMG_check_inputData(PMMG_pParMesh parmesh)
               "  ## Error: lagrangian mode unavailable (MMG3D_IPARAM_lag):\n");
       return 0;
     } else if ( mesh->info.iso ) {
-      fprintf(stderr,"  ## Error: level-set discretisation unavailable"
-              " (MMG3D_IPARAM_iso):\n");
-      return 0;
+      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");
       return 0;
     }
+
     /* specific meshing */
     if ( met->np ) {
       if ( mesh->info.optim ) {
@@ -128,18 +128,23 @@ int PMMG_check_inputData(PMMG_pParMesh parmesh)
 int PMMG_preprocessMesh( PMMG_pParMesh parmesh )
 {
   MMG5_pMesh mesh;
-  MMG5_pSol  met;
+  MMG5_pSol  met,ls;
+  int8_t     tim;
+  char       stim[32];
+  mytime     ctim[TIMEMAX];
+
 
   mesh = parmesh->listgrp[0].mesh;
   met  = parmesh->listgrp[0].met;
+  ls   = parmesh->listgrp[0].ls;
 
   assert ( ( mesh != NULL ) && ( met != NULL ) && "Preprocessing empty args");
 
   /** Function setters (must be assigned before quality computation) */
   MMG3D_Set_commonFunc();
 
   /** Mesh scaling and quality histogram */
-  if ( !MMG5_scaleMesh(mesh,met,NULL) ) {
+  if ( !MMG5_scaleMesh(mesh,met,ls) ) {
     return PMMG_LOWFAILURE;
   }
 
@@ -177,17 +182,45 @@ int PMMG_preprocessMesh( PMMG_pParMesh parmesh )
     return PMMG_STRONGFAILURE;
   }
 
+  /* Discretization of the isovalue  */
+  if (mesh->info.iso) {
+    tim = 1;
+    chrono(ON,&(ctim[tim]));
+    if ( parmesh->info.imprim > PMMG_VERB_VERSION ) {
+      fprintf(stdout,"\n  -- PHASE 1a: ISOVALUE DISCRETIZATION     \n");
+    }
+    if ( !MMG3D_mmg3d2(mesh,ls,met) ) {
+      return PMMG_STRONGFAILURE;
+    }
+    chrono(OFF,&(ctim[tim]));
+    printim(ctim[tim].gdif,stim);
+    if ( parmesh->info.imprim > PMMG_VERB_VERSION ) {
+      fprintf(stdout,"  -- PHASE 1a COMPLETED     %s\n",stim);
+    }
+  }
+
   /** Mesh analysis */
   if ( !MMG3D_analys(mesh) ) {
     return PMMG_STRONGFAILURE;
   }
 
+  /* Check if the LS has led to a non-manifold topology */
+  if ( mesh->info.iso && !MMG3D_chkmani(mesh) ) {
+    fprintf(stderr,"\n  ## LS discretization: non-manifold initial topology. Exit program.\n");
+    return PMMG_STRONGFAILURE;
+  }
+  else {
+    if ( parmesh->info.imprim > PMMG_VERB_VERSION ) {
+      fprintf(stdout,"       LS discretization OK: no non-manifold topology.\n");
+    }
+  }
+
   if ( parmesh->info.imprim0 > PMMG_VERB_ITWAVES && (!mesh->info.iso) && met->m ) {
     PMMG_prilen(parmesh,0,1);
   }
 
   /** Mesh unscaling */
-  if ( !MMG5_unscaleMesh(mesh,met,NULL) ) {
+  if ( !MMG5_unscaleMesh(mesh,met,ls) ) {
     return PMMG_STRONGFAILURE;
   }
 
@@ -337,7 +370,7 @@ int PMMG_preprocessMesh_distributed( PMMG_pParMesh parmesh )
 
 int PMMG_distributeMesh_centralized_timers( PMMG_pParMesh parmesh,mytime *ctim ) {
   MMG5_pMesh    mesh;
-  MMG5_pSol     met;
+  MMG5_pSol     met,ls;
   int           ier,iresult;
   int8_t        tim;
   char          stim[32];
@@ -378,7 +411,8 @@ int PMMG_distributeMesh_centralized_timers( PMMG_pParMesh parmesh,mytime *ctim )
 
     mesh = parmesh->listgrp[0].mesh;
     met  = parmesh->listgrp[0].met;
-    if ( (ier==PMMG_STRONGFAILURE) && MMG5_unscaleMesh( mesh, met, NULL ) ) {
+    ls   = parmesh->listgrp[0].ls;
+    if ( (ier==PMMG_STRONGFAILURE) && MMG5_unscaleMesh( mesh, met, ls ) ) {
       ier = PMMG_LOWFAILURE;
     }
 
@@ -1580,9 +1614,17 @@ int PMMG_parmmglib_post(PMMG_pParMesh parmesh) {
 }
 
 int PMMG_parmmglib_centralized(PMMG_pParMesh parmesh) {
+  return PMMG_parmmg_centralized(parmesh);
+}
+
+int PMMG_parmmgls_centralized(PMMG_pParMesh parmesh) {
+  return PMMG_parmmg_centralized(parmesh);
+}
+
+int PMMG_parmmg_centralized(PMMG_pParMesh parmesh) {
   PMMG_pGrp     grp;
   MMG5_pMesh    mesh;
-  MMG5_pSol     met;
+  MMG5_pSol     met,ls;
   int           ier;
   int           ierlib;
   mytime        ctim[TIMEMAX];
@@ -1619,6 +1661,7 @@ int PMMG_parmmglib_centralized(PMMG_pParMesh parmesh) {
   grp    = &parmesh->listgrp[0];
   mesh   = grp->mesh;
   met    = grp->met;
+  ls     = grp->ls;
 
   /** Remeshing */
   tim = 3;

src/libparmmg.h

@@ -162,6 +162,46 @@ int PMMG_parmmglib_distributed(PMMG_pParMesh parmesh);
  **/
 int PMMG_parmmglib_centralized(PMMG_pParMesh parmesh);
 
+/**
+ * \param parmesh pointer toward the parmesh structure (boundary entities are
+ * stored into MMG5_Tria, MMG5_Edge... structures)
+ *
+ * \return \ref PMMG_SUCCESS if success, \ref PMMG_LOWFAILURE if fail but we can
+ * return a centralized and unscaled mesh or \ref PMMG_STRONGFAILURE if fail and
+ * we can't return a centralized and unscaled mesh.
+ *
+ * Main program for the parallel isovalue discretisation library for centralized
+ * meshes
+ *
+ * \remark Fortran interface:
+ * >   SUBROUTINE PMMG_parmmgls_centralized(parmesh,retval)\n
+ * >     MMG5_DATA_PTR_T,INTENT(INOUT) :: parmesh\n
+ * >     INTEGER, INTENT(OUT)          :: retval\n
+ * >   END SUBROUTINE\n
+ *
+ **/
+int PMMG_parmmgls_centralized(PMMG_pParMesh parmesh);
+
+/**
+ * \param parmesh pointer toward the parmesh structure (boundary entities are
+ * stored into MMG5_Tria, MMG5_Edge... structures)
+ *
+ * \return \ref PMMG_SUCCESS if success, \ref PMMG_LOWFAILURE if fail but we can
+ * return a centralized and unscaled mesh or \ref PMMG_STRONGFAILURE if fail and
+ * we can't return a centralized and unscaled mesh.
+ *
+ * Main program for the parallel isovalue discretisation library and remesh
+ * library for centralized meshes
+ *
+ * \remark Fortran interface:
+ * >   SUBROUTINE PMMG_parmmg_centralized(parmesh,retval)\n
+ * >     MMG5_DATA_PTR_T,INTENT(INOUT) :: parmesh\n
+ * >     INTEGER, INTENT(OUT)          :: retval\n
+ * >   END SUBROUTINE\n
+ *
+ **/
+int PMMG_parmmg_centralized(PMMG_pParMesh parmesh);
+
 /* init file names */
 /**
  * \param parmesh pointer toward a parmesh structure.

src/libparmmg1.c

@@ -707,7 +707,6 @@ int PMMG_parmmglib1( PMMG_pParMesh parmesh )
 
         /* renumerotation if available: no need to renum the field here (they
          * will be interpolated) */
-        assert ( mesh->npi==mesh->np );
         if ( permNodGlob ) {
           if ( !MMG5_scotchCall(mesh,met,NULL,permNodGlob) )
           {

src/mpiunpack_pmmg.c

@@ -92,10 +92,12 @@ int PMMG_mpiunpack_meshSizes ( PMMG_pParMesh parmesh,PMMG_pGrp listgrp,int igrp,
   ier_grp = MMG3D_Init_mesh(MMG5_ARG_start,
                             MMG5_ARG_ppMesh,&(grp->mesh),
                             MMG5_ARG_ppMet ,&(grp->met),
+                            MMG5_ARG_ppLs  ,&(grp->ls),
                             MMG5_ARG_end);
 
   mesh = grp->mesh;
   met  = grp->met;
+  ls   = grp->ls;
 
   /* Set maximum memory */
   mesh->memMax = parmesh->memGloMax;

src/parmmg.c

@@ -275,20 +275,35 @@ int main( int argc, char *argv[] )
       }
     }
     /* In iso mode: read metric if any */
-    if ( grp->mesh->info.iso && parmesh->metin ) {
-      if ( !distributedInput ) {
-        iermesh = PMMG_loadMet_centralized( parmesh, parmesh->metin );
+    if ( grp->mesh->info.iso) {
+      if ( parmesh->metin ) {
+        if ( !distributedInput ) {
+          iermesh = PMMG_loadMet_centralized( parmesh, parmesh->metin );
+        }
+        else {
+          int ier_loc = PMMG_loadMet_distributed( parmesh, parmesh->metin );
+          MPI_Allreduce( &ier_loc, &iermesh, 1, MPI_INT, MPI_MIN, parmesh->comm);
+        }
+        if ( -1 == iermesh ) {
+          if ( rank == parmesh->info.root ) {
+            fprintf(stderr,"\n  ## ERROR: UNABLE TO LOAD METRIC.\n");
+          }
+          ier = 0;
+          goto check_mesh_loading;
+        }
       }
       else {
-        int ier_loc = PMMG_loadMet_distributed( parmesh, parmesh->metin );
-        MPI_Allreduce( &ier_loc, &iermesh, 1, MPI_INT, MPI_MIN, parmesh->comm);
-      }
-      if ( -1 == iermesh ) {
-        if ( rank == parmesh->info.root ) {
-          fprintf(stderr,"\n  ## ERROR: UNABLE TO LOAD METRIC.\n");
+        /* Give a name to the metric if not provided for distributed metric output */
+        if ( !MMG5_Set_inputSolName(grp->mesh,grp->met,"") ) {
+          fprintf(stdout,"  ## WARNING: Unable to give a name to the metric.\n");
+        }
+        else {
+          ier = PMMG_Set_name(parmesh,&parmesh->metin,grp->met->namein,"mesh.sol");
+          if (!ier) {
+            fprintf(stdout,"  ## ERROR: Unable to give a name to the metric.\n");
+            PMMG_RETURN_AND_FREE( parmesh, PMMG_LOWFAILURE );
+          }
         }
-        ier = 0;
-        goto check_mesh_loading;
       }
     }
 
@@ -343,7 +358,7 @@ int main( int argc, char *argv[] )
   }
   else if ( !distributedInput ) {
     /* Parallel remeshing starting from a centralized mesh */
-    ier = PMMG_parmmglib_centralized(parmesh);
+    ier = PMMG_parmmg_centralized(parmesh);
   }
   else {
     /* Parallel remeshing starting from a distributed mesh */