From e17584d83b5a229d4c762c13a2212d3a882e734f Mon Sep 17 00:00:00 2001 From: Fabien Salmon Date: Fri, 12 Sep 2025 14:38:04 +0200 Subject: [PATCH 1/2] Modifications dedicated to Nextsim remeshing --- src/common/libmmgtypes.h | 3 + src/common/scalem.c | 2 +- src/mmg2d/API_functions_2d.c | 24 +- src/mmg2d/API_functionsf_2d.c | 85 ++----- src/mmg2d/delone_2d.c | 26 +- src/mmg2d/libmmg2d.c | 13 +- src/mmg2d/libmmg2d.h | 54 +++-- src/mmg2d/libmmg2d_private.h | 21 +- src/mmg2d/libmmg2d_tools.c | 19 ++ src/mmg2d/libmmg2df.c | 2 +- src/mmg2d/locate_2d.c | 4 +- src/mmg2d/mmg2d.c | 3 +- src/mmg2d/mmg2d1.c | 249 ++++++++++++------- src/mmg2d/variadic_2d.c | 434 ++++++++++++++++++++++++++++++++++ 14 files changed, 744 insertions(+), 195 deletions(-) diff --git a/src/common/libmmgtypes.h b/src/common/libmmgtypes.h index fe75d6193..5dfcdb121 100644 --- a/src/common/libmmgtypes.h +++ b/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 */ diff --git a/src/common/scalem.c b/src/common/scalem.c index 33155d91b..17da1b2e8 100644 --- a/src/common/scalem.c +++ b/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; diff --git a/src/mmg2d/API_functions_2d.c b/src/mmg2d/API_functions_2d.c index 3205094b3..44574688f 100644 --- a/src/mmg2d/API_functions_2d.c +++ b/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; diff --git a/src/mmg2d/API_functionsf_2d.c b/src/mmg2d/API_functionsf_2d.c index 930c6e156..9f22558b3 100644 --- a/src/mmg2d/API_functionsf_2d.c +++ b/src/mmg2d/API_functionsf_2d.c @@ -42,23 +42,13 @@ #include "libmmg2d_private.h" /** - * See \ref MMG2D_Init_mesh function in common/libmmgcommon_private.h file. + * See \ref MMG2D_Init_mesh function in mmg2d/libmmg2d.h file. */ -FORTRAN_VARIADIC ( MMG2D_INIT_MESH, mmg2d_init_mesh, - (const int starter, ... ), - va_list argptr; - int ier; - - va_start(argptr, starter); - - ier = MMG2D_Init_mesh_var(argptr); - - va_end(argptr); - - if ( !ier ) exit(EXIT_FAILURE); - - return; - ) +FORTRAN_NAME(MMG2D_INIT_MESH, mmg2d_init_mesh,(void **arglist, int* retval), + (arglist, retval)) { + *retval = MMG2D_Init_mesh_fortran_var(arglist); + return; +} /** * See \ref MMG2D_Init_fileNames function in mmg2d/libmmg2d.h file. @@ -732,63 +722,32 @@ FORTRAN_NAME(MMG2D_FREE_ALLSOLS,mmg2d_free_allsols, return; } - /** - * See \ref MMG2D_Free_all function in \ref mmg2d/libmmg2d.h file. + * See \ref MMG2D_Init_mesh function in mmg2d/libmmg2d.h file. */ -FORTRAN_VARIADIC(MMG2D_FREE_ALL,mmg2d_free_all, - (const int starter,...), - va_list argptr; - int ier; - - va_start(argptr, starter); - - ier = MMG2D_Free_all_var(argptr); - - va_end(argptr); - - if ( !ier ) exit(EXIT_FAILURE); - - return; - ) +FORTRAN_NAME(MMG2D_FREE_ALL, mmg2d_free_all,(void **arglist, int* retval), + (arglist,retval)) { + *retval = MMG2D_Free_all_fortran_var(arglist); + return; +} /** * See \ref MMG2D_Free_structures function in \ref mmg2d/libmmg2d.h file. */ -FORTRAN_VARIADIC(MMG2D_FREE_STRUCTURES,mmg2d_free_structures, - (const int starter,...), - va_list argptr; - int ier; - - va_start(argptr, starter); - - ier = MMG2D_Free_structures_var(argptr); - - va_end(argptr); - - if ( !ier ) exit(EXIT_FAILURE); - - return; - ) +FORTRAN_NAME(MMG2D_FREE_STRUCTURES, mmg2d_free_structures,(void **arglist, int* retval), + (arglist,retval)) { + *retval = MMG2D_Free_structures_fortran_var(arglist); + return; +} /** * See \ref MMG2D_Free_names function in \ref mmg2d/libmmg2d.h file. */ -FORTRAN_VARIADIC(MMG2D_FREE_NAMES,mmg2d_free_names, - (const int starter,...), - va_list argptr; - int ier; - - va_start(argptr, starter); - - ier = MMG2D_Free_names_var(argptr); - - va_end(argptr); - - if ( !ier ) exit(EXIT_FAILURE); - - return; - ) +FORTRAN_NAME(MMG2D_FREE_NAMES, mmg2d_free_names,(void **arglist, int* retval), + (arglist,retval)) { + *retval = MMG2D_Free_names_fortran_var(arglist); + return; +} /** * See \ref MMG2D_loadMesh function in \ref mmg2d/libmmg2d.h file. diff --git a/src/mmg2d/delone_2d.c b/src/mmg2d/delone_2d.c index 45e258800..c70bc56d7 100644 --- a/src/mmg2d/delone_2d.c +++ b/src/mmg2d/delone_2d.c @@ -260,18 +260,24 @@ int MMG2D_delone(MMG5_pMesh mesh,MMG5_pSol sol,MMG5_int ip,MMG5_int *list,int il MMG5_Hash hedg; static int8_t mmgWarn0=0,mmgWarn1=0; - /* Reset tagdel field */ - for (k=1; kpoint[k].tagdel = 0; + /* Initialize tagdel field */ + for (k=0; ktria[iel]; + for (i=0; i<3; i++) { + ppt = &mesh->point[ pt1->v[i] ]; + ppt->tagdel = 0; + } + } /* Triangles in the cavity are those s.t. pt->base == base */ base = mesh->base; /* Count the number of external faces in the cavity, and tag the corresponding vertices */ size = 0; for (k=0; ktria[old]; - iadr = (old-1)*3 + 1; + iel = list[k]; + pt1 = &mesh->tria[iel]; + iadr = (iel-1)*3 + 1; adja = &mesh->adja[iadr]; nei[0] = adja[0]/3 ; nei[1] = adja[1]/3 ; @@ -292,8 +298,8 @@ int MMG2D_delone(MMG5_pMesh mesh,MMG5_pSol sol,MMG5_int ip,MMG5_int *list,int il /* Check for an isolated vertex (the cavity should ne contain any internal vertex) */ alert = 0; for (k=0; ktria[old]; + iel = list[k]; + pt1 = &mesh->tria[iel]; for (i=0; i<3; i++) { ppt = &mesh->point[ pt1->v[i] ]; if ( !ppt->tagdel ) { @@ -303,8 +309,8 @@ int MMG2D_delone(MMG5_pMesh mesh,MMG5_pSol sol,MMG5_int ip,MMG5_int *list,int il } /* Reset tagdel field */ for (k=0; ktria[old]; + iel = list[k]; + pt1 = &mesh->tria[iel]; for (i=0; i<3; i++) { ppt = &mesh->point[ pt1->v[i] ]; ppt->tagdel = 0; diff --git a/src/mmg2d/libmmg2d.c b/src/mmg2d/libmmg2d.c index 892b95297..2df7224f8 100644 --- a/src/mmg2d/libmmg2d.c +++ b/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); } diff --git a/src/mmg2d/libmmg2d.h b/src/mmg2d/libmmg2d.h index afb57ef79..45775cdb9 100644 --- a/src/mmg2d/libmmg2d.h +++ b/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 SUBROUTINE MMG2D_INIT_MESH(arglist,retval)\n + * > MMG5_DATA_PTR_T,DIMENSION(*),INTENT(IN) :: arglist\n + * > INTEGER, INTENT(OUT) :: retval\n + * > END SUBROUTINE\n */ LIBMMG2D_EXPORT int MMG2D_Init_mesh(const int starter,...); @@ -1619,10 +1629,14 @@ LIBMMG2D_EXPORT int MMG2D_Set_lsBaseReference(MMG5_pMesh mesh, MMG5_pSol sol,MM * * \return 0 on failure, 1 on success * - * \remark we pass the structures by reference in order to have argument - * compatibility between the library call from a Fortran code and a C code. + * \remark Fortran users should provide a MMG5_DATA_PTR_T array, where every pointer to + * a MMG structure should be passed by reference (using Fortran LOC function). * - * \remark no Fortran interface to allow variadic args. + * \remark Fortran interface: + * > SUBROUTINE MMG2D_FREE_ALL(arglist,retval)\n + * > MMG5_DATA_PTR_T,DIMENSION(*),INTENT(IN) :: arglist\n + * > INTEGER, INTENT(OUT) :: retval\n + * > END SUBROUTINE\n * */ LIBMMG2D_EXPORT int MMG2D_Free_all(const int starter,...); @@ -1652,12 +1666,15 @@ LIBMMG2D_EXPORT int MMG2D_Free_all(const int starter,...); * * \return 0 on failure, 1 on success * - * \remark we pass the structures by reference in order to have argument - * compatibility between the library call from a Fortran code and a C code. - * - * \remark No fortran interface to allow variadic arguments. + * \remark Fortran users should provide a MMG5_DATA_PTR_T array, where every + * pointer to a MMG structure should be passed by reference + * (using Fortran LOC function). * - * \remark no Fortran interface to allow variadic args. + * \remark Fortran interface: + * > SUBROUTINE MMG2D_FREE_STRUCTURES(arglist,retval)\n + * > MMG5_DATA_PTR_T,DIMENSION(*),INTENT(IN) :: arglist\n + * > INTEGER, INTENT(OUT) :: retval\n + * > END SUBROUTINE\n * */ LIBMMG2D_EXPORT int MMG2D_Free_structures(const int starter,...); @@ -1687,13 +1704,16 @@ LIBMMG2D_EXPORT int MMG2D_Free_all(const int starter,...); * * \return 0 on failure, 1 otherwise * - * \remark we pass the structures by reference in order to have argument - * compatibility between the library call from a Fortran code and a C code. - * - * \remark No fortran interface to allow variadic arguments. - * - * \remark no Fortran interface to allow variadic args. + * \remark Fortran users should provide a MMG5_DATA_PTR_T array, where every + * pointer to a MMG structure should be passed by reference + * (using Fortran LOC function). * + * \remark Fortran interface: + * > SUBROUTINE MMG2D_FREE_NAMES(arglist,retval)\n + * > MMG5_DATA_PTR_T,DIMENSION(*),INTENT(IN) :: arglist\n + * > INTEGER, INTENT(OUT) :: retval\n + * > END SUBROUTINE\n + * */ LIBMMG2D_EXPORT int MMG2D_Free_names(const int starter,...); @@ -2278,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. diff --git a/src/mmg2d/libmmg2d_private.h b/src/mmg2d/libmmg2d_private.h index a700cf5f5..cab514541 100644 --- a/src/mmg2d/libmmg2d_private.h +++ b/src/mmg2d/libmmg2d_private.h @@ -200,9 +200,13 @@ int MMG2D_outqua(MMG5_pMesh ,MMG5_pSol ); int MMG2D_mmg2d1(MMG5_pMesh ,MMG5_pSol ); int MMG2D_Init_mesh_var( va_list argptr ); +int MMG2D_Init_mesh_fortran_var(void ** arglist); int MMG2D_Free_all_var( va_list argptr ); +int MMG2D_Free_all_fortran_var(void ** arglist); int MMG2D_Free_structures_var( va_list argptr ); +int MMG2D_Free_structures_fortran_var( void ** arglist ); int MMG2D_Free_names_var( va_list argptr ); +int MMG2D_Free_names_fortran_var( void ** arglist ); int MMG2D_mmg2d2(MMG5_pMesh , MMG5_pSol); int MMG2D_mmg2d6(MMG5_pMesh ,MMG5_pSol,MMG5_pSol ); @@ -250,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 ); @@ -261,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 ); -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 ); @@ -283,9 +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 ); -int MMG2D_adpcol(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 ); diff --git a/src/mmg2d/libmmg2d_tools.c b/src/mmg2d/libmmg2d_tools.c index 0ebf7e9ab..c48b53af4 100644 --- a/src/mmg2d/libmmg2d_tools.c +++ b/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") ) { diff --git a/src/mmg2d/libmmg2df.c b/src/mmg2d/libmmg2df.c index 955ccf4c7..953eb7325 100644 --- a/src/mmg2d/libmmg2df.c +++ b/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; } diff --git a/src/mmg2d/locate_2d.c b/src/mmg2d/locate_2d.c index 3b95a91ed..0a734b1e2 100644 --- a/src/mmg2d/locate_2d.c +++ b/src/mmg2d/locate_2d.c @@ -102,11 +102,11 @@ int MMG2D_cutEdge(MMG5_pMesh mesh,MMG5_pTria pt,MMG5_pPoint ppa,MMG5_pPoint ppb) /* Check whether ppa or ppb is a vertex of pt */ for (i=0; i<3; i++) { - if ( fabs(la[i]-1.0) < 1.0e-12 ) { + if (( fabs(la[i]-1.0) < 1.0e-12 ) && ( fabs(la[MMG5_inxt2[i]]) < 1.0e-12 ) && ( fabs(la[MMG5_iprv2[i]]) < 1.0e-12 )) { if ( lb[i] < 0.0 ) return i+1; else return 0; } - if ( fabs(lb[i]-1.0) < 1.0e-12) { + if (( fabs(lb[i]-1.0) < 1.0e-12 ) && ( fabs(lb[MMG5_inxt2[i]]) < 1.0e-12 ) && ( fabs(lb[MMG5_iprv2[i]]) < 1.0e-12 )) { if ( la[i] < 0.0 ) return i+1; else return 0; } diff --git a/src/mmg2d/mmg2d.c b/src/mmg2d/mmg2d.c index 65f8c7d89..b4b5e1d32 100644 --- a/src/mmg2d/mmg2d.c +++ b/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 ) { diff --git a/src/mmg2d/mmg2d1.c b/src/mmg2d/mmg2d1.c index 3f1c7dd38..a7de6734b 100644 --- a/src/mmg2d/mmg2d1.c +++ b/src/mmg2d/mmg2d1.c @@ -37,7 +37,7 @@ based on patterns, collapses and swaps. typchk = 1 -> adaptation based on edge lengths typchk = 2 -> adaptation based on lengths calculated in metric met */ -int MMG2D_anatri(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) { +int MMG2D_anatri(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk,double* velocity) { int it,maxit; MMG5_int ns,nc,nsw,nns,nnc,nnsw; @@ -50,13 +50,13 @@ int MMG2D_anatri(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) { if ( typchk == 2 && it == 0 ) { // #warning Luca: check consistency with 3D } - if ( !mesh->info.noinsert ) { + if ( !mesh->info.noinsert && mesh->info.limit_angle >= 4.*atan(1.)) { /* Memory free */ MMG5_DEL_MEM(mesh,mesh->adja); mesh->adja = 0; /* Split long edges according to patterns */ - ns = MMG2D_anaelt(mesh,met,typchk); + ns = MMG2D_anaelt(mesh,met,typchk,velocity); if ( ns < 0 ) { fprintf(stderr," ## Unable to complete surface mesh. Exit program.\n"); return 0; @@ -69,7 +69,7 @@ int MMG2D_anatri(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) { } /* Collapse short edges */ - nc = MMG2D_colelt(mesh,met,typchk); + nc = MMG2D_colelt(mesh,met,typchk,MMG2D_LSHRT,velocity); if ( nc < 0 ) { fprintf(stderr," ## Unable to collapse mesh. Exiting.\n"); return 0; @@ -82,7 +82,7 @@ int MMG2D_anatri(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) { } /* Swap edges */ if ( !mesh->info.noswap ) { - nsw = MMG2D_swpmsh(mesh,met,typchk); + nsw = MMG2D_swpmsh(mesh,met,typchk,velocity); if ( nsw < 0 ) { fprintf(stderr," ## Unable to improve mesh. Exiting.\n"); return 0; @@ -107,8 +107,134 @@ int MMG2D_anatri(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) { return 1; } +// Minimum angle function +double MMG2D_minangle(MMG5_pMesh mesh, int k, double* velocity) { + + MMG5_pTria pt = &mesh->tria[k]; + MMG5_pPoint pp1, pp2, pp3; + int* adja = &mesh->adja[3*(k-1)+1]; + int p1, p2, p3; + double length[6]; + double c1 = 1.35*mesh->info.limit_angle; + if (c1 > 1.) c1 = 1.; + double coef = 1.; + double factor_max = mesh->info.max[2]; // 1.75; + double factor_min = mesh->info.min[2]; // 1.33; + + p1 = pt->v[0]; + p2 = pt->v[1]; + p3 = pt->v[2]; + + pp1 = &mesh->point[p1]; + pp2 = &mesh->point[p2]; + pp3 = &mesh->point[p3]; + + double square_root_area = sqrt(0.5*fabs((pp2->c[0]-pp1->c[0])*(pp3->c[1]-pp1->c[1]) - (pp2->c[1]-pp1->c[1])*(pp3->c[0]-pp1->c[0]))); + double mean_velocity = 0.; + + // Velocity is only considered for already existing points + if (velocity != NULL && p1 <= mesh->mark && p2 <= mesh->mark && p3 <= mesh->mark && pp1->tmp == -1 && pp2->tmp == -1 && pp3->tmp == -1) + mean_velocity = ( sqrt(pow(velocity[p1-1],2) + pow(velocity[p1-1+mesh->mark],2)) + sqrt(pow(velocity[p2-1],2) + pow(velocity[p2-1+mesh->mark],2)) + sqrt(pow(velocity[p3-1],2) + pow(velocity[p3-1+mesh->mark],2)) )/3.; + + double dt, x1, x2, x3, y1, y2, y3; + if (mean_velocity < 1.e12) { + dt = 0.; + x1 = pp1->c[0]; + x2 = pp2->c[0]; + x3 = pp3->c[0]; + y1 = pp1->c[1]; + y2 = pp2->c[1]; + y3 = pp3->c[1]; + } + else { + dt = 0.25*square_root_area / mean_velocity; + x1 = pp1->c[0]+velocity[p1-1]*dt; + x2 = pp2->c[0]+velocity[p2-1]*dt; + x3 = pp3->c[0]+velocity[p3-1]*dt; + y1 = pp1->c[1]+velocity[p1-1+mesh->mark]*dt; + y2 = pp2->c[1]+velocity[p2-1+mesh->mark]*dt; + y3 = pp3->c[1]+velocity[p3-1+mesh->mark]*dt; + } + + length[0] = pow(x1-x2,2.) + pow(y1-y2,2.); + length[1] = pow(x2-x3,2.) + pow(y2-y3,2.); + length[2] = pow(x3-x1,2.) + pow(y3-y1,2.); + + length[3] = pow(pp1->c[0] - pp2->c[0],2.) + pow(pp1->c[1] - pp2->c[1],2.); + length[4] = pow(pp2->c[0] - pp3->c[0],2.) + pow(pp2->c[1] - pp3->c[1],2.); + length[5] = pow(pp3->c[0] - pp1->c[0],2.) + pow(pp3->c[1] - pp1->c[1],2.); + + double min = length[0]; + double max = length[0]; + int min_index = 0; + + for (int i = 1; i < 3; i++) { + if (length[i] < min) { + min = length[i]; + min_index = i; + } + if (length[i] > max) max = length[i]; + } + + if (sqrt(min) < mesh->info.hmin/factor_min) return 0.; + if (sqrt(max) > factor_max*mesh->info.hmax) return 0.; + + double cosA = ( length[(min_index+1)%3] + length[(min_index+2)%3] - length[min_index] ) / (2*pow(length[(min_index+1)%3],0.5)*pow(length[(min_index+2)%3],0.5)); + + min = length[3]; + max = length[3]; + min_index = 0; + for (int i = 1; i < 3; i++) { + if (length[i+3] < min) { + min = length[i+3]; + min_index = i; + } + if (length[i+3] > max) max = length[i+3]; + } + + if (sqrt(min) < mesh->info.hmin/factor_min) return 0.; + if (sqrt(max) > factor_max*mesh->info.hmax) return 0.; + + // If the triangle is on a boundary, the tolerance is multiplied by two + if (mesh->info.bdy_adaptation) { + if ((pp1->tag & MG_BDY) || (pp2->tag & MG_BDY) || (pp3->tag & MG_BDY)) { + coef = c1; + } + else { + for(int i = 0; i < 3; i++ ) { + int k1 = adja[i]/3; + + if (!k1) continue; + + MMG5_pTria pt_neigh = &mesh->tria[k1]; + + for (int j = 0; j < 3; j++) { + MMG5_pPoint pp = &mesh->point[pt_neigh->v[j]]; + if (pp->tag & MG_BDY) { + coef = c1; + break; + } + } + + if (coef < 0.999) break; + } + } + } + + double cosA2 = ( length[(min_index+1)%3+3] + length[(min_index+2)%3+3] - length[min_index+3] ) / (2*pow(length[(min_index+1)%3+3],0.5)*pow(length[(min_index+2)%3+3],0.5)); + if (cosA2 > cosA) cosA = cosA2; + + if (cosA <= -1.) { + return 3.14; + } + else if (cosA >= 1) { + return 0.; + } + else return coef*acos(cosA); +} + /* Travel triangles and split long edges according to patterns */ -MMG5_int MMG2D_anaelt(MMG5_pMesh mesh,MMG5_pSol met,int typchk) { +MMG5_int MMG2D_anaelt(MMG5_pMesh mesh,MMG5_pSol met,int typchk,double* velocity) { MMG5_pTria pt; MMG5_pPoint ppt,p1,p2; MMG5_Hash hash; @@ -130,6 +256,8 @@ MMG5_int MMG2D_anaelt(MMG5_pMesh mesh,MMG5_pSol met,int typchk) { if ( !MG_EOK(pt) || (pt->ref < 0) ) continue; if ( MG_SIN(pt->tag[0]) || MG_SIN(pt->tag[1]) || MG_SIN(pt->tag[2]) ) continue; + if (mesh->info.limit_angle <= 4*atan(1.) && MMG2D_minangle(mesh,k,velocity) > mesh->info.limit_angle) continue; + /* Check if pt should be cut */ pt->flag = 0; @@ -434,10 +562,10 @@ int MMG2D_dichoto(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int *vx) { } /* Travel triangles and collapse short edges */ -MMG5_int MMG2D_colelt(MMG5_pMesh mesh,MMG5_pSol met,int typchk) { +MMG5_int MMG2D_colelt(MMG5_pMesh mesh,MMG5_pSol met,int typchk, double lmax,double* velocity) { MMG5_pTria pt; MMG5_pPoint p1,p2; - double ux,uy,ll,hmin2; + double ll; MMG5_int k; int ilist; MMG5_int nc; @@ -445,12 +573,13 @@ MMG5_int MMG2D_colelt(MMG5_pMesh mesh,MMG5_pSol met,int typchk) { MMG5_int list[MMG5_TRIA_LMAX+2]; nc = 0; - hmin2 = mesh->info.hmin * mesh->info.hmin; for (k=1; k<=mesh->nt; k++) { pt = &mesh->tria[k]; if ( !MG_EOK(pt) || pt->ref < 0 ) continue; + if (mesh->info.limit_angle <= 4*atan(1.) && MMG2D_minangle(mesh,k,velocity) > mesh->info.limit_angle) continue; + /* Travel 3 edges of the triangle and decide whether to collapse p1->p2, based on length criterion */ pt->flag = 0; // was here before, but probably serves for nothing @@ -476,14 +605,16 @@ MMG5_int MMG2D_colelt(MMG5_pMesh mesh,MMG5_pSol met,int typchk) { /* Check length */ if ( typchk == 1 ) { + double ux, uy, hmin2; ux = p2->c[0] - p1->c[0]; uy = p2->c[1] - p1->c[1]; ll = ux*ux + uy*uy; + hmin2 = mesh->info.hmin * mesh->info.hmin; if ( ll > hmin2 ) continue; } else { ll = MMG2D_lencurv(mesh,met,pt->v[i1],pt->v[i2]); - if ( ll > MMG2D_LSHRT ) continue; + if ( ll > lmax ) continue; } /* Check whether the geometry is preserved */ @@ -511,7 +642,7 @@ MMG5_int MMG2D_colelt(MMG5_pMesh mesh,MMG5_pSol met,int typchk) { } /* Travel triangles and swap edges to improve quality */ -MMG5_int MMG2D_swpmsh(MMG5_pMesh mesh,MMG5_pSol met,int typchk) { +MMG5_int MMG2D_swpmsh(MMG5_pMesh mesh,MMG5_pSol met,int typchk,double* velocity) { MMG5_pTria pt; int it,maxit; MMG5_int ns,nns; @@ -528,6 +659,8 @@ MMG5_int MMG2D_swpmsh(MMG5_pMesh mesh,MMG5_pSol met,int typchk) { pt = &mesh->tria[k]; if ( !MG_EOK(pt) || pt->ref < 0 ) continue; + if (mesh->info.limit_angle <= 4*atan(1.) && MMG2D_minangle(mesh,k,velocity) > mesh->info.limit_angle) continue; + for (i=0; i<3; i++) { if ( MG_SIN(pt->tag[i]) || MG_EDG(pt->tag[i]) ) continue; else if ( MMG2D_chkswp(mesh,met,k,i,typchk) ) { @@ -547,7 +680,7 @@ MMG5_int MMG2D_swpmsh(MMG5_pMesh mesh,MMG5_pSol met,int typchk) { /* Mesh adaptation routine for the final stage of the algorithm: intertwine splitting based on patterns, collapses, swaps and vertex relocations.*/ -int MMG2D_adptri(MMG5_pMesh mesh,MMG5_pSol met) { +int MMG2D_adptri(MMG5_pMesh mesh,MMG5_pSol met,double* velocity) { int maxit,it; MMG5_int nns,ns,nnc,nc,nnsw,nsw,nnm,nm; @@ -557,16 +690,16 @@ int MMG2D_adptri(MMG5_pMesh mesh,MMG5_pSol met) { do { if ( !mesh->info.noinsert ) { - ns = MMG2D_adpspl(mesh,met); + ns = MMG2D_adpspl(mesh,met,velocity); if ( ns < 0 ) { fprintf(stderr," ## Problem in function adpspl." " Unable to complete mesh. Exit program.\n"); return 0; } - nc = MMG2D_adpcol(mesh,met); + nc = MMG2D_colelt(mesh,met,2,MMG2D_LOPTS,velocity); if ( nc < 0 ) { - fprintf(stderr," ## Problem in function adpcol." + fprintf(stderr," ## Problem in function colelt." " Unable to complete mesh. Exit program.\n"); return 0; } @@ -577,7 +710,7 @@ int MMG2D_adptri(MMG5_pMesh mesh,MMG5_pSol met) { } if ( !mesh->info.noswap ) { - nsw = MMG2D_swpmsh(mesh,met,2); + nsw = MMG2D_swpmsh(mesh,met,2,velocity); if ( nsw < 0 ) { fprintf(stderr," ## Problem in function swpmsh." " Unable to complete mesh. Exit program.\n"); @@ -588,7 +721,7 @@ int MMG2D_adptri(MMG5_pMesh mesh,MMG5_pSol met) { nsw = 0; if ( !mesh->info.nomove ) { - nm = MMG2D_movtri(mesh,met,1,0); + nm = MMG2D_movtri(mesh,met,1,0,velocity); if ( nm < 0 ) { fprintf(stderr," ## Problem in function movtri. " "Unable to complete mesh. Exit program.\n"); @@ -612,7 +745,7 @@ int MMG2D_adptri(MMG5_pMesh mesh,MMG5_pSol met) { /* Last iterations of vertex relocation only */ if ( !mesh->info.nomove ) { - nm = MMG2D_movtri(mesh,met,5,1); + nm = MMG2D_movtri(mesh,met,5,1,velocity); if ( nm < 0 ) { fprintf(stderr," ## Problem in function movtri. Unable to complete mesh." " Exit program.\n"); @@ -637,7 +770,7 @@ int MMG2D_adptri(MMG5_pMesh mesh,MMG5_pSol met) { * edges are only splitted on a one-by-one basis * */ -MMG5_int MMG2D_adpspl(MMG5_pMesh mesh,MMG5_pSol met) { +MMG5_int MMG2D_adpspl(MMG5_pMesh mesh,MMG5_pSol met,double* velocity) { MMG5_pTria pt; double lmax,len; MMG5_int ip,ns,k,nt; @@ -652,6 +785,8 @@ MMG5_int MMG2D_adpspl(MMG5_pMesh mesh,MMG5_pSol met) { pt = &mesh->tria[k]; if ( !MG_EOK(pt) || pt->ref < 0 ) continue; + if (mesh->info.limit_angle <= 4*atan(1.) && MMG2D_minangle(mesh,k,velocity) > mesh->info.limit_angle) continue; + imax = -1; lmax = -1.0; for (i=0; i<3; i++) { @@ -692,67 +827,8 @@ MMG5_int MMG2D_adpspl(MMG5_pMesh mesh,MMG5_pSol met) { return ns; } -/* Analysis and collapse routine for edges in the final step of the algorithm */ -int MMG2D_adpcol(MMG5_pMesh mesh,MMG5_pSol met) { - MMG5_pTria pt; - MMG5_pPoint p1,p2; - double len; - MMG5_int k,nc; - int ilist; - int8_t i,i1,i2,open; - MMG5_int list[MMG5_TRIA_LMAX+2]; - - nc = 0; - for (k=1; k<=mesh->nt; k++) { - pt = &mesh->tria[k]; - if ( !MG_EOK(pt) || pt->ref < 0 ) continue; - - /* Check edge length, and attempt collapse */ - pt->flag = 0; - for (i=0; i<3; i++) { - if ( MG_SIN(pt->tag[i]) ) continue; - - open = ( mesh->adja[3*(k-1)+1+i] == 0 ) ? 1 : 0; - - i1 = MMG5_inxt2[i]; - i2 = MMG5_iprv2[i]; - p1 = &mesh->point[pt->v[i1]]; - p2 = &mesh->point[pt->v[i2]]; - - if ( MG_SIN(p1->tag) || p1->tag & MG_NOM ) continue; - else if ( p1->tag & MG_GEO ) { - if ( ! (p2->tag & MG_GEO) || !(pt->tag[i] & MG_GEO) ) continue; - } - else if ( p1->tag & MG_REF ) { - if ( ! (p2->tag & MG_GEO || p2->tag & MG_REF) || !(pt->tag[i] & MG_REF) ) continue; - } - - len = MMG2D_lencurv(mesh,met,pt->v[i1],pt->v[i2]); - - if ( len > MMG2D_LOPTS ) continue; - - ilist = MMG2D_chkcol(mesh,met,k,i,list,2); - - if ( ilist > 3 || ( ilist==3 && open ) ) { - nc += MMG2D_colver(mesh,ilist,list); - break; - } - else if ( ilist == 3 ) { - nc += MMG2D_colver3(mesh,list); - break; - } - else if ( ilist == 2 ) { - nc += MMG2D_colver2(mesh,list); - break; - } - } - } - - return nc; -} - /* Analyze points to relocate them according to a quality criterion */ -MMG5_int MMG2D_movtri(MMG5_pMesh mesh,MMG5_pSol met,int maxit,int8_t improve) { +MMG5_int MMG2D_movtri(MMG5_pMesh mesh,MMG5_pSol met,int maxit,int8_t improve,double* velocity) { MMG5_pTria pt; MMG5_pPoint p0; MMG5_int nnm,nm,ns,k; @@ -773,6 +849,8 @@ MMG5_int MMG2D_movtri(MMG5_pMesh mesh,MMG5_pSol met,int maxit,int8_t improve) { pt = &mesh->tria[k]; if ( !MG_EOK(pt) || pt->ref < 0 ) continue; + if (mesh->info.limit_angle <= 4*atan(1.) && MMG2D_minangle(mesh,k,velocity) > mesh->info.limit_angle) continue; + for (i=0; i<3; i++) { p0 = &mesh->point[pt->v[i]]; if ( p0->flag == base || MG_SIN(p0->tag) || p0->tag & MG_NOM ) continue; @@ -785,7 +863,7 @@ MMG5_int MMG2D_movtri(MMG5_pMesh mesh,MMG5_pSol met,int maxit,int8_t improve) { if ( ier ) ns++; } else { - if ( met->size == 3 && met->m ) + if ( met->size == 3 && met->m && !mesh->info.isotropic) ier = MMG2D_movintpt_ani(mesh,met,ilist,list,improve); else ier = MMG2D_movintpt(mesh,met,ilist,list,improve); @@ -816,13 +894,16 @@ MMG5_int MMG2D_movtri(MMG5_pMesh mesh,MMG5_pSol met,int maxit,int8_t improve) { * Mesh adaptation -- new version of mmg2d1.c * **/ -int MMG2D_mmg2d1n(MMG5_pMesh mesh,MMG5_pSol met) { +int MMG2D_mmg2d1n(MMG5_pMesh mesh,MMG5_pSol met,double* velocity) { + + mesh->mark = mesh->np; + for (int i = 1; i <= mesh->np; i++) mesh->point[i].tmp = -1; /* Stage 1: creation of a geometric mesh */ if ( abs(mesh->info.imprim) > 4 || mesh->info.ddebug ) fprintf(stdout," ** GEOMETRIC MESH\n"); - if ( !MMG2D_anatri(mesh,met,1) ) { + if ( !MMG2D_anatri(mesh,met,1,velocity) ) { fprintf(stderr," ## Unable to split mesh-> Exiting.\n"); return 0; } @@ -871,7 +952,7 @@ int MMG2D_mmg2d1n(MMG5_pMesh mesh,MMG5_pSol met) { return 1; } - if ( !MMG2D_anatri(mesh,met,2) ) { + if ( !MMG2D_anatri(mesh,met,2,velocity) ) { fprintf(stderr," ## Unable to proceed adaptation. Exit program.\n"); return 0; } @@ -885,7 +966,7 @@ int MMG2D_mmg2d1n(MMG5_pMesh mesh,MMG5_pSol met) { /* Stage 3: fine mesh improvements */ - if ( !MMG2D_adptri(mesh,met) ) { + if ( !MMG2D_adptri(mesh,met,velocity) ) { fprintf(stderr," ## Unable to make fine improvements. Exit program.\n"); return 0; } diff --git a/src/mmg2d/variadic_2d.c b/src/mmg2d/variadic_2d.c index 8053ec901..8fac72d3f 100644 --- a/src/mmg2d/variadic_2d.c +++ b/src/mmg2d/variadic_2d.c @@ -227,6 +227,90 @@ int MMG2D_Init_mesh_var( va_list argptr ) { return 1; } +/** + * \param arglist list of the mmg structures that must be initialized. Each + * structure must follow one of the \a MMG5_ARG* preprocessor variable that allow to identify + * it. + * + * \a arglist contains at least a pointer to a \a MMG5_pMesh structure + * (that will contain the mesh and identified by the MMG5_ARG_ppMesh keyword) + * + * To call the \a MMG2D_mmg2dlib function, you must also provide + * a pointer to a \a MMG5_pSol structure (that will contain the ouput + * metric (and the input one, if provided) and identified by the MMG5_ARG_ppMet + * keyword). + * + * To call the \a MMG2D_mmg2dls function, you must also provide a pointer + * toward a \a MMG5_pSol structure (that will contain the level-set function and + * identified by the MMG5_ARG_ppLs keyword). + * + * To call the \a MMG2D_mmg2dmov library, you must also provide a + * pointer to a \a MMG5_pSol structure storing the displacement (and + * identified by the MMG5_ARG_ppDisp keyword). + * + * \return 0 if fail, 1 otherwise + * + * Internal function for structure allocations. + * Fortran users should provide a MMG5_DATA_PTR_T array, where every pointer to + * a MMG structure should be passed by reference (using Fortran LOC function). + */ +int MMG2D_Init_mesh_fortran_var(void **arglist) { + + MMG5_pMesh *mesh; + MMG5_pSol *sol,*disp,*ls; + int typArg,i; + int meshCount; + + meshCount = 0; + mesh = NULL; + disp = sol = ls = NULL; + i = 1; + + while ( (typArg = (intptr_t)arglist[i]) != MMG5_ARG_end ) + { + switch ( typArg ) + { + case(MMG5_ARG_ppMesh): + mesh = (MMG5_pMesh*)arglist[i+1]; + ++meshCount; + break; + case(MMG5_ARG_ppMet): + sol = (MMG5_pSol*)arglist[i+1]; + break; + case(MMG5_ARG_ppLs): + ls = (MMG5_pSol*)arglist[i+1]; + break; + case(MMG5_ARG_ppDisp): + disp = (MMG5_pSol*)arglist[i+1]; + break; + default: + fprintf(stderr,"\n ## Error: %s: MMG2D_Init_mesh:\n" + " unexpected argument type: %d\n",__func__,typArg); + fprintf(stderr," Argument type must be one of the following" + " preprocessor variable: MMG5_ARG_ppMesh, MMG5_ARG_ppMet," + " MMG5_ARG_ppLs, MMG5_ARG_ppDisp\n"); + return 0; + } + i += 2; + } + + if ( meshCount !=1 ) { + fprintf(stderr,"\n ## Error: %s: MMG2D_Init_mesh:\n" + " you need to initialize the mesh structure that" + " will contain your mesh.\n",__func__); + return 0; + } + + /* allocations */ + if ( !MMG2D_Alloc_mesh(mesh,sol,ls,disp) ) { + return 0; + } + + /* initialisations */ + MMG2D_Init_woalloc_mesh(mesh,sol,ls,disp); + return 1; +} + /** * \param argptr list of the mmg structures that must be deallocated. Each * structure must follow one of the \a MMG5_ARG preprocessor variable that allow to @@ -339,6 +423,119 @@ int MMG2D_Free_all_var(va_list argptr) return ier; } +/** + * \param argptr list of the mmg structures that must be deallocated. Each + * structure must follow one of the \a MMG5_ARG preprocessor variable that allow to + * identify it. + * + * \a argptr contains at least a pointer to a \a MMG5_pMesh structure + * (that will contain the mesh and identified by the MMG5_ARG_ppMesh keyword) + * + * To call the \a MMG2D_mmg2dlib function, you must also provide + * a pointer to a \a MMG5_pSol structure (that will contain the ouput + * metric (and the input one, if provided) and identified by the MMG5_ARG_ppMet + * keyword). + * + * To call the \a MMG2D_mmg2dls function, you must also provide a pointer + * toward a \a MMG5_pSol structure (that will contain the level-set function and + * identified by the MMG5_ARG_ppLs keyword). + * + * To call the \a MMG2D_mmg2dmov library, you must also provide a + * pointer to a \a MMG5_pSol structure storing the displacement (and + * identified by the MMG5_ARG_ppDisp keyword). + * + * \return 0 if fail, 1 if success + * + * Internal function for deallocations before return. + * Fortran users should provide a MMG5_DATA_PTR_T array, where every pointer to + * a MMG structure should be passed by reference (using Fortran LOC function). + * + */ +int MMG2D_Free_all_fortran_var(void **arglist) +{ + + MMG5_pMesh *mesh; + MMG5_pSol *sol,*disp,*ls,*sols; + int typArg; + int meshCount,metCount,lsCount,dispCount,fieldsCount; + int ier, i; + + meshCount = metCount = lsCount = dispCount = fieldsCount = 0; + disp = sol = sols = ls = NULL; + i = 1; + + while ( (typArg = (intptr_t)arglist[i]) != MMG5_ARG_end ) + { + switch ( typArg ) + { + case(MMG5_ARG_ppMesh): + mesh = (MMG5_pMesh*)arglist[i+1]; + ++meshCount; + break; + case(MMG5_ARG_ppMet): + ++metCount; + sol = (MMG5_pSol*)arglist[i+1]; + break; + case(MMG5_ARG_ppLs): + ++lsCount; + ls = (MMG5_pSol*)arglist[i+1]; + break; + case(MMG5_ARG_ppDisp): + ++dispCount; + disp = (MMG5_pSol*)arglist[i+1]; + break; + case(MMG5_ARG_ppSols): + ++fieldsCount; + sols = (MMG5_pSol*)arglist[i+1]; + break; + default: + fprintf(stderr,"\n ## Error: %s: MMG2D_Free_all:\n" + " unexpected argument type: %d\n",__func__,typArg); + fprintf(stderr," Argument type must be one of the following" + " preprocessor variable: MMG5_ARG_ppMesh or MMG5_ARG_ppMet\n"); + return 0; + } + i+=2; + } + + if ( meshCount !=1 ) { + fprintf(stderr,"\n ## Error: %s: MMG2D_Free_all:\n" + " you need to provide your mesh structure" + " to allow to free the associated memory.\n",__func__); + return 0; + } + + if ( metCount > 1 || lsCount > 1 || dispCount > 1 || fieldsCount > 1 ) { + fprintf(stdout,"\n ## Warning: %s: MMG2D_Free_all:\n" + " This function can free only one structure of each type.\n" + " Probable memory leak.\n", + __func__); + } + + ier = MMG2D_Free_structures(MMG5_ARG_start, + MMG5_ARG_ppMesh, mesh, MMG5_ARG_ppMet, sol, + MMG5_ARG_ppLs,ls ,MMG5_ARG_ppDisp, disp, + MMG5_ARG_ppSols, sols, + MMG5_ARG_end); + + if ( sol ) + MMG5_SAFE_FREE(*sol); + + if ( disp ) + MMG5_SAFE_FREE(*disp); + + if ( ls ) + MMG5_SAFE_FREE(*ls); + + if ( sols ) { + MMG5_DEL_MEM(*mesh,*sols); + } + + MMG5_SAFE_FREE(*mesh); + + return ier; +} + /** * \param argptr list of the mmg structures that must be deallocated. Each * structure must follow one of the \a MMG5_ARG* preprocessor variable that allow @@ -462,6 +659,131 @@ int MMG2D_Free_structures_var(va_list argptr) return 1; } +/** + * \param argptr list of the mmg structures that must be deallocated. Each + * structure must follow one of the \a MMG5_ARG* preprocessor variable that allow + * to identify it. \a argptr contains at least a pointer to a \a MMG5_pMesh + * structure (that will contain the mesh and identified by the MMG5_ARG_ppMesh + * keyword) and a pointer to a \a MMG5_pSol structure (that will contain the + * ouput metric (and the input one, if provided) and identified by the + * MMG5_ARG_ppMet keyword). + * + * To call the \a MMG2D_mmg2dls function, you must also provide a pointer + * toward a \a MMG5_pSol structure (that will contain the level-set function and + * identified by the MMG5_ARG_ppLs keyword). + * + * To call the \a MMG2D_mmg2dmov library, you must also provide a + * pointer to a \a MMG5_pSol structure storing the displacement (and + * identified by the MMG5_ARG_ppDisp keyword). + * + * \return 1 if success, 0 if fail + * + * Internal function for structures deallocations before return (taking a + * va_list as argument). + * + * Fortran users should provide a MMG5_DATA_PTR_T array, where every pointer to + * a MMG structure should be passed by reference (using Fortran LOC function). + * + */ +int MMG2D_Free_structures_fortran_var(void **arglist) +{ + + MMG5_pMesh *mesh; + MMG5_pSol *sol,*disp,*ls,*sols; + int typArg,i; + int meshCount; + + meshCount = 0; + mesh = NULL; + disp = sol = ls = sols = NULL; + i = 1; + + while ( (typArg = (intptr_t)arglist[i]) != MMG5_ARG_end ) + { + switch ( typArg ) + { + case(MMG5_ARG_ppMesh): + mesh = (MMG5_pMesh*)arglist[i+1]; + ++meshCount; + break; + case(MMG5_ARG_ppMet): + sol = (MMG5_pSol*)arglist[i+1]; + break; + case(MMG5_ARG_ppLs): + ls = (MMG5_pSol*)arglist[i+1]; + break; + case(MMG5_ARG_ppDisp): + disp = (MMG5_pSol*)arglist[i+1]; + break; + case(MMG5_ARG_ppSols): + sols = (MMG5_pSol*)arglist[i+1]; + break; + default: + fprintf(stderr,"\n ## Error: %s: MMG2D_Free_structures:\n" + " unexpected argument type: %d\n",__func__,typArg); + fprintf(stderr," Argument type must be one of the following" + " preprocessor variable: MMG5_ARG_ppMesh or MMG5_ARG_ppMet\n"); + return 0; + } + i+=2; + } + + if ( meshCount !=1 ) { + fprintf(stderr,"\n ## Error: %s: MMG2D_Free_structures:\n" + " you need to provide your mesh structure" + " to allow to free the associated memory.\n",__func__); + return 0; + } + + if ( !MMG2D_Free_names(MMG5_ARG_start, + MMG5_ARG_ppMesh, mesh, MMG5_ARG_ppMet, sol, + MMG5_ARG_ppLs, ls ,MMG5_ARG_ppDisp, disp, + MMG5_ARG_ppSols, sols, + MMG5_ARG_end) ) + return 0; + + /* mesh */ + assert(mesh && *mesh); + + if ( (*mesh)->edge ) + MMG5_DEL_MEM((*mesh),(*mesh)->edge); + + if ( (*mesh)->adja ) + MMG5_DEL_MEM((*mesh),(*mesh)->adja); + + if ( (*mesh)->adjq ) + MMG5_DEL_MEM((*mesh),(*mesh)->adjq); + + if ( (*mesh)->tria ) + MMG5_DEL_MEM((*mesh),(*mesh)->tria); + + if ( (*mesh)->quadra ) + MMG5_DEL_MEM((*mesh),(*mesh)->quadra); + + /* disp */ + if ( disp && (*disp) && (*disp)->m ) + MMG5_DEL_MEM((*mesh),(*disp)->m); + + /* ls */ + if ( ls && (*ls) && (*ls)->m ) + MMG5_DEL_MEM((*mesh),(*ls)->m); + + /* met */ + if ( sol && (*sol) && (*sol)->m ) + MMG5_DEL_MEM((*mesh),(*sol)->m); + + /* field */ + if ( sols && (*mesh)->nsols ) { + for ( i=0; i<(*mesh)->nsols; ++i ) { + MMG5_DEL_MEM((*mesh),(*sols)[i].m); + } + } + + MMG5_Free_structures(*mesh,NULL); + + return 1; +} + /** * \param argptr list of the mmg structures for whose we want to deallocate the * name. Each structure must follow one of the \a MMG5_ARG* preprocessor variable @@ -572,3 +894,115 @@ int MMG2D_Free_names_var(va_list argptr) return 1; } +/** + * \param argptr list of the mmg structures for whose we want to deallocate the + * name. Each structure must follow one of the \a MMG5_ARG* preprocessor variable + * that allow to identify it. \a argptr contains at least a pointer to a \a + * MMG5_pMesh structure (that will contain the mesh and identified by the + * MMG5_ARG_ppMesh keyword) and a pointer to a \a MMG5_pSol structure (that + * will contain the ouput metric (and the input one, if provided) and identified + * by the MMG5_ARG_ppMet keyword). + * + * Internal function for name deallocations before return (taking a va_list as + * argument). + * + * \return 0 if fail, 1 if success + * + * \remark Fortran users should provide a MMG5_DATA_PTR_T array, where every pointer to + * a MMG structure should be passed by reference (using Fortran LOC function). + * + */ +int MMG2D_Free_names_fortran_var(void ** arglist) +{ + + MMG5_pMesh *mesh; + MMG5_pSol psl,*sol,*disp,*ls,*sols; + int typArg,i; + int meshCount; + + meshCount = 0; + disp = sol = ls = sols = NULL; + i = 1; + + while ( (typArg = (intptr_t)arglist[i]) != MMG5_ARG_end ) + { + switch ( typArg ) + { + case(MMG5_ARG_ppMesh): + mesh = (MMG5_pMesh*)arglist[i+1]; + ++meshCount; + break; + case(MMG5_ARG_ppMet): + sol = (MMG5_pSol*)arglist[i+1]; + break; + case(MMG5_ARG_ppLs): + ls = (MMG5_pSol*)arglist[i+1]; + break; + case(MMG5_ARG_ppDisp): + disp = (MMG5_pSol*)arglist[i+1]; + break; + case(MMG5_ARG_ppSols): + sols = (MMG5_pSol*)arglist[i+1]; + break; + default: + fprintf(stderr,"\n ## Error: %s: MMG2D_Free_names:\n" + " unexpected argument type: %d\n",__func__,typArg); + fprintf(stderr," Argument type must be one of the following" + " preprocessor variable: MMG5_ARG_ppMesh or MMG5_ARG_ppMet\n"); + return 0; + } + i+=2; + } + + if ( meshCount !=1 ) { + fprintf(stderr,"\n ## Error: %s: MMG2D_Free_names:\n" + " you need to provide your mesh structure" + " to allow to free the associated memory.\n",__func__); + return 0; + } + + /* mesh & met */ + if ( sol ) { + MMG5_mmgFree_names(*mesh,*sol); + } + else { + MMG5_mmgFree_names(*mesh,NULL); + } + + /* disp */ + if ( disp && *disp ) { + if ( (*disp)->namein ) { + MMG5_DEL_MEM(*mesh,(*disp)->namein); + } + + if ( (*disp)->nameout ) { + MMG5_DEL_MEM(*mesh,(*disp)->nameout); + } + } + + /* ls */ + if ( ls && *ls ) { + if ( (*ls)->namein ) { + MMG5_DEL_MEM(*mesh,(*ls)->namein); + } + + if ( (*ls)->nameout ) { + MMG5_DEL_MEM(*mesh,(*ls)->nameout); + } + } + + /* Fields */ + if ( sols ) { + for ( i=0; i<(*mesh)->nsols; ++i ) { + psl = (*sols) + i; + if ( psl->namein ) { + MMG5_DEL_MEM(*mesh,psl->namein); + } + if ( psl->nameout ) { + MMG5_DEL_MEM(*mesh,psl->nameout); + } + } + } + + return 1; +} From bffb562bd69cc196cd5f986fa7bcd698bae02d17 Mon Sep 17 00:00:00 2001 From: Fabien Salmon Date: Thu, 4 Dec 2025 16:53:51 +0100 Subject: [PATCH 2/2] Improvement for the swap operations in particular cases --- src/mmg2d/mmg2d1.c | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/mmg2d/mmg2d1.c b/src/mmg2d/mmg2d1.c index a7de6734b..d0e957ac7 100644 --- a/src/mmg2d/mmg2d1.c +++ b/src/mmg2d/mmg2d1.c @@ -653,6 +653,8 @@ MMG5_int MMG2D_swpmsh(MMG5_pMesh mesh,MMG5_pSol met,int typchk,double* velocity) maxit = 2; mesh->base++; + if (mesh->info.limit_angle <= 4*atan(1.)) mesh->info.fem = 0; + do { ns = 0; for (k=1; k<=mesh->nt; k++) { @@ -674,6 +676,9 @@ MMG5_int MMG2D_swpmsh(MMG5_pMesh mesh,MMG5_pSol met,int typchk,double* velocity) while ( ns > 0 && ++itinfo.imprim) > 5 || mesh->info.ddebug) && nns > 0 ) fprintf(stdout," %8" MMG5_PRId " edge swapped\n",nns); + + if (mesh->info.limit_angle <= 4*atan(1.)) mesh->info.fem = mesh->info.setfem; + return nns; }