From 73c879d0728a7922cd088a479da14a178587501d Mon Sep 17 00:00:00 2001 From: Fabien Salmon Date: Tue, 27 Jan 2026 15:44:37 +0100 Subject: [PATCH 1/4] Additional option to remesh only the elements of bad quality --- src/common/libmmgtypes.h | 4 + src/common/scalem.c | 2 +- src/mmg2d/API_functions_2d.c | 24 +++++- src/mmg2d/libmmg2d.h | 5 ++ src/mmg2d/libmmg2d_tools.c | 19 +++++ src/mmg2d/mmg2d1.c | 153 ++++++++++++++++++++++++++++++++++- 6 files changed, 202 insertions(+), 5 deletions(-) diff --git a/src/common/libmmgtypes.h b/src/common/libmmgtypes.h index fe75d6193..23340efeb 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 */ @@ -646,6 +649,7 @@ typedef struct { \f$k^th\f$ quadrilaterals are adjacent and share their edges \a j and \a l (resp.) */ int *ipar; /*!< Store indices of the local parameters */ + double *velocity; /*!< Velocity of the vertices when Lagrangian resolution */ MMG5_pPoint point; /*!< Pointer toward the \ref MMG5_Point structure */ MMG5_pxPoint xpoint; /*!< Pointer toward the \ref MMG5_xPoint structure */ MMG5_pTetra tetra; /*!< Pointer toward the \ref MMG5_Tetra structure */ 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..29d0e1735 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; + /* to adapat more thoroughly close to boundaries */ + 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/libmmg2d.h b/src/mmg2d/libmmg2d.h index f0097d566..c0881676d 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 (0mark=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/mmg2d1.c b/src/mmg2d/mmg2d1.c index 4b4f7704e..30c385674 100644 --- a/src/mmg2d/mmg2d1.c +++ b/src/mmg2d/mmg2d1.c @@ -50,7 +50,7 @@ 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; @@ -107,6 +107,137 @@ 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) { + + 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]; + double factor_min = mesh->info.min[2]; + + 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 (mesh->velocity && p1 <= mesh->mark && p2 <= mesh->mark && p3 <= mesh->mark && pp1->tmp == -1 && pp2->tmp == -1 && pp3->tmp == -1) + mean_velocity = ( sqrt(pow(mesh->velocity[p1-1],2) + pow(mesh->velocity[p1-1+mesh->mark],2)) + + sqrt(pow(mesh->velocity[p2-1],2) + pow(mesh->velocity[p2-1+mesh->mark],2)) + + sqrt(pow(mesh->velocity[p3-1],2) + pow(mesh->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] + mesh->velocity[p1-1]*dt; + x2 = pp2->c[0] + mesh->velocity[p2-1]*dt; + x3 = pp3->c[0] + mesh->velocity[p3-1]*dt; + y1 = pp1->c[1] + mesh->velocity[p1-1+mesh->mark]*dt; + y2 = pp2->c[1] + mesh->velocity[p2-1+mesh->mark]*dt; + y3 = pp3->c[1] + mesh->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_pTria pt; @@ -130,6 +261,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) > mesh->info.limit_angle) continue; + /* Check if pt should be cut */ pt->flag = 0; @@ -450,6 +583,8 @@ MMG5_int MMG2D_colelt(MMG5_pMesh mesh,MMG5_pSol met,int typchk, double lmax) { pt = &mesh->tria[k]; if ( !MG_EOK(pt) || pt->ref < 0 ) continue; + if (mesh->info.limit_angle <= 4*atan(1.) && MMG2D_minangle(mesh,k) > 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 @@ -523,12 +658,16 @@ MMG5_int MMG2D_swpmsh(MMG5_pMesh mesh,MMG5_pSol met,int typchk) { 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++) { pt = &mesh->tria[k]; if ( !MG_EOK(pt) || pt->ref < 0 ) continue; + if (mesh->info.limit_angle <= 4*atan(1.) && MMG2D_minangle(mesh,k) > 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) ) { @@ -542,6 +681,9 @@ MMG5_int MMG2D_swpmsh(MMG5_pMesh mesh,MMG5_pSol met,int typchk) { 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; } @@ -653,6 +795,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) > mesh->info.limit_angle) continue; + imax = -1; lmax = -1.0; for (i=0; i<3; i++) { @@ -715,6 +859,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) > 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; @@ -727,7 +873,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); @@ -760,6 +906,9 @@ MMG5_int MMG2D_movtri(MMG5_pMesh mesh,MMG5_pSol met,int maxit,int8_t improve) { **/ int MMG2D_mmg2d1n(MMG5_pMesh mesh,MMG5_pSol met) { + 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"); From 98b6d5e23d6d5fa04626ab1c2a3dcc6d017966ff Mon Sep 17 00:00:00 2001 From: Fabien Salmon Date: Thu, 29 Jan 2026 10:21:45 +0100 Subject: [PATCH 2/4] Cosmetic correction of the pull request #321 --- src/common/libmmgtypes.h | 4 ++-- src/mmg2d/API_functions_2d.c | 10 +++++----- src/mmg2d/libmmg2d.h | 2 +- src/mmg2d/libmmg2d_tools.c | 12 ++++++------ src/mmg2d/mmg2d1.c | 33 ++++++++++++++++++++++++++++----- 5 files changed, 42 insertions(+), 19 deletions(-) diff --git a/src/common/libmmgtypes.h b/src/common/libmmgtypes.h index 23340efeb..c39719e94 100644 --- a/src/common/libmmgtypes.h +++ b/src/common/libmmgtypes.h @@ -528,8 +528,8 @@ typedef struct { 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 isotropic_pt_relocation; /*!< force the use of the isotropic point relocation */ + int bdy_adaptation; /*!< extend the remeshing close to boundaries when limit_angle is activated */ 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/mmg2d/API_functions_2d.c b/src/mmg2d/API_functions_2d.c index 29d0e1735..7d3b6c103 100644 --- a/src/mmg2d/API_functions_2d.c +++ b/src/mmg2d/API_functions_2d.c @@ -103,13 +103,13 @@ void MMG2D_Init_parameters(MMG5_pMesh mesh) { mesh->info.ls = MMG5_LS; /* xreg relaxation parameter value */ mesh->info.lxreg = MMG5_XREG; - /* [0/1] ,avoid/enforce istropic remeshing even with anisotropic metric */ - mesh->info.isotropic = MMG5_OFF; + /* [0/1] ,avoid/enforce isotropic smoothing even with anisotropic metric */ + mesh->info.isotropic_pt_relocation = MMG5_OFF; /* limit angle to avoid remeshing some good triangles */ mesh->info.limit_angle = 5.*atan(1.); /* Ridge detection */ mesh->info.dhd = MMG5_ANGEDG; - /* to adapat more thoroughly close to boundaries */ + /* to adapt more thoroughly close to boundaries */ mesh->info.bdy_adaptation = MMG5_OFF; } @@ -170,8 +170,8 @@ 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; + case MMG2D_IPARAM_isotropic_smoothing : + mesh->info.isotropic_pt_relocation = val; break; case MMG2D_IPARAM_lag : #ifdef USE_ELAS diff --git a/src/mmg2d/libmmg2d.h b/src/mmg2d/libmmg2d.h index c0881676d..9ca331ade 100644 --- a/src/mmg2d/libmmg2d.h +++ b/src/mmg2d/libmmg2d.h @@ -145,7 +145,7 @@ extern "C" { 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_IPARAM_isotropic_smoothing, /*!< [0/1], Avoid/enforce isotropic smoothing 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) */ }; diff --git a/src/mmg2d/libmmg2d_tools.c b/src/mmg2d/libmmg2d_tools.c index c48b53af4..f69964962 100644 --- a/src/mmg2d/libmmg2d_tools.c +++ b/src/mmg2d/libmmg2d_tools.c @@ -149,7 +149,7 @@ int MMG2D_parsar(int argc,char *argv[],MMG5_pMesh mesh,MMG5_pSol met,MMG5_pSol s return 0; break; case 'b': - if ( !strcmp(argv[i],"-bdy-adaptation") ) { + if ( !strcmp(argv[i],"-bdyadaptation") ) { if ( !MMG2D_Set_iparameter(mesh,met,MMG2D_IPARAM_bdy_adaptation,1) ) return 0; } @@ -199,10 +199,10 @@ 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") ) { + else if ( !strcmp(argv[i],"-hminfactor") ) { param = MMG2D_DPARAM_hmin_factor; } - else if ( !strcmp(argv[i],"-hmax_factor") ) { + else if ( !strcmp(argv[i],"-hmaxfactor") ) { param = MMG2D_DPARAM_hmax_factor; } else { @@ -246,8 +246,8 @@ 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) ) + else if ( !strcmp(argv[i],"-isotropicsmoothing") ) { + if ( !MMG2D_Set_iparameter(mesh,met,MMG2D_IPARAM_isotropic_smoothing,1) ) return 0; } else { @@ -293,7 +293,7 @@ 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 ) { + else if ( !strcmp(argv[i],"-limitangle") && ++i < argc ) { if ( !MMG2D_Set_dparameter(mesh,met,MMG2D_DPARAM_limit_angle,atof(argv[i])) ) return 0; } diff --git a/src/mmg2d/mmg2d1.c b/src/mmg2d/mmg2d1.c index 30c385674..471411a5e 100644 --- a/src/mmg2d/mmg2d1.c +++ b/src/mmg2d/mmg2d1.c @@ -107,7 +107,18 @@ int MMG2D_anatri(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) { return 1; } -// Minimum angle function +/** + * \param mesh pointer to the mesh structure. + * \param k triangle index. + * \return the minimal angle of a triangle or 0 if the triangle size + * is not conform to the tolerance sizes. + * + * Compute the minimal angle of triangle k. If the velocity of each + * node is given, the extrapolated minimal angle of the triangle + * is also calculated. If the edge lengths are too small of too + * large, 0 is returned as we want this triangle to be remeshed. + * + */ double MMG2D_minangle(MMG5_pMesh mesh, int k) { MMG5_pTria pt = &mesh->tria[k]; @@ -129,11 +140,13 @@ double MMG2D_minangle(MMG5_pMesh mesh, int k) { 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 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 (mesh->velocity && p1 <= mesh->mark && p2 <= mesh->mark && p3 <= mesh->mark && pp1->tmp == -1 && pp2->tmp == -1 && pp3->tmp == -1) + if (mesh->velocity && p1 <= mesh->mark && p2 <= mesh->mark && p3 <= mesh->mark + && pp1->tmp == -1 && pp2->tmp == -1 && pp3->tmp == -1) mean_velocity = ( sqrt(pow(mesh->velocity[p1-1],2) + pow(mesh->velocity[p1-1+mesh->mark],2)) + sqrt(pow(mesh->velocity[p2-1],2) + pow(mesh->velocity[p2-1+mesh->mark],2)) + sqrt(pow(mesh->velocity[p3-1],2) + pow(mesh->velocity[p3-1+mesh->mark],2)) )/3.; @@ -149,6 +162,7 @@ double MMG2D_minangle(MMG5_pMesh mesh, int k) { y3 = pp3->c[1]; } else { + // Extrapolated triangle location based on a fictitious time step dt = 0.25*square_root_area / mean_velocity; x1 = pp1->c[0] + mesh->velocity[p1-1]*dt; x2 = pp2->c[0] + mesh->velocity[p2-1]*dt; @@ -158,14 +172,17 @@ double MMG2D_minangle(MMG5_pMesh mesh, int k) { y3 = pp3->c[1] + mesh->velocity[p3-1+mesh->mark]*dt; } + // Compute the length of the triangle edges 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.); + // Compute the length of the fictitious triangle edges 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.); + // Check if the edge lengths satisfy the extended tolerances double min = length[0]; double max = length[0]; int min_index = 0; @@ -198,10 +215,12 @@ double MMG2D_minangle(MMG5_pMesh mesh, int k) { 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 the triangle is on a boundary, the angle tolerance is divided by two + // to extend remeshing close to boundaries where the mesh can be often squeezed if (mesh->info.bdy_adaptation) { if ((pp1->tag & MG_BDY) || (pp2->tag & MG_BDY) || (pp3->tag & MG_BDY)) coef = c1; else { + // The neighbours of boundary triangles also have a reduced tolerance for(int i = 0; i < 3; i++ ) { int k1 = adja[i]/3; @@ -222,6 +241,7 @@ double MMG2D_minangle(MMG5_pMesh mesh, int k) { } } + // Minimal angle 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)); @@ -872,8 +892,11 @@ MMG5_int MMG2D_movtri(MMG5_pMesh mesh,MMG5_pSol met,int maxit,int8_t improve) { ier = MMG2D_movedgpt(mesh,met,ilist,list,improve); if ( ier ) ns++; } + else if (mesh->info.isotropic_pt_relocation) { + ier = MMG2D_movintpt(mesh,met,ilist,list,improve); + } else { - if ( met->size == 3 && met->m && !mesh->info.isotropic) + if ( met->size == 3 && met->m ) ier = MMG2D_movintpt_ani(mesh,met,ilist,list,improve); else ier = MMG2D_movintpt(mesh,met,ilist,list,improve); From 0d12d47ba78f3d1bfd63db55be349e89e87ab215 Mon Sep 17 00:00:00 2001 From: Fabien Salmon Date: Fri, 30 Jan 2026 10:52:39 +0100 Subject: [PATCH 3/4] Modification of variable names --- src/common/libmmgtypes.h | 3 ++- src/common/scalem.c | 2 +- src/mmg2d/API_functions_2d.c | 6 ++--- src/mmg2d/mmg2d1.c | 47 ++++++++++++++++++++---------------- 4 files changed, 32 insertions(+), 26 deletions(-) diff --git a/src/common/libmmgtypes.h b/src/common/libmmgtypes.h index c39719e94..5b0c76968 100644 --- a/src/common/libmmgtypes.h +++ b/src/common/libmmgtypes.h @@ -525,6 +525,7 @@ typedef struct { 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 */ + double relative_min_tolerance, relative_max_tolerance; /* < Tolerance on the minimal and maximal element sizes */ 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) */ @@ -649,7 +650,7 @@ typedef struct { \f$k^th\f$ quadrilaterals are adjacent and share their edges \a j and \a l (resp.) */ int *ipar; /*!< Store indices of the local parameters */ - double *velocity; /*!< Velocity of the vertices when Lagrangian resolution */ + double *lagrangian_velocity; /*!< Velocity of the vertices when Lagrangian resolution */ MMG5_pPoint point; /*!< Pointer toward the \ref MMG5_Point structure */ MMG5_pxPoint xpoint; /*!< Pointer toward the \ref MMG5_xPoint structure */ MMG5_pTetra tetra; /*!< Pointer toward the \ref MMG5_Tetra structure */ diff --git a/src/common/scalem.c b/src/common/scalem.c index 17da1b2e8..33155d91b 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.; - if (mesh->dim == 3) mesh->info.min[2]= 0.; + 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 7d3b6c103..2d018f400 100644 --- a/src/mmg2d/API_functions_2d.c +++ b/src/mmg2d/API_functions_2d.c @@ -106,7 +106,7 @@ void MMG2D_Init_parameters(MMG5_pMesh mesh) { /* [0/1] ,avoid/enforce isotropic smoothing even with anisotropic metric */ mesh->info.isotropic_pt_relocation = MMG5_OFF; /* limit angle to avoid remeshing some good triangles */ - mesh->info.limit_angle = 5.*atan(1.); + mesh->info.limit_angle = -1.; // Deactivated when negative or > PI/3 /* Ridge detection */ mesh->info.dhd = MMG5_ANGEDG; /* to adapt more thoroughly close to boundaries */ @@ -361,10 +361,10 @@ int MMG2D_Set_dparameter(MMG5_pMesh mesh, MMG5_pSol sol, int dparam, double val) mesh->info.hausd = val; break; case MMG2D_DPARAM_hmin_factor : - mesh->info.min[2] = val; + mesh->info.relative_min_tolerance = val; break; case MMG2D_DPARAM_hmax_factor : - mesh->info.max[2] = val; + mesh->info.relative_max_tolerance = val; break; case MMG2D_DPARAM_limit_angle : mesh->info.limit_angle = val; diff --git a/src/mmg2d/mmg2d1.c b/src/mmg2d/mmg2d1.c index 471411a5e..5cd84a810 100644 --- a/src/mmg2d/mmg2d1.c +++ b/src/mmg2d/mmg2d1.c @@ -50,7 +50,7 @@ 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 && mesh->info.limit_angle >= 4.*atan(1.)) { + if ( !mesh->info.noinsert && (mesh->info.limit_angle >= M_PI/3. || mesh->info.limit_angle <= 0.)) { /* Memory free */ MMG5_DEL_MEM(mesh,mesh->adja); mesh->adja = 0; @@ -129,8 +129,8 @@ double MMG2D_minangle(MMG5_pMesh mesh, int k) { double c1 = 1.35*mesh->info.limit_angle; if (c1 > 1.) c1 = 1.; double coef = 1.; - double factor_max = mesh->info.max[2]; - double factor_min = mesh->info.min[2]; + double factor_max = mesh->info.relative_max_tolerance; + double factor_min = mesh->info.relative_min_tolerance; p1 = pt->v[0]; p2 = pt->v[1]; @@ -143,13 +143,14 @@ double MMG2D_minangle(MMG5_pMesh mesh, int k) { 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.; + double* velocity = &mesh->lagrangian_velocity[0]; // Velocity is only considered for already existing points - if (mesh->velocity && p1 <= mesh->mark && p2 <= mesh->mark && p3 <= mesh->mark - && pp1->tmp == -1 && pp2->tmp == -1 && pp3->tmp == -1) - mean_velocity = ( sqrt(pow(mesh->velocity[p1-1],2) + pow(mesh->velocity[p1-1+mesh->mark],2)) + - sqrt(pow(mesh->velocity[p2-1],2) + pow(mesh->velocity[p2-1+mesh->mark],2)) + - sqrt(pow(mesh->velocity[p3-1],2) + pow(mesh->velocity[p3-1+mesh->mark],2)) )/3.; + if (velocity && 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) { @@ -164,12 +165,12 @@ double MMG2D_minangle(MMG5_pMesh mesh, int k) { else { // Extrapolated triangle location based on a fictitious time step dt = 0.25*square_root_area / mean_velocity; - x1 = pp1->c[0] + mesh->velocity[p1-1]*dt; - x2 = pp2->c[0] + mesh->velocity[p2-1]*dt; - x3 = pp3->c[0] + mesh->velocity[p3-1]*dt; - y1 = pp1->c[1] + mesh->velocity[p1-1+mesh->mark]*dt; - y2 = pp2->c[1] + mesh->velocity[p2-1+mesh->mark]*dt; - y3 = pp3->c[1] + mesh->velocity[p3-1+mesh->mark]*dt; + 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; } // Compute the length of the triangle edges @@ -281,7 +282,7 @@ 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) > mesh->info.limit_angle) continue; + if (mesh->info.limit_angle >= 0. && mesh->info.limit_angle <= M_PI/3. && MMG2D_minangle(mesh,k) > mesh->info.limit_angle) continue; /* Check if pt should be cut */ pt->flag = 0; @@ -603,7 +604,7 @@ MMG5_int MMG2D_colelt(MMG5_pMesh mesh,MMG5_pSol met,int typchk, double lmax) { pt = &mesh->tria[k]; if ( !MG_EOK(pt) || pt->ref < 0 ) continue; - if (mesh->info.limit_angle <= 4*atan(1.) && MMG2D_minangle(mesh,k) > mesh->info.limit_angle) continue; + if (mesh->info.limit_angle >= 0. && mesh->info.limit_angle <= M_PI/3. && MMG2D_minangle(mesh,k) > 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 @@ -678,7 +679,10 @@ MMG5_int MMG2D_swpmsh(MMG5_pMesh mesh,MMG5_pSol met,int typchk) { maxit = 2; mesh->base++; - if (mesh->info.limit_angle <= 4*atan(1.)) mesh->info.fem = 0; + // If some triangles are not adapted using the limiting option on the triangle minimal angle, + // using fem = 1 here can lead to very poor triangles that remain in the final mesh. + // So, in this case, we enforce the use of fem = 0. + if (mesh->info.limit_angle >= 0. && mesh->info.limit_angle <= M_PI/3.) mesh->info.fem = 0; do { ns = 0; @@ -686,7 +690,7 @@ 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) > mesh->info.limit_angle) continue; + if (mesh->info.limit_angle >= 0. && mesh->info.limit_angle <= M_PI/3. && MMG2D_minangle(mesh,k) > mesh->info.limit_angle) continue; for (i=0; i<3; i++) { if ( MG_SIN(pt->tag[i]) || MG_EDG(pt->tag[i]) ) continue; @@ -702,7 +706,8 @@ MMG5_int MMG2D_swpmsh(MMG5_pMesh mesh,MMG5_pSol met,int typchk) { if ( (abs(mesh->info.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; + // The fem value is restored after the swaps + if (mesh->info.limit_angle >= 0. && mesh->info.limit_angle <= M_PI/3.) mesh->info.fem = mesh->info.setfem; return nns; } @@ -815,7 +820,7 @@ 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) > mesh->info.limit_angle) continue; + if (mesh->info.limit_angle >= 0. && mesh->info.limit_angle <= M_PI/3. && MMG2D_minangle(mesh,k) > mesh->info.limit_angle) continue; imax = -1; lmax = -1.0; @@ -879,7 +884,7 @@ 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) > mesh->info.limit_angle) continue; + if (mesh->info.limit_angle >= 0. && mesh->info.limit_angle <= M_PI/3. && MMG2D_minangle(mesh,k) > mesh->info.limit_angle) continue; for (i=0; i<3; i++) { p0 = &mesh->point[pt->v[i]]; From 7787afcb0b507aafc87e544bb1df61d8ff13f28d Mon Sep 17 00:00:00 2001 From: Fabien Salmon Date: Fri, 30 Jan 2026 11:12:39 +0100 Subject: [PATCH 4/4] Initialize the relative min and max tolerances --- src/mmg2d/API_functions_2d.c | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/mmg2d/API_functions_2d.c b/src/mmg2d/API_functions_2d.c index 2d018f400..ac6ec8db0 100644 --- a/src/mmg2d/API_functions_2d.c +++ b/src/mmg2d/API_functions_2d.c @@ -107,6 +107,9 @@ void MMG2D_Init_parameters(MMG5_pMesh mesh) { mesh->info.isotropic_pt_relocation = MMG5_OFF; /* limit angle to avoid remeshing some good triangles */ mesh->info.limit_angle = -1.; // Deactivated when negative or > PI/3 + /* Tolerances on the element sizes when the limit_angle is activated */ + mesh->info.relative_min_tolerance = 1.e10; + mesh->info.relative_max_tolerance = 1.e10; /* Ridge detection */ mesh->info.dhd = MMG5_ANGEDG; /* to adapt more thoroughly close to boundaries */