From 1185009c79ad459ea6f9e39f523dccaa7bae5836 Mon Sep 17 00:00:00 2001 From: Algiane Froehly Date: Thu, 20 Oct 2022 16:53:01 +0200 Subject: [PATCH 1/5] Do not try anymore to insert nodes using dichotomy if surface doesn't come from level-set discretization. --- src/mmg3d/mmg3d1.c | 99 +++++++++++++++++++++++++++++++++++++++------- 1 file changed, 84 insertions(+), 15 deletions(-) diff --git a/src/mmg3d/mmg3d1.c b/src/mmg3d/mmg3d1.c index 96e543ace..cfafbbef2 100644 --- a/src/mmg3d/mmg3d1.c +++ b/src/mmg3d/mmg3d1.c @@ -89,35 +89,73 @@ void MMG5_tet2tri(MMG5_pMesh mesh,MMG5_int k,int8_t ie,MMG5_Tria *ptt) { * * Find acceptable position for splitting. * + * \remark Historically Mmg tried to insert points by dichotomy whatever the + * type of surface considered but inserting points elsewhere than on the ideal + * surface creates a bad approximation of the boundary so we don't want to do it + * on surfaces whose definition we are sure of (general case). The insertion by + * dichotomy can however make sense in the case of the discretization of an + * iso-value during a shape optimization process, where the definition of the + * surface is temporary and will continue to evolve during the iterations. For + * backward compatibility for these applications, we continue to test dichotomy + * insertion along iso surfaces in -ls mode. + * */ int MMG3D_dichoto(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int *vx) { MMG5_pTetra pt; + MMG5_pxTetra pxt; MMG5_pPoint pa,pb,ps; double o[6][3],p[6][3]; float to,tp,t; int ier,it,maxit; MMG5_int ia,ib; - int8_t i; + int8_t i,iso_flag; + + /** Check for insertion by dichotomy only along iso surfaces */ + if ( !mesh->info.iso ) { + return 0; + } ier = 1; pt = &mesh->tetra[k]; + if ( !pt->xt ) { + /** No way to know if edge is along iso surface if we don't work from a + * tetra with a boundary face */ + return 0; + } - /* get point on surface and along segment for edge split */ + pxt = &mesh->xtetra[pt->xt]; + iso_flag = 0; + for (i=0; i<4; i++) { + /* Flag edges belonging to iso surf */ + if ( pxt->ftag[i] & MG_ISO ) { + MG_SET(iso_flag,MMG5_iarf[i][0]); + MG_SET(iso_flag,MMG5_iarf[i][1]); + MG_SET(iso_flag,MMG5_iarf[i][2]); + } + } + + if ( !iso_flag ) { + /* None of the edges is along iso surface: we don't want to compute the + * dichotomy */ + return 0; + } + + /** Get point on surface and along segment for edge split */ for (i=0; i<6; i++) { memset(p[i],0,3*sizeof(double)); memset(o[i],0,3*sizeof(double)); - if ( vx[i] > 0 ) { + if ( vx[i] > 0 && MG_GET(iso_flag,i) ) { + ps = &mesh->point[vx[i]]; + p[i][0] = ps->c[0]; + p[i][1] = ps->c[1]; + p[i][2] = ps->c[2]; ia = pt->v[MMG5_iare[i][0]]; ib = pt->v[MMG5_iare[i][1]]; pa = &mesh->point[ia]; pb = &mesh->point[ib]; - ps = &mesh->point[vx[i]]; o[i][0] = 0.5 * (pa->c[0] + pb->c[0]); o[i][1] = 0.5 * (pa->c[1] + pb->c[1]); o[i][2] = 0.5 * (pa->c[2] + pb->c[2]); - p[i][0] = ps->c[0]; - p[i][1] = ps->c[1]; - p[i][2] = ps->c[2]; } } maxit = 4; @@ -125,10 +163,10 @@ int MMG3D_dichoto(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int *vx) { tp = 1.0; to = 0.0; do { - /* compute new position */ + /** Compute new position */ t = 0.5 * (tp + to); for (i=0; i<6; i++) { - if ( vx[i] > 0 ) { + if ( vx[i] > 0 && MG_GET(iso_flag,i) ) { ps = &mesh->point[vx[i]]; ps->c[0] = o[i][0] + t*(p[i][0] - o[i][0]); ps->c[1] = o[i][1] + t*(p[i][1] - o[i][1]); @@ -176,11 +214,11 @@ int MMG3D_dichoto(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int *vx) { tp = t; } while ( ++it < maxit ); - /* restore coords of last valid pos. */ + /** Restore coords of last valid pos. */ if ( !ier ) { t = to; for (i=0; i<6; i++) { - if ( vx[i] > 0 ) { + if ( vx[i] > 0 && MG_GET(iso_flag,i) ) { ps = &mesh->point[vx[i]]; ps->c[0] = o[i][0] + t*(p[i][0] - o[i][0]); ps->c[1] = o[i][1] + t*(p[i][1] - o[i][1]); @@ -189,7 +227,7 @@ int MMG3D_dichoto(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int *vx) { } } - /* For very ill-shaped elements, we can have no valid position */ + /** For very ill-shaped elements, we can have no valid position */ switch (pt->flag) { case 1: case 2: case 4: case 8: case 16: case 32: ier = MMG3D_split1_sim(mesh,met,k,vx); @@ -242,19 +280,50 @@ int MMG3D_dichoto(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int *vx) { * Find acceptable position for MMG5_split1b, passing the shell of * considered edge, starting from o point. * + * \remark Historically Mmg tried to insert points by dichotomy whatever the + * type of surface considered but inserting points elsewhere than on the ideal + * surface creates a bad approximation of the boundary so we don't want to do it + * on surfaces whose definition we are sure of (general case). The insertion by + * dichotomy can however make sense in the case of the discretization of an + * iso-value during a shape optimization process, where the definition of the + * surface is temporary and will continue to evolve during the iterations. For + * backward compatibility for these applications, we continue to test dichotomy + * insertion along iso surfaces in -ls mode. + * */ int MMG3D_dichoto1b(MMG5_pMesh mesh,MMG5_pSol met,int64_t *list,int ret,MMG5_int ip) { MMG5_pTetra pt; + MMG5_pxTetra pxt; MMG5_pPoint p0,p1,ppt; int it,maxit; MMG5_int iel,np,nq; double m[3],o[3],tp,to,t; - int8_t ia,ier; + int8_t ia,ier,iso_edge; + + if ( !mesh->info.iso ) { + /** Check for insertion by dichotomy only along iso surfaces (so we can + * leave in iso mode) */ + return 0; + } iel = list[0] / 6; ia = list[0] % 6; pt = &mesh->tetra[iel]; + if ( !pt->xt ) { + /** No way to know if edge is along iso surface if we don't work from a + * tetra with a boundary face */ + return 0; + } + + pxt = &mesh->xtetra[pt->xt]; + iso_edge = (pxt->ftag[MMG5_ifar[ia][0]] & MG_ISO) || (pxt->ftag[MMG5_ifar[ia][0]] & MG_ISO); + + if ( !iso_edge ) { + /* Edge is not along iso surface: we don't want to compute the dichotomy */ + return 0; + } + np = pt->v[MMG5_iare[ia][0]]; nq = pt->v[MMG5_iare[ia][1]]; p0 = &mesh->point[np]; @@ -1163,7 +1232,7 @@ MMG5_anatetv(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) { MMG5_pPar par; double ll,o[3],ux,uy,uz,hma2,mincal; int l,memlack,ier; - MMG5_int src,vx[6],ip,ip1,ip2,k,ne,ns,nap; + MMG5_int src,vx[6],ip,ip1,ip2,k,ne,ns,nap; int8_t i,j,ia; /** 1. analysis */ @@ -2046,7 +2115,7 @@ MMG3D_anatets_iso(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) { MMG5_INCREASE_MEM_MESSAGE(); MMG3D_delPatternPts(mesh,hash);return -1; ,o,MG_BDY,src); - // Now pb->p contain a wrong memory address. + // Now pb->p contains a wrong memory address. pb.p[0] = &mesh->point[ptt.v[0]]; pb.p[1] = &mesh->point[ptt.v[1]]; pb.p[2] = &mesh->point[ptt.v[2]]; From e20a19329408e335ebc270c4f8eecf1915bb298b Mon Sep 17 00:00:00 2001 From: Algiane Froehly Date: Fri, 4 Nov 2022 16:38:46 +0100 Subject: [PATCH 2/5] Clean split_3d.c file before modifs. --- src/mmg3d/split_3d.c | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/mmg3d/split_3d.c b/src/mmg3d/split_3d.c index 528e89018..6bad703a5 100644 --- a/src/mmg3d/split_3d.c +++ b/src/mmg3d/split_3d.c @@ -1199,7 +1199,7 @@ void MMG3D_update_qual(MMG5_pMesh mesh,MMG5_pSol met,const int ne, } else if ( (!met) || (!met->m) ) { /* in ls mode + -A option, orcal calls caltet_ani that fails */ - for (i=0; iqual=MMG5_caltet_iso(mesh,met,pt[i]); } } @@ -1839,8 +1839,6 @@ int MMG3D_split3cone_sim(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int vx[6] if ( pt->v[tau[3]] < pt->v[ib] ) { ib = tau[3]; } - else { - } } /* Check orientation of the 4 newly created tets */ From 9e40d3c27560fac10cecd4d614488335f9b74bcf Mon Sep 17 00:00:00 2001 From: Algiane Froehly Date: Fri, 4 Nov 2022 16:40:09 +0100 Subject: [PATCH 3/5] Fix wrong check in split5_sim. --- src/mmg3d/split_3d.c | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/mmg3d/split_3d.c b/src/mmg3d/split_3d.c index 6bad703a5..0a441ffaf 100644 --- a/src/mmg3d/split_3d.c +++ b/src/mmg3d/split_3d.c @@ -3988,6 +3988,8 @@ int MMG3D_split5_sim(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int vx[6]) { memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[0]] = vx[taued[2]]; pt0->v[tau[1]] = vx[taued[4]]; + pt0->v[tau[2]] = vx[taued[3]]; pt0->v[tau[3]] = vx[taued[5]]; + vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; From 71a370e76306a8fbd590a2aa1ae900a649e8db58 Mon Sep 17 00:00:00 2001 From: Algiane Froehly Date: Fri, 18 Nov 2022 09:30:25 +0100 Subject: [PATCH 4/5] MMG3D: Add checks of normal deviation regarding normals along MG_EDG points to non conformal splits to avoid projection issues (when ware are not anymore able to find a positive projection) in mmg3dBezierCP. It fixes issue #167. --- src/mmg3d/split_3d.c | 806 ++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 804 insertions(+), 2 deletions(-) diff --git a/src/mmg3d/split_3d.c b/src/mmg3d/split_3d.c index 0a441ffaf..0c9bda535 100644 --- a/src/mmg3d/split_3d.c +++ b/src/mmg3d/split_3d.c @@ -82,6 +82,156 @@ void MMG3D_split1_cfg(MMG5_int flag,uint8_t *tau,const uint8_t **taued) { } /** + * \param mesh pointer toward the mesh structure. + * \param pt pointer toward tetra used to simulate splitting (possible future tetra) + * \param k index of tetra used to simulate splitting. + * \param ifac index of boundary face for which we want to check normal projection. + * + * \return 1 if split can be performed (no ridge point or positive projection of + * at least one of the normals at each ridge point onto the normal at new + * triangle). 0 otherwise (for at least one ridge point, the 2 normals have + * negative projection onto the normal at new tria). + * + * Check the positivity of the projection of normals at possible ridge points + * onto the normal at new triangle along boundary face \a ifac. + * + */ +static inline +int MMG3D_check_ridpt_normal_proj_fac ( MMG5_pMesh mesh, MMG5_pTetra pt,MMG5_int k, int8_t ifac ) { + + /* Note that the tetra may have no xtetra if element is split by adjacency + * (even for more than one edge) */ + if ( !pt->xt ) { + return 1; + } + + /* Remark: xtetra is not updated regarding to the splitting that is simulated + * but it is not an issue as it is used only to check the face tag */ + MMG5_pxTetra pxt = &mesh->xtetra[pt->xt]; + + if ( !(pxt->ftag[ifac] & MG_BDY) ) { + /* nothing to check */ + return 1; + } + + MMG5_Tria ptt; + MMG5_tet2tri(mesh,k,ifac,&ptt); + double nt[3]; + MMG5_nortri(mesh,&ptt,nt); + + int8_t i; + for ( i=0; i<3; ++i ) { + MMG5_int ip = pt->v[MMG5_idir[ifac][i]]; + MMG5_pPoint ppt = &mesh->point[ip]; + if ( MG_SIN(ppt->tag) || !MG_EDG( ppt->tag ) ) { + continue; + } + + /* Non singular ridge point: check for normal projection */ + assert ( ppt->xp ); + MMG5_pxPoint pxp = &mesh->xpoint[ppt->xp]; + double ps = pxp->n1[0]*nt[0] + pxp->n1[1]*nt[1] + pxp->n1[2]*nt[2]; + if ( ps > 0. ) { + /* Normal projection is ok for this point */ + continue; + } + ps = pxp->n2[0]*nt[0] + pxp->n2[1]*nt[1] + pxp->n2[2]*nt[2]; + if ( ps <= 0. ) { + /* We are here only if ps is <= 0 too so it means that we don't want + * to insert the point */ + return 0; + } + } + return 1; +} + +/** + * \param mesh pointer toward the mesh structure. + * \param pt pointer toward tetra used to simulate splitting (one of the tetra + * that will be created by the splitting). + * \param k index of tetra used to simulate splitting. + * \param ia index of edge shared by the 2 faces to test. + * + * \return 1 if split can be performed (no ridge point or positive projection of + * at least one of the normals at each ridge point onto the normal at new + * triangle). 0 otherwise (for at least one ridge point, the 2 normals have + * negative projection onto the normal at new tria). + * + * Check the positivity of the projection of normals at possible ridge points + * onto the normales at the 2 new triangles belonging to new tetra \a pt and + * created on both sides of edge \a ia. + * + */ +static inline +int MMG3D_check_ridpt_normal_proj_edg ( MMG5_pMesh mesh, MMG5_pTetra pt,MMG5_int k, int8_t ia ) { + int8_t ifac; + + if ( !pt->xt ) { + return 1; + } + + /* If tetra has a boundary face, check the projection of the normals at + * ridge points onto new triangle normal (otherwise we will fail later in + * mmg3dBezierCP) */ + ifac = MMG5_ifar[ia][0]; + if ( !MMG3D_check_ridpt_normal_proj_fac(mesh,pt,k,ifac) ) { + return 0; + } + + ifac = MMG5_ifar[ia][1]; + if ( !MMG3D_check_ridpt_normal_proj_fac(mesh,pt,k,ifac) ) { + return 0; + } + return 1; +} + +/** + * \param mesh pointer toward the mesh structure. + * \param pt pointer toward tetra used to simulate splitting (one of the tetra + * that will be created by the splitting). + * \param k index of tetra used to simulate splitting. + * \param ifac1 index of face to test. + * \param ifac2 index of the second face to test. + * \param ifac3 index of the third face to test. + * + * \return 1 if split can be performed (no ridge point or positive projection of + * at least one of the normals at each ridge point onto the normal at new + * triangle). 0 otherwise (for at least one ridge point, the 2 normals have + * negative projection onto the normal at new tria). + * + * Check the positivity of the projection of normals at possible ridge points + * onto the normald at the 3 new triangles belonging to new tetra \a pt and + * faces with local index \a ifac1, \a ifac2 \a ifac 3. + * + */ +static inline +int MMG3D_check_ridpt_normal_proj_3fac ( MMG5_pMesh mesh,MMG5_pTetra pt,MMG5_int k, + int8_t ifac1,int8_t ifac2,int8_t ifac3 ) { + + + if ( !pt->xt ) { + return 1; + } + + /* If tetra has a boundary face, check the projection of the normals at + * ridge points onto new triangle normal (otherwise we will fail later in + * mmg3dBezierCP) */ + if ( !MMG3D_check_ridpt_normal_proj_fac(mesh,pt,k,ifac1) ) { + return 0; + } + + if ( !MMG3D_check_ridpt_normal_proj_fac(mesh,pt,k,ifac2) ) { + return 0; + } + + if ( !MMG3D_check_ridpt_normal_proj_fac(mesh,pt,k,ifac3) ) { + return 0; + } + + return 1; +} + +/** mmg3d_debug dense.mesh -sol dense.sol -ls 0 -hausd 0.5 -hmin 1 -hmax 5 -ar 80 -nofem * \param mesh pointer toward the mesh structure. * \param met pointer toward the metric structure. * \param k index of element to split. @@ -107,17 +257,30 @@ int MMG3D_split1_sim(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int vx[6]) { MMG3D_split1_cfg(pt->flag,tau,&taued); - /* Test volume of the two created tets */ + /* Test volume and normal deviation with regard to ridge points for the two + * created tets */ memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[1]] = vx[taued[0]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at the + * new triangle (otherwise we will fail later in mmg3dBezierCP) */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[0]) ) { + return 0; + } + memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[0]] = vx[taued[0]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at the + * new triangle (otherwise we will fail later in mmg3dBezierCP) */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[0]) ) { + return 0; + } + return 1; } @@ -1093,6 +1256,10 @@ int MMG3D_split2sf_sim(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int vx[6]){ pt0->v[tau[2]] = vx[taued[5]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + if ( !MMG3D_check_ridpt_normal_proj_3fac(mesh,pt0,0,tau[0],tau[1],tau[2]) ) { + return 0; + } + if ( imin == tau[1] ) { memcpy(pt0,pt,sizeof(MMG5_Tetra)); @@ -1101,22 +1268,53 @@ int MMG3D_split2sf_sim(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int vx[6]){ vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0]-tau[1]-tau[3] and + * tau[0]-tau[1]-tau[5], that are, the two faces on both sides of edge + * taued[0] */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[0]) ) { + return 0; + } + memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[3]] = vx[taued[5]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + + /* Check the projection of the normals at ridge points onto the normal at the + * new triangle (otherwise we will fail later in mmg3dBezierCP) */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[5]) ) { + return 0; + } } else { memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[3]] = vx[taued[4]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0]-tau[1]-tau[3] and + * tau[1]-tau[3]-tau[2], that are, the two faces on both sides of edge + * taued[4] */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[4]) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[1]] = vx[taued[4]]; pt0->v[tau[3]] = vx[taued[5]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0]-tau[1]-tau[3] and + * tau[1]-tau[3]-tau[2], that are, the two faces on both sides of edge + * taued[4] */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[5]) ) { + return 0; + } } return 1; } @@ -1396,21 +1594,53 @@ int MMG3D_split2_sim(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int vx[6]){ pt0->v[tau[1]] = vx[taued[0]]; pt0->v[tau[2]] = vx[taued[5]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0]-tau[1]-tau[3] and + * tau[0]-tau[3]-tau[2], that are, the two faces on both sides of edge + * taued[2] */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[2]) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[1]] = vx[taued[0]]; pt0->v[tau[3]] = vx[taued[5]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0]-tau[1]-tau[2] and + * tau[0]-tau[3]-tau[2], that are, the two faces on both sides of edge + * taued[1] */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[1]) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[0]] = vx[taued[0]]; pt0->v[tau[2]] = vx[taued[5]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0]-tau[1]-tau[3] and + * tau[1]-tau[3]-tau[2], that are, the two faces on both sides of edge + * taued[4] */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[2]) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[0]] = vx[taued[0]]; pt0->v[tau[3]] = vx[taued[5]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0]-tau[1]-tau[2] and + * tau[1]-tau[3]-tau[2], that are, the two faces on both sides of edge + * taued[4] */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[3]) ) { + return 0; + } return 1; } @@ -1604,18 +1834,33 @@ int MMG3D_split3_sim(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int vx[6]) { pt0->v[tau[2]] = vx[taued[1]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangles (otherwise we will fail later in mmg3dBezierCP) */ + if ( !MMG3D_check_ridpt_normal_proj_3fac(mesh,pt0,0,tau[2],tau[1],tau[3]) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[0]] = vx[taued[0]]; pt0->v[tau[2]] = vx[taued[3]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangles (otherwise we will fail later in mmg3dBezierCP) */ + if ( !MMG3D_check_ridpt_normal_proj_3fac(mesh,pt0,0,tau[2],tau[0],tau[3]) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[0]] = vx[taued[1]]; pt0->v[tau[1]] = vx[taued[3]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangles (otherwise we will fail later in mmg3dBezierCP) */ + if ( !MMG3D_check_ridpt_normal_proj_3fac(mesh,pt0,0,tau[3],tau[0],tau[1]) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[0]] = vx[taued[0]]; @@ -1623,6 +1868,11 @@ int MMG3D_split3_sim(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int vx[6]) { pt0->v[tau[2]] = vx[taued[1]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (useful only if some of the new points are ridges) */ + if ( !MMG3D_check_ridpt_normal_proj_fac(mesh,pt0,0,tau[3]) ) { + return 0; + } return 1; } @@ -1848,6 +2098,11 @@ int MMG3D_split3cone_sim(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int vx[6] pt0->v[tau[3]] = vx[taued[2]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangles (otherwise we will fail later in mmg3dBezierCP) */ + if ( !MMG3D_check_ridpt_normal_proj_3fac(mesh,pt0,0,tau[2],tau[1],tau[3]) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); if ( ia == tau[3] ) { @@ -1856,6 +2111,14 @@ int MMG3D_split3cone_sim(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int vx[6] pt0->v[tau[2]] = vx[taued[1]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0]-tau[2]-tau[3] and + * tau[0]-tau[1]-tau[3], that are, the two faces on both sides of edge + * taued[2] */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[2]) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); if ( ib == tau[1] ) { @@ -1863,11 +2126,27 @@ int MMG3D_split3cone_sim(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int vx[6] pt0->v[tau[2]] = vx[taued[1]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0]-tau[1]-tau[2] and + * tau[0]-tau[1]-tau[3], that are, the two faces on both sides of edge + * taued[0] */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[0]) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[0]] = vx[taued[1]] ; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0]-tau[2]-tau[3] and + * tau[0]-tau[1]-tau[2], that are, the two faces on both sides of edge + * taued[1] */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[1]) ) { + return 0; + } } else { assert(ib == tau[2]); @@ -1875,11 +2154,27 @@ int MMG3D_split3cone_sim(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int vx[6] pt0->v[tau[1]] = vx[taued[0]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0]-tau[2]-tau[3] and + * tau[0]-tau[1]-tau[2], that are, the two faces on both sides of edge + * taued[1] */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[1]) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[0]] = vx[taued[0]] ; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0]-tau[1]-tau[3] and + * tau[0]-tau[1]-tau[2], that are, the two faces on both sides of edge + * taued[0] */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[0]) ) { + return 0; + } } } else if (ia == tau[2] ) { @@ -1889,6 +2184,15 @@ int MMG3D_split3cone_sim(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int vx[6] vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0]-tau[3]-tau[2] and + * tau[0]-tau[1]-tau[2], that are, the two faces on both sides of edge + * taued[1] */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[1]) ) { + return 0; + } + memcpy(pt0,pt,sizeof(MMG5_Tetra)); if ( ib == tau[3] ) { pt0->v[tau[0]] = vx[taued[2]]; @@ -1896,10 +2200,28 @@ int MMG3D_split3cone_sim(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int vx[6] vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0]-tau[3]-tau[2] and + * tau[0]-tau[1]-tau[2], that are, the two faces on both sides of edge + * taued[1] */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[1]) ) { + return 0; + } + memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[0]] = vx[taued[0]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0]-tau[1]-tau[3] and + * tau[0]-tau[1]-tau[2], that are, the two faces on both sides of edge + * taued[0] */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[0]) ) { + return 0; + } } else { assert(ib == tau[1]); @@ -1908,10 +2230,28 @@ int MMG3D_split3cone_sim(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int vx[6] vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0]-tau[1]-tau[3] and + * tau[0]-tau[1]-tau[2], that are, the two faces on both sides of edge + * taued[0] */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[0]) ) { + return 0; + } + memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[0]] = vx[taued[2]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0]-tau[1]-tau[3] and + * tau[0]-tau[3]-tau[2], that are, the two faces on both sides of edge + * taued[2] */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[2]) ) { + return 0; + } } } else { @@ -1923,17 +2263,42 @@ int MMG3D_split3cone_sim(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int vx[6] vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0]-tau[1]-tau[3] and + * tau[0]-tau[1]-tau[2], that are, the two faces on both sides of edge + * taued[1] */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[0]) ) { + return 0; + } + memcpy(pt0,pt,sizeof(MMG5_Tetra)); if ( ib == tau[2] ) { pt0->v[tau[0]] = vx[taued[1]]; pt0->v[tau[3]] = vx[taued[2]] ; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0]-tau[1]-tau[2] and + * tau[0]-tau[3]-tau[2], that are, the two faces on both sides of edge + * taued[1] */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[1]) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[0]] = vx[taued[2]] ; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0]-tau[1]-tau[3] and + * tau[0]-tau[3]-tau[2], that are, the two faces on both sides of edge + * taued[2] */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[2]) ) { + return 0; + } } else { @@ -1944,10 +2309,28 @@ int MMG3D_split3cone_sim(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int vx[6] vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0]-tau[1]-tau[3] and + * tau[0]-tau[3]-tau[2], that are, the two faces on both sides of edge + * taued[2] */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[2]) ) { + return 0; + } + memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[0]] = vx[taued[1]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0]-tau[1]-tau[2] and + * tau[0]-tau[3]-tau[2], that are, the two faces on both sides of edge + * taued[1] */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[1]) ) { + return 0; + } } } @@ -2459,77 +2842,191 @@ int MMG3D_split3op_sim(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int vx[6]) pt0->v[ip0] = vx[ie1] ; pt0->v[ip1] = vx[ie0] ; pt0->v[ip3] = vx[ie5] ; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles ip0-ip1-ip2 and + * ip0-ip3-ip2, that are, the two faces on both sides of edge + * ie1 */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,ie1) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[ip0] = vx[ie0] ; pt0->v[ip3] = vx[ie5] ; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles ip0-ip1-ip2 and + * ip1-ip3-ip2, that are, the two faces on both sides of edge + * ie3 */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,ie3) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[ip0] = vx[ie0] ; pt0->v[ip2] = vx[ie5] ; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles ip3-ip1-ip2 and + * ip1-ip3-ip0, that are, the two faces on both sides of edge + * ie4 */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,ie4) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[ip1] = vx[ie0] ; pt0->v[ip2] = vx[ie1] ; pt0->v[ip3] = vx[ie5] ; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles ip0-ip1-ip2 and + * ip0-ip3-ip2, that are, the two faces on both sides of edge + * ie1 */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,ie1) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[ip1] = vx[ie0] ; pt0->v[ip2] = vx[ie5]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles ip0-ip3-ip2 and + * ip0-ip3-ip1, that are, the two faces on both sides of edge + * ie2 */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,ie2) ) { + return 0; + } } else if ( (imin12 == ip1) && (imin03 == ip0) ) { pt0->v[ip0] = vx[ie1] ; pt0->v[ip3] = vx[ie5] ; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles ip0, ip1 and ip2 */ + if ( !MMG3D_check_ridpt_normal_proj_3fac(mesh,pt0,0,ip0,ip1,ip2) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[ip0] = vx[ie0] ; pt0->v[ip2] = vx[ie1] ; pt0->v[ip3] = vx[ie5]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangle ip3 */ + if ( !MMG3D_check_ridpt_normal_proj_fac(mesh,pt0,0,ip3) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[ip0] = vx[ie0] ; pt0->v[ip2] = vx[ie5] ; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles ip1 ip2 ip3 and + * ip0 ip1 ip3, that are, the two faces on both sides of edge + * ie4 */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,ie4) ) { + return 0; + } + memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[ip1] = vx[ie0] ; pt0->v[ip2] = vx[ie5]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles ip0 ip2 ip3 and + * ip0 ip1 ip3, that are, the two faces on both sides of edge + * ie2 */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,ie2) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[ip1] = vx[ie0] ; pt0->v[ip2] = vx[ie1]; pt0->v[ip3] = vx[ie5]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles ip0 ip2 ip3 and + * ip0 ip1 ip2, that are, the two faces on both sides of edge + * ie1 */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,ie1) ) { + return 0; + } } else if ( (imin12 == ip2) && (imin03 == ip3) ) { pt0->v[ip1] = vx[ie0] ; pt0->v[ip2] = vx[ie1] ; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles ip1, ip2 and ip3 */ + if ( !MMG3D_check_ridpt_normal_proj_3fac(mesh,pt0,0,ip3,ip1,ip2) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[ip0] = vx[ie1] ; pt0->v[ip1] = vx[ie0] ; pt0->v[ip2] = vx[ie5]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangle ip1 */ + if ( !MMG3D_check_ridpt_normal_proj_fac(mesh,pt0,0,ip1) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[ip0] = vx[ie0] ; pt0->v[ip2] = vx[ie5] ; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles ip1 ip2 ip3 and + * ip0 ip1 ip3, that are, the two faces on both sides of edge + * ie4 */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,ie4) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[ip0] = vx[ie1] ; pt0->v[ip1] = vx[ie0]; pt0->v[ip3] = vx[ie5]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles ip0 ip2 ip3 and + * ip0 ip1 ip2, that are, the two faces on both sides of edge + * ie1 */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,ie1) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[ip0] = vx[ie0] ; pt0->v[ip3] = vx[ie5]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles ip0-ip1-ip2 and + * ip1-ip3-ip2, that are, the two faces on both sides of edge + * ie3 */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,ie3) ) { + return 0; + } + } else { assert((imin12 == ip1) && (imin03 == ip3)) ; @@ -2537,21 +3034,50 @@ int MMG3D_split3op_sim(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int vx[6]) pt0->v[ip1] = vx[ie0] ; pt0->v[ip2] = vx[ie1] ; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles ip1, ip2 and ip3 */ + if ( !MMG3D_check_ridpt_normal_proj_3fac(mesh,pt0,0,ip3,ip1,ip2) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[ip0] = vx[ie1] ; pt0->v[ip3] = vx[ie5] ; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles ip0, ip1 and ip3 */ + if ( !MMG3D_check_ridpt_normal_proj_3fac(mesh,pt0,0,ip0,ip1,ip3) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[ip0] = vx[ie0] ; pt0->v[ip2] = vx[ie1] ; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles ip0-ip1-ip2 and + * ip0-ip1-ip2, that are, the two faces on both sides of edge + * ie0 */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,ie0) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[ip0] = vx[ie1] ; pt0->v[ip2] = vx[ie5] ; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles ip0-ip2-ip3 and + * ip1-ip3-ip2, that are, the two faces on both sides of edge + * ie4 */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,ie4) ) { + return 0; + } + } return 1; @@ -3269,6 +3795,12 @@ int MMG3D_split4sf_sim(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int vx[6]) pt0->v[tau[3]] = vx[taued[2]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[1], tau[2] and tau[3] */ + if ( !MMG3D_check_ridpt_normal_proj_3fac(mesh,pt0,0,tau[1],tau[2],tau[3]) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[0]] = vx[taued[2]]; @@ -3277,6 +3809,12 @@ int MMG3D_split4sf_sim(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int vx[6]) pt0->v[tau[3]] = vx[taued[4]] ; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normal at triangle tau[0]-tau[1]-tau[3] */ + if ( !MMG3D_check_ridpt_normal_proj_fac(mesh,pt0,0,tau[2]) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); if ( imin12 == tau[1] ) { @@ -3285,12 +3823,29 @@ int MMG3D_split4sf_sim(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int vx[6]) pt0->v[tau[3]] = vx[taued[4]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0]-tau[1]-tau[3] and + * tau[0]-tau[1]-tau[2], that are, the two faces on both sides of edge + * taued[0] */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[0]) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[0]] = vx[taued[1]]; pt0->v[tau[3]] = vx[taued[4]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[1]-tau[2]-tau[3] and + * tau[0]-tau[1]-tau[2], that are, the two faces on both sides of edge + * taued[3] */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[3]) ) { + return 0; + } + } else { pt0->v[tau[0]] = vx[taued[1]]; @@ -3298,12 +3853,25 @@ int MMG3D_split4sf_sim(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int vx[6]) pt0->v[tau[3]] = vx[taued[4]] ; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normal at triangle tau[0]-tau[1]-tau[2] */ + if ( !MMG3D_check_ridpt_normal_proj_fac(mesh,pt0,0,tau[3]) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[0]] = vx[taued[0]]; pt0->v[tau[3]] = vx[taued[4]] ; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0],tau[2] and + * tau[3] */ + if ( !MMG3D_check_ridpt_normal_proj_3fac(mesh,pt0,0,tau[0],tau[2],tau[3]) ) { + return 0; + } } memcpy(pt0,pt,sizeof(MMG5_Tetra)); @@ -3313,12 +3881,26 @@ int MMG3D_split4sf_sim(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int vx[6]) pt0->v[tau[3]] = vx[taued[2]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normal at triangle tau[0]-tau[2]-tau[3] */ + if ( !MMG3D_check_ridpt_normal_proj_fac(mesh,pt0,0,tau[1]) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[0]] = vx[taued[2]]; pt0->v[tau[1]] = vx[taued[4]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0],tau[1] and + * tau[2] */ + if ( !MMG3D_check_ridpt_normal_proj_3fac(mesh,pt0,0,tau[0],tau[1],tau[2]) ) { + return 0; + } + } else { pt0->v[tau[0]] = vx[taued[2]]; @@ -3326,12 +3908,29 @@ int MMG3D_split4sf_sim(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int vx[6]) pt0->v[tau[2]] = vx[taued[1]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0]-tau[2]-tau[3] and + * tau[0]-tau[1]-tau[3], that are, the two faces on both sides of edge + * taued[2] */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[2]) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[0]] = vx[taued[1]]; pt0->v[tau[1]] = vx[taued[4]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0]-tau[2]-tau[3] and + * tau[1]-tau[2]-tau[3], that are, the two faces on both sides of edge + * taued[5] */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[5]) ) { + return 0; + } + } return 1; @@ -3589,35 +4188,77 @@ int MMG3D_split4op_sim(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int vx[6]) pt0->v[tau[2]] = vx[taued[3]]; pt0->v[tau[3]] = vx[taued[4]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0],tau[2] and + * tau[3] */ + if ( !MMG3D_check_ridpt_normal_proj_3fac(mesh,pt0,0,tau[0],tau[2],tau[3]) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[1]] = vx[taued[4]]; pt0->v[tau[2]] = vx[taued[3]]; pt0->v[tau[3]] = vx[taued[2]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normal at triangle tau[0]-tau[1]-tau[3] */ + if ( !MMG3D_check_ridpt_normal_proj_fac(mesh,pt0,0,tau[2]) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[1]] = vx[taued[3]]; pt0->v[tau[2]] = vx[taued[1]]; pt0->v[tau[3]] = vx[taued[2]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0]-tau[1]-tau[2] and + * tau[0]-tau[2]-tau[3], that are, the two faces on both sides of edge + * taued[1] */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[1]) ) { + return 0; + } } else { pt0->v[tau[2]] = vx[taued[1]]; pt0->v[tau[3]] = vx[taued[2]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[01,tau[2] and + * tau[3] */ + if ( !MMG3D_check_ridpt_normal_proj_3fac(mesh,pt0,0,tau[1],tau[2],tau[3]) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[0]] = vx[taued[1]]; pt0->v[tau[2]] = vx[taued[3]]; pt0->v[tau[3]] = vx[taued[2]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normal at triangle tau[0]-tau[1]-tau[2] */ + if ( !MMG3D_check_ridpt_normal_proj_fac(mesh,pt0,0,tau[3]) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[0]] = vx[taued[2]]; pt0->v[tau[2]] = vx[taued[3]]; pt0->v[tau[3]] = vx[taued[4]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0]-tau[1]-tau[3] and + * tau[1]-tau[2]-tau[3], that are, the two faces on both sides of edge + * taued[4] */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[4]) ) { + return 0; + } } memcpy(pt0,pt,sizeof(MMG5_Tetra)); @@ -3625,35 +4266,77 @@ int MMG3D_split4op_sim(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int vx[6]) pt0->v[tau[0]] = vx[taued[2]]; pt0->v[tau[1]] = vx[taued[4]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0],tau[1] and + * tau[2] */ + if ( !MMG3D_check_ridpt_normal_proj_3fac(mesh,pt0,0,tau[0],tau[1],tau[2]) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[0]] = vx[taued[2]]; pt0->v[tau[1]] = vx[taued[3]]; pt0->v[tau[3]] = vx[taued[4]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normal at triangle tau[1]-tau[2]-tau[3] */ + if ( !MMG3D_check_ridpt_normal_proj_fac(mesh,pt0,0,tau[0]) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[0]] = vx[taued[1]]; pt0->v[tau[1]] = vx[taued[3]]; pt0->v[tau[3]] = vx[taued[2]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0]-tau[1]-tau[2] and + * tau[0]-tau[2]-tau[3], that are, the two faces on both sides of edge + * taued[1] */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[1]) ) { + return 0; + } } else { pt0->v[tau[0]] = vx[taued[1]]; pt0->v[tau[1]] = vx[taued[3]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0],tau[1] and + * tau[3] */ + if ( !MMG3D_check_ridpt_normal_proj_3fac(mesh,pt0,0,tau[0],tau[1],tau[3]) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[0]] = vx[taued[2]]; pt0->v[tau[1]] = vx[taued[3]]; pt0->v[tau[2]] = vx[taued[1]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normal at triangle tau[0]-tau[2]-tau[3] */ + if ( !MMG3D_check_ridpt_normal_proj_fac(mesh,pt0,0,tau[1]) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[0]] = vx[taued[2]]; pt0->v[tau[1]] = vx[taued[4]]; pt0->v[tau[2]] = vx[taued[3]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0]-tau[1]-tau[3] and + * tau[1]-tau[2]-tau[3], that are, the two faces on both sides of edge + * taued[4] */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[4]) ) { + return 0; + } } return 1; @@ -3979,60 +4662,127 @@ int MMG3D_split5_sim(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int vx[6]) { pt0->v[tau[2]] = vx[taued[5]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0],tau[1] and + * tau[2] */ + if ( !MMG3D_check_ridpt_normal_proj_3fac(mesh,pt0,0,tau[0],tau[1],tau[2]) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[0]] = vx[taued[1]]; pt0->v[tau[1]] = vx[taued[3]]; pt0->v[tau[3]] = vx[taued[5]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0],tau[1] and + * tau[3] */ + if ( !MMG3D_check_ridpt_normal_proj_3fac(mesh,pt0,0,tau[0],tau[1],tau[3]) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[0]] = vx[taued[2]]; pt0->v[tau[1]] = vx[taued[4]]; pt0->v[tau[2]] = vx[taued[3]]; pt0->v[tau[3]] = vx[taued[5]]; - vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normal at triangle tau[1]-tau[2]-tau[3] */ + if ( !MMG3D_check_ridpt_normal_proj_fac(mesh,pt0,0,tau[0]) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[0]] = vx[taued[2]]; pt0->v[tau[1]] = vx[taued[3]]; pt0->v[tau[2]] = vx[taued[1]]; pt0->v[tau[3]] = vx[taued[5]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normal at triangle tau[0]-tau[2]-tau[3] */ + if ( !MMG3D_check_ridpt_normal_proj_fac(mesh,pt0,0,tau[1]) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); if ( imin == tau[0] ) { pt0->v[tau[2]] = vx[taued[3]]; pt0->v[tau[3]] = vx[taued[4]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0],tau[2] and + * tau[3] */ + if ( !MMG3D_check_ridpt_normal_proj_3fac(mesh,pt0,0,tau[0],tau[2],tau[3]) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[1]] = vx[taued[4]]; pt0->v[tau[2]] = vx[taued[3]]; pt0->v[tau[3]] = vx[taued[2]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normal at triangle tau[0]-tau[1]-tau[3] */ + if ( !MMG3D_check_ridpt_normal_proj_fac(mesh,pt0,0,tau[2]) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[1]] = vx[taued[3]]; pt0->v[tau[2]] = vx[taued[1]]; pt0->v[tau[3]] = vx[taued[2]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0]-tau[2]-tau[3] and + * tau[0]-tau[1]-tau[2], that are, the two faces on both sides of edge + * taued[1] */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[1]) ) { + return 0; + } } else { pt0->v[tau[2]] = vx[taued[1]]; pt0->v[tau[3]] = vx[taued[2]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[1],tau[2] and + * tau[3] */ + if ( !MMG3D_check_ridpt_normal_proj_3fac(mesh,pt0,0,tau[1],tau[2],tau[3]) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[0]] = vx[taued[2]]; pt0->v[tau[2]] = vx[taued[3]]; pt0->v[tau[3]] = vx[taued[4]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles tau[0]-tau[1]-tau[3] and + * tau[1]-tau[2]-tau[3], that are, the two faces on both sides of edge + * taued[4] */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[4]) ) { + return 0; + } memcpy(pt0,pt,sizeof(MMG5_Tetra)); pt0->v[tau[0]] = vx[taued[1]]; pt0->v[tau[2]] = vx[taued[3]]; pt0->v[tau[3]] = vx[taued[2]]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normal at triangle tau[0]-tau[1]-tau[2] */ + if ( !MMG3D_check_ridpt_normal_proj_fac(mesh,pt0,0,tau[3]) ) { + return 0; + } } return 1; @@ -4269,45 +5019,97 @@ int MMG3D_split6_sim(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int vx[6]) { pt0->v[1] = vx[0]; pt0->v[2] = vx[1]; pt0->v[3] = vx[2]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles vx[1],vx[2] and + * vx[3] */ + if ( !MMG3D_check_ridpt_normal_proj_3fac(mesh,pt0,0,vx[1],vx[2],vx[3]) ) { + return 0; + } /* Modify second tetra */ pt0->v[0] = vx[0]; pt0->v[2] = vx[3]; pt0->v[3] = vx[4]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles vx[0],vx[2] and + * vx[3] */ + if ( !MMG3D_check_ridpt_normal_proj_3fac(mesh,pt0,0,vx[0],vx[2],vx[3]) ) { + return 0; + } /* Modify 3rd tetra */ pt0->v[0] = vx[1]; pt0->v[1] = vx[3]; pt0->v[3] = vx[5]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles vx[0],vx[1] and + * vx[3] */ + if ( !MMG3D_check_ridpt_normal_proj_3fac(mesh,pt0,0,vx[0],vx[1],vx[3]) ) { + return 0; + } /* Modify 4th tetra */ pt0->v[0] = vx[2]; pt0->v[1] = vx[4]; pt0->v[2] = vx[5]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normals at triangles vx[0],vx[1] and + * vx[2] */ + if ( !MMG3D_check_ridpt_normal_proj_3fac(mesh,pt0,0,vx[0],vx[1],vx[2]) ) { + return 0; + } /* Modify 5th tetra */ pt0->v[0] = vx[0]; pt0->v[1] = vx[3]; pt0->v[2] = vx[1]; pt0->v[3] = vx[2]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normal at triangle vx[0]-vx[1]-vx[2] */ + if ( !MMG3D_check_ridpt_normal_proj_fac(mesh,pt0,0,vx[3]) ) { + return 0; + } /* Modify 6th tetra */ pt0->v[0] = vx[2]; pt0->v[1] = vx[0]; pt0->v[2] = vx[3]; pt0->v[3] = vx[4]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normal at triangle vx[0]-vx[1]-vx[3] */ + if ( !MMG3D_check_ridpt_normal_proj_fac(mesh,pt0,0,vx[2]) ) { + return 0; + } /* Modify 7th tetra */ pt0->v[0] = vx[2]; pt0->v[1] = vx[3]; pt0->v[2] = vx[1]; pt0->v[3] = vx[5]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normal at triangle vx[0]-vx[2]-vx[3] */ + if ( !MMG3D_check_ridpt_normal_proj_fac(mesh,pt0,0,vx[1]) ) { + return 0; + } /* Modify last tetra */ pt0->v[0] = vx[2]; pt0->v[1] = vx[3]; pt0->v[2] = vx[5]; pt0->v[3] = vx[4]; vnew = MMG5_orvol(mesh->point,pt0->v); if ( vnew < MMG5_EPSOK ) return 0; + /* Check the projection of the normals at ridge points onto the normal at + * the new triangle (otherwise we will fail later in mmg3dBezierCP): for + * this tetra we want to check normal at triangle vx[1]-vx[2]-vx[3] */ + if ( !MMG3D_check_ridpt_normal_proj_fac(mesh,pt0,0,vx[0]) ) { + return 0; + } return 1; } From b3fc68f5b9ce463fd7c07b11084a2e0fb532c8b5 Mon Sep 17 00:00:00 2001 From: Algiane Froehly Date: Wed, 23 Nov 2022 19:16:02 +0100 Subject: [PATCH 5/5] Fix wrong normal check in split2sf. --- src/mmg3d/split_3d.c | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/mmg3d/split_3d.c b/src/mmg3d/split_3d.c index 0c9bda535..90f59d438 100644 --- a/src/mmg3d/split_3d.c +++ b/src/mmg3d/split_3d.c @@ -1271,9 +1271,9 @@ int MMG3D_split2sf_sim(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int vx[6]){ /* Check the projection of the normals at ridge points onto the normal at * the new triangle (otherwise we will fail later in mmg3dBezierCP): for * this tetra we want to check normals at triangles tau[0]-tau[1]-tau[3] and - * tau[0]-tau[1]-tau[5], that are, the two faces on both sides of edge - * taued[0] */ - if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[0]) ) { + * tau[1]-tau[3]-tau[2], that are, the two faces on both sides of edge + * taued[4] */ + if ( !MMG3D_check_ridpt_normal_proj_edg(mesh,pt0,0,taued[4]) ) { return 0; }