Modifications dedicated to Nextsim remeshing

src/common/libmmgtypes.h

@@ -524,9 +524,12 @@ typedef struct {
   MMG5_pPar     par;
   double        dhd,hmin,hmax,hsiz,hgrad,hgradreq,hausd;
   double        min[3],max[3],delta,ls,lxreg,rmc;
+  double        limit_angle; /*< Angle threshold for modifying triangles or not */
   MMG5_int      *br; /*!< list of based references to which an implicit surface can be attached */
   MMG5_int      isoref; /*!< isovalue reference in ls mode */
   MMG5_int      nsd; /*!< index of subdomain to save (0 by default == all subdomains are saved) */
+  int           isotropic; /*!< force the use of some isotropic functions */
+  int           bdy_adaptation;
   int           mem,npar,npari;
   int           nbr,nbri; /*!< number of based references for level-set (BC to which a material can be attached) */
   int           opnbdy; /*!< floating surfaces */

src/common/scalem.c

@@ -745,7 +745,7 @@ int MMG5_unscaleMesh(MMG5_pMesh mesh,MMG5_pSol met,MMG5_pSol sol) {
   mesh->info.delta = 1.;
   mesh->info.min[0]= 0.;
   mesh->info.min[1]= 0.;
-  mesh->info.min[2]= 0.;
+  if (mesh->dim == 3) mesh->info.min[2]= 0.;
 
   /* de-normalize metric */
   if ( !(met && met->np && met->m) )  return 1;

src/mmg2d/API_functions_2d.c

@@ -101,11 +101,16 @@ void MMG2D_Init_parameters(MMG5_pMesh mesh) {
   /* default values for doubles */
   /* level set value */
   mesh->info.ls       = MMG5_LS;
-/* xreg relaxation parameter value */
+  /* xreg relaxation parameter value */
   mesh->info.lxreg    = MMG5_XREG;
-
+  /* [0/1]   ,avoid/enforce istropic remeshing even with anisotropic metric */
+  mesh->info.isotropic = MMG5_OFF;
+  /* limit angle to avoid remeshing some good triangles */
+  mesh->info.limit_angle = 5.*atan(1.);
   /* Ridge detection */
   mesh->info.dhd      = MMG5_ANGEDG;
+
+  mesh->info.bdy_adaptation = MMG5_OFF;
 }
 
 
@@ -147,6 +152,9 @@ int MMG2D_Set_iparameter(MMG5_pMesh mesh, MMG5_pSol sol, int iparam, MMG5_int va
       mesh->info.dhd    = MMG5_ANGEDG;
     }
     break;
+  case MMG2D_IPARAM_bdy_adaptation:
+    mesh->info.bdy_adaptation = val;
+    break;
   case MMG2D_IPARAM_nofem :
     mesh->info.setfem = (val==1)? 0 : 1;
     break;
@@ -162,6 +170,9 @@ int MMG2D_Set_iparameter(MMG5_pMesh mesh, MMG5_pSol sol, int iparam, MMG5_int va
   case MMG2D_IPARAM_isosurf :
     mesh->info.isosurf = val;
     break;
+  case MMG2D_IPARAM_isotropic :
+    mesh->info.isotropic = val;
+    break;
   case MMG2D_IPARAM_lag :
 #ifdef USE_ELAS
     if ( val < 0 || val > 2 )
@@ -349,6 +360,15 @@ int MMG2D_Set_dparameter(MMG5_pMesh mesh, MMG5_pSol sol, int dparam, double val)
     else
       mesh->info.hausd    = val;
     break;
+  case MMG2D_DPARAM_hmin_factor :
+    mesh->info.min[2] = val;
+    break;
+  case MMG2D_DPARAM_hmax_factor :
+    mesh->info.max[2] = val;
+    break;
+  case MMG2D_DPARAM_limit_angle :
+    mesh->info.limit_angle = val;
+    break;
   case MMG2D_DPARAM_ls :
     mesh->info.ls       = val;
     break;

src/mmg2d/libmmg2d.c

@@ -60,7 +60,7 @@ void MMG2D_Set_commonFunc(void) {
     return;
 }
 
-int MMG2D_mmg2dlib(MMG5_pMesh mesh,MMG5_pSol met) {
+int MMG2D_mmg2dlib(MMG5_pMesh mesh,MMG5_pSol met,double* velocity) {
   MMG5_pSol sol=NULL; // unused
   mytime    ctim[TIMEMAX];
   char      stim[32];
@@ -219,7 +219,7 @@ int MMG2D_mmg2dlib(MMG5_pMesh mesh,MMG5_pSol met) {
     fprintf(stdout,"\n  -- PHASE 2 : %s MESHING\n",met->size < 3 ? "ISOTROPIC" : "ANISOTROPIC");
 
   /* Mesh improvement */
-  if ( !MMG2D_mmg2d1n(mesh,met) ) {
+  if ( !MMG2D_mmg2d1n(mesh,met,velocity) ) {
     if ( !MMG5_unscaleMesh(mesh,met,NULL) )  _LIBMMG5_RETURN(mesh,met,sol,MMG5_STRONGFAILURE);
     MMG2D_RETURN_AND_PACK(mesh,met,sol,MMG5_LOWFAILURE);
   }
@@ -306,6 +306,7 @@ int MMG2D_mmg2dmesh(MMG5_pMesh mesh,MMG5_pSol met) {
   MMG5_pSol sol=NULL; // unused
   mytime    ctim[TIMEMAX];
   char      stim[32];
+  double* velocity;
 
   assert ( mesh );
   assert ( met );
@@ -479,7 +480,7 @@ int MMG2D_mmg2dmesh(MMG5_pMesh mesh,MMG5_pSol met) {
     fprintf(stdout,"\n  -- PHASE 3 : MESH IMPROVEMENT (%s)\n",
             met->size < 3 ? "ISOTROPIC" : "ANISOTROPIC");
 
-  if ( !MMG2D_mmg2d1n(mesh,met) ) {
+  if ( !MMG2D_mmg2d1n(mesh,met,velocity) ) {
     if ( !MMG5_unscaleMesh(mesh,met,NULL) )  _LIBMMG5_RETURN(mesh,met,sol,MMG5_STRONGFAILURE);
     MMG2D_RETURN_AND_PACK(mesh,met,sol,MMG5_LOWFAILURE);
   }
@@ -526,6 +527,7 @@ int MMG2D_mmg2dls(MMG5_pMesh mesh,MMG5_pSol sol,MMG5_pSol umet) {
   mytime    ctim[TIMEMAX];
   char      stim[32];
   int8_t    mettofree = 0;
+  double*   velocity;
 
   assert ( mesh );
   assert ( sol );
@@ -758,7 +760,7 @@ int MMG2D_mmg2dls(MMG5_pMesh mesh,MMG5_pSol sol,MMG5_pSol umet) {
     fprintf(stdout,"\n  -- PHASE 3 : MESH IMPROVEMENT\n");
   }
 
-  if ( !MMG2D_mmg2d1n(mesh,met) ) {
+  if ( !MMG2D_mmg2d1n(mesh,met,velocity) ) {
     if ( mettofree ) {
       MMG5_DEL_MEM(mesh,met->m);
       MMG5_SAFE_FREE (met);
@@ -822,6 +824,7 @@ int MMG2D_mmg2dmov(MMG5_pMesh mesh,MMG5_pSol met,MMG5_pSol disp) {
   char      stim[32];
   int       ier;
   MMG5_int  k,*invalidTris;
+  double*   velocity;
 
   assert ( mesh );
   assert ( disp );
@@ -984,7 +987,7 @@ int MMG2D_mmg2dmov(MMG5_pMesh mesh,MMG5_pSol met,MMG5_pSol disp) {
       fprintf(stdout,"\n  -- PHASE 3 : MESH IMPROVEMENT\n");
     }
 
-    if ( !MMG2D_mmg2d1n(mesh,met) ) {
+    if ( !MMG2D_mmg2d1n(mesh,met,velocity) ) {
       if ( !MMG5_unscaleMesh(mesh,met,NULL) )  _LIBMMG5_RETURN(mesh,met,disp,MMG5_STRONGFAILURE);
       MMG2D_RETURN_AND_PACK(mesh,met,disp,MMG5_LOWFAILURE);
     }

src/mmg2d/libmmg2d.h

@@ -138,11 +138,16 @@ extern "C" {
     MMG2D_DPARAM_hausd,             /*!< [val], Global Hausdorff distance (on all the boundary surfaces of the mesh) */
     MMG2D_DPARAM_hgrad,             /*!< [val], Global gradation */
     MMG2D_DPARAM_hgradreq,          /*!< [val], Global gradation on required entites (advanced usage) */
+    MMG2D_DPARAM_hmin_factor,       /*!< [val], Factor on minimal edge length when limit_angle is active */
+    MMG2D_DPARAM_hmax_factor,       /*!< [val], Factor on maximal edge length when limit_angle is active */
     MMG2D_DPARAM_ls,                /*!< [val], Function value where the level set is to be discretized */
     MMG2D_DPARAM_xreg,              /*!< [val], Relaxation parameter for coordinate regularization (0<val<1) */
     MMG2D_DPARAM_rmc,               /*!< [-1/val], Remove small disconnected components in level-set mode */
     MMG2D_IPARAM_nofem,             /*!< [1/0], Do not attempt to make the mesh suitable for finite-element computations */
     MMG2D_IPARAM_isoref,            /*!< [0/n], Iso-surface boundary material reference */
+    MMG2D_IPARAM_isotropic,         /*!< [0/1], Avoid/enforce isotropic remeshing even with anisotropic metric */
+    MMG2D_DPARAM_limit_angle,       /*!< [val], Minimal angle in triangles under which remeshing is achieved */
+    MMG2D_IPARAM_bdy_adaptation,    /*!< [1/0], Enable thorough adaptation close to the boundaries (if limit_angle < Pi) */
   };
 
 /*----------------------------- function headers -----------------------------*/
@@ -2293,7 +2298,7 @@ LIBMMG2D_EXPORT int MMG2D_Free_all(const int starter,...);
  * >   END SUBROUTINE\n
  *
  */
-  LIBMMG2D_EXPORT int MMG2D_mmg2dlib(MMG5_pMesh mesh,MMG5_pSol sol);
+  LIBMMG2D_EXPORT int MMG2D_mmg2dlib(MMG5_pMesh mesh,MMG5_pSol sol,double* velocity);
 
 /**
  * \brief Main "program" for the mesh generation library.

src/mmg2d/libmmg2d_private.h

@@ -254,9 +254,9 @@ int MMG2D_norver(MMG5_pMesh,MMG5_int );
 int MMG2D_regnor(MMG5_pMesh );
 int MMG2D_regver(MMG5_pMesh );
 int MMG2D_boulen(MMG5_pMesh , MMG5_int ,int8_t ,MMG5_int *,MMG5_int *,double *);
-int MMG2D_mmg2d1n(MMG5_pMesh ,MMG5_pSol );
-int MMG2D_anatri(MMG5_pMesh ,MMG5_pSol ,int8_t );
-int MMG2D_adptri(MMG5_pMesh ,MMG5_pSol );
+int MMG2D_mmg2d1n(MMG5_pMesh ,MMG5_pSol, double*);
+int MMG2D_anatri(MMG5_pMesh ,MMG5_pSol ,int8_t , double*);
+int MMG2D_adptri(MMG5_pMesh ,MMG5_pSol, double*);
 int MMG2D_defsiz_iso(MMG5_pMesh ,MMG5_pSol );
 int MMG2D_defsiz_ani(MMG5_pMesh ,MMG5_pSol );
 int MMG2D_defmetbdy_2d(MMG5_pMesh ,MMG5_pSol ,MMG5_int ,int8_t );
@@ -265,9 +265,9 @@ MMG5_int MMG2D_grad2met_ani(MMG5_pMesh ,MMG5_pSol ,MMG5_pTria,MMG5_int,MMG5_int)
 int MMG2D_grad2metreq_ani(MMG5_pMesh ,MMG5_pSol ,MMG5_pTria,MMG5_int,MMG5_int);
 int MMG2D_gradsiz_ani(MMG5_pMesh ,MMG5_pSol );
 int MMG2D_gradsizreq_ani(MMG5_pMesh ,MMG5_pSol );
-MMG5_int MMG2D_anaelt(MMG5_pMesh ,MMG5_pSol ,int );
-MMG5_int MMG2D_colelt(MMG5_pMesh ,MMG5_pSol ,int, double );
-MMG5_int MMG2D_swpmsh(MMG5_pMesh ,MMG5_pSol ,int );
+MMG5_int MMG2D_anaelt(MMG5_pMesh ,MMG5_pSol ,int ,double*);
+MMG5_int MMG2D_colelt(MMG5_pMesh ,MMG5_pSol ,int ,double, double*);
+MMG5_int MMG2D_swpmsh(MMG5_pMesh ,MMG5_pSol ,int ,double*);
 double MMG2D_lencurv_iso(MMG5_pMesh ,MMG5_pSol ,MMG5_int ,MMG5_int );
 double MMG2D_lencurv_ani(MMG5_pMesh ,MMG5_pSol ,MMG5_int ,MMG5_int );
 int MMG2D_chkedg(MMG5_pMesh ,MMG5_int );
@@ -287,8 +287,8 @@ int MMG2D_swapar(MMG5_pMesh ,MMG5_int ,int8_t );
 int MMG5_interpmet22(MMG5_pMesh ,double *,double *,double ,double *);
 int MMG2D_intmet_iso(MMG5_pMesh ,MMG5_pSol ,MMG5_int ,int8_t ,MMG5_int ,double );
 int MMG2D_intmet_ani(MMG5_pMesh ,MMG5_pSol ,MMG5_int ,int8_t ,MMG5_int ,double );
-MMG5_int MMG2D_adpspl(MMG5_pMesh ,MMG5_pSol );
-MMG5_int MMG2D_movtri(MMG5_pMesh ,MMG5_pSol ,int ,int8_t );
+MMG5_int MMG2D_adpspl(MMG5_pMesh ,MMG5_pSol ,double*);
+MMG5_int MMG2D_movtri(MMG5_pMesh ,MMG5_pSol ,int ,int8_t ,double*);
 MMG5_int MMG2D_chkspl(MMG5_pMesh ,MMG5_pSol ,MMG5_int ,int8_t );
 int MMG2D_split1b(MMG5_pMesh ,MMG5_int ,int8_t ,MMG5_int );
 int MMG2D_movedgpt(MMG5_pMesh ,MMG5_pSol ,int ,MMG5_int *,int8_t );

src/mmg2d/libmmg2d_tools.c

@@ -148,6 +148,11 @@ int MMG2D_parsar(int argc,char *argv[],MMG5_pMesh mesh,MMG5_pSol met,MMG5_pSol s
         if ( !MMG2D_Set_solSize(mesh,met,MMG5_Vertex,0,MMG5_Tensor) )
           return 0;
         break;
+      case 'b':
+        if ( !strcmp(argv[i],"-bdy-adaptation") ) {
+          if ( !MMG2D_Set_iparameter(mesh,met,MMG2D_IPARAM_bdy_adaptation,1) )
+            return 0;
+        }
       case 'd':
         if ( !strcmp(argv[i],"-default") ) {
           mesh->mark=1;
@@ -194,6 +199,12 @@ int MMG2D_parsar(int argc,char *argv[],MMG5_pMesh mesh,MMG5_pSol met,MMG5_pSol s
           else if ( !strcmp(argv[i],"-hgrad") ) {
             param = MMG2D_DPARAM_hgrad;
           }
+          else if ( !strcmp(argv[i],"-hmin_factor") ) {
+            param = MMG2D_DPARAM_hmin_factor;
+          }
+          else if ( !strcmp(argv[i],"-hmax_factor") ) {
+            param = MMG2D_DPARAM_hmax_factor;
+          }
           else {
             /* Arg unknown by Mmg: arg starts with -h but is not known */
             MMG2D_usage(argv[0]);
@@ -235,6 +246,10 @@ int MMG2D_parsar(int argc,char *argv[],MMG5_pMesh mesh,MMG5_pSol met,MMG5_pSol s
                                      atoi(argv[i])) )
             return 0;
         }
+        else if ( !strcmp(argv[i],"-isotropic") ) {
+          if ( !MMG2D_Set_iparameter(mesh,met,MMG2D_IPARAM_isotropic,1) )
+            return 0;
+        }
         else {
           MMG2D_usage(argv[0]);
           return 0;
@@ -278,6 +293,10 @@ int MMG2D_parsar(int argc,char *argv[],MMG5_pMesh mesh,MMG5_pSol met,MMG5_pSol s
             }
           }
         }
+        else if ( !strcmp(argv[i],"-limit_angle") && ++i < argc ) {
+          if ( !MMG2D_Set_dparameter(mesh,met,MMG2D_DPARAM_limit_angle,atof(argv[i])) )
+            return 0;
+        }
         break;
       case 'm':  /* memory */
         if ( !strcmp(argv[i],"-met") ) {

src/mmg2d/libmmg2df.c

@@ -51,7 +51,7 @@ FORTRAN_NAME(MMG2D_MMG2DLIB,mmg2d_mmg2dlib,(MMG5_pMesh *mesh,MMG5_pSol *met
                                             ,int* retval),(mesh,met
                                                            ,retval)){
 
-  *retval = MMG2D_mmg2dlib(*mesh,*met);
+  *retval = MMG2D_mmg2dlib(*mesh,*met,NULL);
 
   return;
 }

src/mmg2d/mmg2d.c

@@ -275,6 +275,7 @@ int MMG2D_defaultOption(MMG5_pMesh mesh,MMG5_pSol met,MMG5_pSol sol) {
 int main(int argc,char *argv[]) {
   MMG5_pMesh    mesh;
   MMG5_pSol     sol,met,disp,ls;
+  double*       velocity;
   int           ier,ierSave,fmtin,fmtout;
   char          stim[32],*ptr;
 
@@ -469,7 +470,7 @@ int main(int argc,char *argv[]) {
       MMG2D_RETURN_AND_FREE(mesh,met,ls,disp,MMG5_STRONGFAILURE);
     }
 
-    ier = MMG2D_mmg2dlib(mesh,met);
+    ier = MMG2D_mmg2dlib(mesh,met,velocity);
   }
 
   if ( ier != MMG5_STRONGFAILURE ) {