From 23d3673ef27547d18dc6ff5bffe29c398decc9a9 Mon Sep 17 00:00:00 2001 From: Algiane Date: Wed, 13 Nov 2019 16:39:51 +0100 Subject: [PATCH 01/10] Copy the mmgHashTria function into hash_s.c (in order to modify the adjacency relations. --- src/mmgs/hash_s.c | 144 +++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 143 insertions(+), 1 deletion(-) diff --git a/src/mmgs/hash_s.c b/src/mmgs/hash_s.c index 9d6a2ae85..ddf0d1747 100644 --- a/src/mmgs/hash_s.c +++ b/src/mmgs/hash_s.c @@ -75,7 +75,12 @@ static int paktri(MMG5_pMesh mesh) { * */ int MMGS_hashTria(MMG5_pMesh mesh) { - MMG5_Hash hash; + MMG5_Hash hash; + MMG5_pTria pt,pt1; + MMG5_hedge *ph; + int *adja,*adjt,k,kk,jel,lel,hmax,dup,nmf,ia,ib; + char i,i1,i2,j,l; + unsigned int key; if ( mesh->adja ) return 1; if ( abs(mesh->info.imprim) > 5 || mesh->info.ddebug ) @@ -91,6 +96,143 @@ int MMGS_hashTria(MMG5_pMesh mesh) { if ( !MMG5_mmgHashTria(mesh, mesh->adja, &hash, 0) ) return 0; + adjt = mesh->adja; + + /* adjust hash table params */ + hmax = 3.71*mesh->np; + hash.siz = mesh->np; + hash.max = hmax + 1; + hash.nxt = hash.siz; + MMG5_ADD_MEM(mesh,(hash.max+1)*sizeof(MMG5_hedge),"hash table",return 0); + MMG5_SAFE_CALLOC(hash.item,hash.max+1,MMG5_hedge,return 0); + + for (k=hash.siz; kinfo.ddebug ) fprintf(stdout," h- stage 1: init\n"); + + /* hash triangles */ + mesh->base = 1; + dup = nmf = 0; + for (k=1; k<=mesh->nt; k++) { + pt = &mesh->tria[k]; + if ( !MG_EOK(pt) ) continue; + + pt->flag = 0; + pt->base = mesh->base; + adja = &adjt[3*(k-1)+1]; + for (i=0; i<3; i++) { + i1 = MMG5_inxt2[i]; + i2 = MMG5_iprv2[i]; + + /* compute key */ + ia = MG_MIN(pt->v[i1],pt->v[i2]); + ib = MG_MAX(pt->v[i1],pt->v[i2]); + key = (MMG5_KA*ia + MMG5_KB*ib) % hash.siz; + ph = &hash.item[key]; + + /* store edge */ + if ( ph->a == 0 ) { + ph->a = ia; + ph->b = ib; + ph->k = 3*k + i; + ph->nxt = 0; + ++ph->s; + continue; + } + /* update info about adjacent */ +#ifndef NDEBUG + char ok = 0; +#endif + while ( ph->a ) { + if ( ph->a == ia && ph->b == ib ) { + jel = ph->k / 3; + j = ph->k % 3; + pt1 = &mesh->tria[jel]; + /* discard duplicate face */ + if ( pt1->v[j] == pt->v[i] ) { + pt1->v[0] = 0; + dup++; + } + /* update adjacent */ + else if ( !adjt[3*(jel-1)+1+j] ) { + adja[i] = 3*jel + j; + adjt[3*(jel-1)+1+j] = 3*k + i; + ++ph->s; + } + /* non-manifold case */ + else if ( adja[i] != 3*jel+j ) { + lel = adjt[3*(jel-1)+1+j]/3; + l = adjt[3*(jel-1)+1+j]%3; + (mesh->tria[lel]).tag[l] |= MG_GEO + MG_NOM; + pt1->tag[j] |= MG_GEO + MG_NOM; + pt->tag[i] |= MG_GEO + MG_NOM; + nmf++; + ++ph->s; + } +#ifndef NDEBUG + ok = 1; +#endif + break; + } + else if ( !ph->nxt ) { + ph->nxt = hash.nxt; + ph = &hash.item[ph->nxt]; + assert(ph); + + if ( hash.nxt >= hash.max-1 ) { + if ( mesh->info.ddebug ) { + fprintf(stderr,"\n ## Warning: %s: memory alloc problem (edge):" + " %d\n",__func__,hash.max); + } + MMG5_TAB_RECALLOC(mesh,hash.item,hash.max,MMG5_GAP,MMG5_hedge, + "MMG5_edge", + MMG5_DEL_MEM(mesh,hash.item); + return 0); + + ph = &hash.item[hash.nxt]; + + for (kk=ph->nxt; kknxt; + ph->a = ia; + ph->b = ib; + ph->k = 3*k + i; + ph->nxt = 0; + ++ph->s; +#ifndef NDEBUG + ok = 1; +#endif + break; + } + else + ph = &hash.item[ph->nxt]; + } + assert(ok); + } + } + + /* set tag */ + for (k=1; k<=mesh->nt; k++) { + pt = &mesh->tria[k]; + for (i=0; i<3; i++) { + if ( pt->tag[i] & MG_NOM ) { + mesh->point[pt->v[MMG5_inxt2[i]]].tag |= MG_NOM; + mesh->point[pt->v[MMG5_iprv2[i]]].tag |= MG_NOM; + } + } + } + + if ( abs(mesh->info.imprim) > 5 && dup+nmf > 0 ) { + fprintf(stdout," ## "); fflush(stdout); + if ( nmf > 0 ) fprintf(stdout,"[non-manifold model] "); + if ( dup > 0 ) fprintf(stdout," %d duplicate removed",dup); + fprintf(stdout,"\n"); + } + if ( mesh->info.ddebug ) fprintf(stdout," h- completed.\n"); + MMG5_DEL_MEM(mesh,hash.item); return 1; From 563ec4dbfb0ea1d25596fe3d08fc7778e25e65fe Mon Sep 17 00:00:00 2001 From: Algiane Date: Tue, 19 Nov 2019 11:09:56 +0100 Subject: [PATCH 02/10] forbid the call of bouletrid on nm points. --- src/mmgs/boulep_s.c | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/mmgs/boulep_s.c b/src/mmgs/boulep_s.c index 15cb5ca1b..54b467ea1 100644 --- a/src/mmgs/boulep_s.c +++ b/src/mmgs/boulep_s.c @@ -256,6 +256,7 @@ int boulechknm(MMG5_pMesh mesh,int start,int ip,int *list) { * to normal \a n1's side. \a ip0 and \a ip1 are the indices of the 2 ending * point of the ridge. Both lists are returned enumerated in direct order. * + * \warning can't be called from a non-manifold point */ int bouletrid(MMG5_pMesh mesh,int start,int ip,int *il1,int *l1,int *il2,int *l2,int *ip0,int *ip1) { MMG5_pTria pt; @@ -269,7 +270,9 @@ int bouletrid(MMG5_pMesh mesh,int start,int ip,int *il1,int *l1,int *il2,int *l2 idp = pt->v[ip]; ppt = &mesh->point[idp]; + assert( ppt->tag & MG_GEO ); + assert( !(ppt->tag & MG_NOM) ); /* set pointers: first manifold is on side of triangle */ if ( !MMG5_nortri(mesh,pt,nt) ) return 0; From d0a31c9347c2256fb1f129b2bc31f3788f090944 Mon Sep 17 00:00:00 2001 From: Algiane Date: Tue, 19 Nov 2019 11:11:01 +0100 Subject: [PATCH 03/10] DOxygen corrections + comment. --- src/common/libmmgtypes.h | 2 +- src/mmgs/analys_s.c | 2 +- src/mmgs/colver_s.c | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/common/libmmgtypes.h b/src/common/libmmgtypes.h index 5fae7e3e1..c80e50ffe 100644 --- a/src/common/libmmgtypes.h +++ b/src/common/libmmgtypes.h @@ -427,7 +427,7 @@ typedef struct { typedef MMG5_xPrism * MMG5_pxPrism; /** - * \struc MMG5_Mat + * \struct MMG5_Mat * \brief To store user-defined references in the mesh (useful in LS mode) */ typedef struct { diff --git a/src/mmgs/analys_s.c b/src/mmgs/analys_s.c index 6b7e3b4e9..925ed8d0f 100644 --- a/src/mmgs/analys_s.c +++ b/src/mmgs/analys_s.c @@ -520,7 +520,7 @@ static int MMG5_singul(MMG5_pMesh mesh) { } } - /* check for handle */ + /* check for handle (point connecting two surfaces non connected by edges) */ for (k=1; k<=mesh->nt; k++) { pt = &mesh->tria[k]; if ( !MG_EOK(pt) ) continue; diff --git a/src/mmgs/colver_s.c b/src/mmgs/colver_s.c index 3d5d5b6b0..d55cee155 100644 --- a/src/mmgs/colver_s.c +++ b/src/mmgs/colver_s.c @@ -345,7 +345,7 @@ int colver(MMG5_pMesh mesh,int *list,int ilist) { * \param list pointer toward the ball of the point to collapse. * \return 1 if success, 0 if fail. * - * Collapse edge \f$list[0]%3\f$ in tet \f$list[0]/3\f$ (\f$ ip->i1\f$ ) for a + * Collapse edge \f$list[0]\%3\f$ in tet \f$list[0]/3\f$ (\f$ ip->i1\f$ ) for a * ball of the collapsed point of size 3: the collapsed point is removed. * */ From fe417740c4d08bf9af5f160d601678a9ac75f2b7 Mon Sep 17 00:00:00 2001 From: Algiane Date: Thu, 21 Nov 2019 10:42:41 +0100 Subject: [PATCH 04/10] Add circular adjacency to allow the management of non manifold surface meshes in mmgs. For now the non manifold points remains freezed. --- src/mmgs/analys_s.c | 45 ++++++++++++- src/mmgs/boulep_s.c | 159 ++++++++++++++++++-------------------------- src/mmgs/chkmsh_s.c | 149 ++++++++++++++++++++++++++++++++++------- src/mmgs/hash_s.c | 4 ++ src/mmgs/libmmgs.c | 51 ++++++-------- src/mmgs/mmgs.h | 2 + src/mmgs/mmgs1.c | 16 ++++- 7 files changed, 269 insertions(+), 157 deletions(-) diff --git a/src/mmgs/analys_s.c b/src/mmgs/analys_s.c index 925ed8d0f..6a50f8e89 100644 --- a/src/mmgs/analys_s.c +++ b/src/mmgs/analys_s.c @@ -35,6 +35,42 @@ #include "mmgs.h" +/** + * \param mesh pointer toward the mesh structure + * \param k starting tria + * \param edge in the tria + * \return 1 if k is the lowest tria that contains the edge i, 0 otherwise + * + * Travel by adjacency through the list of trias that contains the edge \a i of + * \a k and return 1 if k is the tria of lowest index. + * + */ +int MMGS_is_smallestAdja(MMG5_pMesh mesh, int k, int i) { + int kmin = k; + int imin = i; + int kcur = k; + int icur = i; + int next = mesh->adja[ 3*(kcur-1) + icur+1]; + + if ( next ) { + while ( next != 3*k + i ) { + kcur = next/3; + icur = next%3; + next = mesh->adja[ 3*(kcur-1) + icur+1 ]; + + if ( kcur < kmin ) { + kmin = kcur; + imin = icur; + } + } + } + if ( kmin != k ) { + return 0; + } + + return 1; +} + /** * \param mesh pointer toward the mesh * @@ -96,6 +132,9 @@ static int setadj(MMG5_pMesh mesh){ kk = adja[i] / 3; ii = adja[i] % 3; + + /* Due to the order in which adjacency relationships are created + * (1->4->3->2->1), nm edges are counted only once */ if ( kk > k ) ned++; /* store adjacent */ @@ -205,8 +244,7 @@ static int setadj(MMG5_pMesh mesh){ for (i=0; i<3; i++) { if ( ( !MG_EDG(pt->tag[i]) ) && ( !(pt->tag[i] & MG_REQ) ) ) continue; - jel = adja[i] / 3; - if ( !jel || jel > k ) { + if ( MMGS_is_smallestAdja(mesh, k,i) ) { if ( pt->tag[i] & MG_GEO ) nr++; if ( pt->tag[i] & MG_REF ) nre++; if ( pt->tag[i] & MG_REQ ) nreq++; @@ -239,7 +277,7 @@ static void nmpoints(MMG5_pMesh mesh) { char i0,i1,i,jp; nmp = 0; - /* Initialize point flags */ + /* Initialize point flags with the index of a tria that contains the point */ for (k=1; k<=mesh->np; k++) mesh->point[k].s = 0; @@ -392,6 +430,7 @@ static void nmpoints(MMG5_pMesh mesh) { /* break; */ /* } */ /* else if ( adja[i] ) { */ +/* #warning does it work on nm points? */ /* iel = adja[i] / 3; */ /* j = adja[i] % 3; */ /* if ( litcol(mesh,iel,j,kal) ) { */ diff --git a/src/mmgs/boulep_s.c b/src/mmgs/boulep_s.c index 54b467ea1..8934ea30d 100644 --- a/src/mmgs/boulep_s.c +++ b/src/mmgs/boulep_s.c @@ -40,59 +40,66 @@ * \param start index of triangle to start. * \param ip index of point for wich we compute the ball. * \param list pointer toward the computed ball of \a ip. - * \return the size of the computed ball or 0 if fail. + * \return the size of the computed ball (ilist) or -ilist if fail. * * Find all triangles sharing \a ip, \f$list[0] =\f$ \a start do not stop when * crossing ridge. * */ int boulet(MMG5_pMesh mesh,int start,int ip,int *list) { - MMG5_pTria pt; + MMG5_pTria pt,pt1; MMG5_pPoint ppt; - int *adja,k,ilist; - char i,i1,i2; + int *adja,cur,k,ilist,nump,k1,j,base; + char i,l; - pt = &mesh->tria[start]; + base = ++mesh->base; + pt = &mesh->tria[start]; - ppt = &mesh->point[pt->v[ip]]; - ilist = 0; + nump = pt->v[ip]; + ppt = &mesh->point[nump]; - /* store neighbors */ - k = start; - i = ip; - do { - if ( ilist > MMGS_LMAX-2 ) return 0; - list[ilist] = 3*k + i; - ++ilist; + /* store initial tria */ + pt->flag = base; + list[0] = 3*start + ip; + ilist = 1; - adja = &mesh->adja[3*(k-1)+1]; - i1 = MMG5_inxt2[i]; - k = adja[i1] / 3; - i = adja[i1] % 3; - i = MMG5_inxt2[i]; - } - while ( k && k != start ); - if ( k > 0 ) return ilist; + /* Explore list and travel by adjacency through elements sharing nump */ + cur = 0; - if ( ppt->tag & MG_NOM ) - return 0; + while ( cur < ilist ) { + k = list[cur] / 3; + i = list[cur] % 3; // index of point nump in tria k - /* check if boundary hit */ - k = start; - i = ip; - do { - adja = &mesh->adja[3*(k-1)+1]; - i2 = MMG5_iprv2[i]; - k = adja[i2] / 3; - if ( k == 0 ) break; - i = adja[i2] % 3; - i = MMG5_iprv2[i]; + adja = &mesh->adja[ 3*(k-1)+ 1 ]; + + for (l=0; l<2; l++) { + i = MMG5_inxt2[i]; + k1 = adja[i]; + if ( !k1 ) continue; + + j = k1%3; + k1 /= 3; + + pt1 = &mesh->tria[k1]; + if ( pt1->flag == base ) continue; + pt1->flag = base; + + /* Find the local index of nump */ + j = MMG5_inxt2[j]; + if ( pt1->v[MMG5_inxt2[j]] == nump ) { + j = MMG5_inxt2[j]; + } + assert ( pt1->v[j] == nump ); + + /* Check for overflow */ + if ( ilist > MMGS_LMAX-2 ) return 0; - if ( ilist > MMGS_LMAX-2 ) return 0; - list[ilist] = 3*k + i; - ilist++; + /* Store adjacent */ + list[ilist] = 3*k1 + j; + ++ilist; + } + ++cur; } - while ( k ); return ilist; } @@ -104,6 +111,9 @@ int boulet(MMG5_pMesh mesh,int start,int ip,int *list) { * the ball. * \param list pointer toward the computed ball of point. * + * \return -ilist (number of tria in ball) if buffer overflow, ilist if the + * collapse is ok, 0 if the collapse lead to a non manifold situation. + * * Find all triangles sharing \a ip, \f$list[0] = start\f$. Do not stop when * crossing ridge. Check whether resulting configuration is manifold. * @@ -114,70 +124,26 @@ int boulechknm(MMG5_pMesh mesh,int start,int ip,int *list) { int *adja,k,ilist,base,iel; char i,i1,i2,ia,iq,voy; - base = ++mesh->base; - pt = &mesh->tria[start]; - ia = MMG5_iprv2[ip]; - iq = MMG5_inxt2[ip]; if ( !MG_EOK(pt) ) return 0; - ppt = &mesh->point[pt->v[ip]]; - if ( ppt->tag & MG_NOM ) return 0; - ilist = 0; - /* store neighbors */ - k = start; - i = ip; - do { - if ( ilist > MMGS_LMAX-2 ) return -ilist; - list[ilist] = 3*k + i; - ++ilist; - pt = &mesh->tria[k]; +#warning ajeter : on, fait quoi dans ce cas? + if ( mesh->point[pt->v[ip]].tag & MG_NOM ) return 0; - i1 = MMG5_inxt2[i]; - i2 = MMG5_iprv2[i]; - ppt = &mesh->point[pt->v[i1]]; - ppt->s = base; - ppt = &mesh->point[pt->v[i2]]; - ppt->s = base; + /** Stage 1: fill the ball of point */ + ilist = boulet(mesh,start,ip,list); - adja = &mesh->adja[3*(k-1)+1]; - k = adja[i1] / 3; - i = adja[i1] % 3; - i = MMG5_inxt2[i]; - } - while ( k && k != start ); - - /* check if boundary hit */ - if ( k <= 0 ) { - k = start; - i = ip; - do { - adja = &mesh->adja[3*(k-1)+1]; - i1 = MMG5_inxt2[i]; - i2 = MMG5_iprv2[i]; - - pt = &mesh->tria[k]; - ppt = &mesh->point[pt->v[i1]]; - ppt->s = base; - ppt = &mesh->point[pt->v[i2]]; - ppt->s = base; - - k = adja[i2] / 3; - if ( k == 0 ) break; - i = adja[i2] % 3; - i = MMG5_iprv2[i]; - - if ( ilist > MMGS_LMAX-2 ) return -ilist; - list[ilist] = 3*k + i; - ilist++; - } - while ( k ); - } + if ( ilist < 1 ) return ilist; + /** Stage 2: check if a collapse may lead to a non-convex situation */ + /* Flags initialization */ pt = &mesh->tria[start]; - i1 = MMG5_inxt2[ip]; - i2 = MMG5_iprv2[ip]; + ia = MMG5_iprv2[ip]; + iq = MMG5_inxt2[ip]; + + i1 = ia; + i2 = iq; ppt = &mesh->point[pt->v[i1]]; ppt->s = 0; ppt = &mesh->point[pt->v[i2]]; @@ -193,7 +159,8 @@ int boulechknm(MMG5_pMesh mesh,int start,int ip,int *list) { ppt->s = 0; } - /* check if a collapse may lead to a non-convex situation */ + /* Unfold the ball af iq and check that it doesn't contains common tria with + * the ball of ip (except for the ip-iq shell) */ k = start; i = iq; do { @@ -214,7 +181,7 @@ int boulechknm(MMG5_pMesh mesh,int start,int ip,int *list) { while ( k && k != start ); if( k > 0 ) return ilist; - /* check if boundary hit */ + /* Unfold the ball in opposite direction if boundary hit */ k = start; i = iq; do { diff --git a/src/mmgs/chkmsh_s.c b/src/mmgs/chkmsh_s.c index 111f484a4..9356dc7d1 100644 --- a/src/mmgs/chkmsh_s.c +++ b/src/mmgs/chkmsh_s.c @@ -35,6 +35,99 @@ #include "mmgs.h" +/** + * \param mesh pointer toward the mesh structure. + * \return 0 if fail, 1 if success. + * + * Check the adjacency relationship validity + * + */ +int MMGS_chkadja(MMG5_pMesh mesh) { + MMG5_pTria pt,pt1,pt2; + int k,i,key,jel,j,jkey,lel,l,kcur,icur,next,count; + +#define MMGS_MAXNM 20 + + if ( !mesh->adja ) return 1; + + for ( k=1; k<=mesh->nt; ++k ) { + pt = &mesh->tria[k]; + + for ( i=0; i<3; ++i) { + if ( pt->tag[i] & MG_NOM ) { + /* Check point tag */ + if ( ! (mesh->point[pt->v[MMG5_inxt2[i]]].tag & MG_NOM) ) { + fprintf(stderr," ## Error: %s: nm edge: wrong tag at point %d\n", + __func__,pt->v[MMG5_inxt2[i]]); + } + /* Check point tag */ + if ( ! (mesh->point[pt->v[MMG5_iprv2[i]]].tag & MG_NOM) ) { + fprintf(stderr," ## Error: %s: nm edge: wrong tag at point %d\n", + __func__,pt->v[MMG5_iprv2[i]]); + } + + /* Check circulary adjacency relation along nm edge */ + kcur = k; + icur = i; + next = mesh->adja[ 3*(kcur-1) + icur+1]; + + count = 0; + assert ( next ); + while ( next != 3*k + i && ++count < MMGS_MAXNM ) { + kcur = next/3; + icur = next%3; + next = mesh->adja[ 3*(kcur-1) + icur+1 ]; + + pt1 = &mesh->tria[kcur]; + if ( (pt->v[MMG5_inxt2[i]] != pt1->v[MMG5_iprv2[icur]] || + pt->v[MMG5_iprv2[i]] != pt1->v[MMG5_inxt2[icur]]) && + (pt->v[MMG5_iprv2[i]] != pt1->v[MMG5_iprv2[icur]] || + pt->v[MMG5_inxt2[i]] != pt1->v[MMG5_inxt2[icur]]) ) { + fprintf(stderr," ## Error: %s: nm edge: wrong adjacency relationship:\n" + " elt %d: %d %d %d (edg %d, tag %d) adjacent to elt %d: %d %d %d (edg %d, tag %d)\n", + __func__,k,pt->v[0],pt->v[1],pt->v[2],i,pt->tag[i],kcur,pt1->v[0],pt1->v[1],pt1->v[2],icur,pt1->tag[icur]); + return 0; + } + if ( !(pt1->tag[icur] & MG_NOM) ) { + fprintf(stderr," ## Error: %s: non symmetric tag on nm edge:\n" + " elt %d: %d %d %d (edg %d, tag %d) adjacent to elt %d: %d %d %d (edg %d, tag %d)\n", + __func__,k,pt->v[0],pt->v[1],pt->v[2],i,pt->tag[i],kcur,pt1->v[0],pt1->v[1],pt1->v[2],icur,pt1->tag[icur]); + } + } + } + else { + /* Classic edge */ + key = mesh->adja[3*(k-1)+i+1]; + + if ( !key ) continue; + + jel = key/3; + j = key%3; + + jkey = mesh->adja[3*(jel-1) + j+1]; + lel = jkey/3; + l = jkey%3; + + pt1 = &mesh->tria[jel]; + if ( (3*k+i != jkey) || + ( pt->v[MMG5_inxt2[i]] != pt1->v[MMG5_iprv2[j]] || + pt->v[MMG5_iprv2[i]] != pt1->v[MMG5_inxt2[j]] ) ) { + + pt2 = &mesh->tria[lel]; + fprintf(stderr," ## Error: %s: wrong adjacency relationship:\n" + " elt %d: %d %d %d (edg %d, tag %d) adjacent to elt %d: %d %d %d (edg %d, tag %d)\n" + " elt %d: %d %d %d (edg %d, tag %d) adjacent to elt %d: %d %d %d (edg %d, tag %d)\n", + __func__,k,pt->v[0],pt->v[1],pt->v[2],i,pt->tag[i],jel,pt1->v[0],pt1->v[1],pt1->v[2],j,pt1->tag[j], + jel,pt1->v[0],pt1->v[1],pt1->v[2],j,pt1->tag[j],lel,pt2->v[0],pt2->v[1],pt2->v[2],l,pt2->tag[l]); + return 0; + } + } + + } + } + + return 1; +} /** * \param mesh pointer toward the mesh structure. @@ -54,6 +147,9 @@ int MMG5_mmgsChkmsh(MMG5_pMesh mesh,int severe,int base) { static char mmgErr0=0,mmgErr1=0,mmgErr2=0,mmgErr3=0,mmgErr4=0; static char mmgErr5=0,mmgErr6=0,mmgErr7=0; + if ( !mesh->adja ) return 1; + MMGS_chkadja(mesh); + for (k=1; k<=mesh->nt; k++) { pt1 = &mesh->tria[k]; if ( !MG_EOK(pt1) ) continue; @@ -117,32 +213,34 @@ int MMG5_mmgsChkmsh(MMG5_pMesh mesh,int severe,int base) { MMGS_indElt(mesh,adj),pt1->tag[i],pt2->tag[voy]); return 0; } - adjb = &mesh->adja[3*(adj-1)+1]; - adj1 = adjb[voy] / 3; - voy1 = adjb[voy] % 3; - if ( adj1 != k || voy1 != i ) { - if ( !mmgErr3 ) { - mmgErr3 = 1; - fprintf(stderr,"\n ## Error: %s: 2. at least 1 wrong" - " adjacency (%d %d)\n",__func__, MMGS_indElt(mesh,k), - MMGS_indElt(mesh,adj1)); - fprintf(stderr,"vertices of %d: %d %d %d \n",MMGS_indElt(mesh,k), - MMGS_indPt(mesh,pt1->v[0]),MMGS_indPt(mesh,pt1->v[1]), - MMGS_indPt(mesh,pt1->v[2])); - fprintf(stderr,"vertices of adj %d: %d %d %d \n", - MMGS_indElt(mesh,adj), - MMGS_indPt(mesh,pt2->v[0]),MMGS_indPt(mesh,pt2->v[1]), - MMGS_indPt(mesh,pt2->v[2])); - fprintf(stderr,"adj(%d): %d %d %d\n",MMGS_indElt(mesh,k), - MMGS_indElt(mesh,adja[0]/3),MMGS_indElt(mesh,adja[1]/3), - MMGS_indElt(mesh,adja[2]/3)); - fprintf(stderr,"adj(%d): %d %d %d\n",MMGS_indElt(mesh,adj), - MMGS_indElt(mesh,adjb[0]/3),MMGS_indElt(mesh,adjb[1]/3), - MMGS_indElt(mesh,adjb[2]/3)); + + if ( !(pt1->tag[i] & MG_NOM) ) { + adjb = &mesh->adja[3*(adj-1)+1]; + adj1 = adjb[voy] / 3; + voy1 = adjb[voy] % 3; + if ( adj1 != k || voy1 != i ) { + if ( !mmgErr3 ) { + mmgErr3 = 1; + fprintf(stderr,"\n ## Error: %s: 2. at least 1 wrong" + " adjacency (%d %d)\n",__func__, MMGS_indElt(mesh,k), + MMGS_indElt(mesh,adj1)); + fprintf(stderr,"vertices of %d: %d %d %d \n",MMGS_indElt(mesh,k), + MMGS_indPt(mesh,pt1->v[0]),MMGS_indPt(mesh,pt1->v[1]), + MMGS_indPt(mesh,pt1->v[2])); + fprintf(stderr,"vertices of adj %d: %d %d %d \n", + MMGS_indElt(mesh,adj), + MMGS_indPt(mesh,pt2->v[0]),MMGS_indPt(mesh,pt2->v[1]), + MMGS_indPt(mesh,pt2->v[2])); + fprintf(stderr,"adj(%d): %d %d %d\n",MMGS_indElt(mesh,k), + MMGS_indElt(mesh,adja[0]/3),MMGS_indElt(mesh,adja[1]/3), + MMGS_indElt(mesh,adja[2]/3)); + fprintf(stderr,"adj(%d): %d %d %d\n",MMGS_indElt(mesh,adj), + MMGS_indElt(mesh,adjb[0]/3),MMGS_indElt(mesh,adjb[1]/3), + MMGS_indElt(mesh,adjb[2]/3)); + } + return 0; } - return 0; - } - if ( !MS_SIN(pt1->tag[i]) ) { + if ( !MS_SIN(pt1->tag[i]) ) { j1 = MMG5_inxt2[voy]; j2 = MMG5_iprv2[voy]; if ( pt2->v[j2] != pt1->v[i1] || pt2->v[j1] != pt1->v[i2] ) { @@ -154,6 +252,7 @@ int MMG5_mmgsChkmsh(MMG5_pMesh mesh,int severe,int base) { } return 0; } + } } } } diff --git a/src/mmgs/hash_s.c b/src/mmgs/hash_s.c index ddf0d1747..db062abb6 100644 --- a/src/mmgs/hash_s.c +++ b/src/mmgs/hash_s.c @@ -164,6 +164,10 @@ int MMGS_hashTria(MMG5_pMesh mesh) { else if ( adja[i] != 3*jel+j ) { lel = adjt[3*(jel-1)+1+j]/3; l = adjt[3*(jel-1)+1+j]%3; + + /* Circle adjacency relation jel -> k -> lel -> ... -> jel */ + adjt[3*(jel-1)+1+j] = 3*k + i; + adjt[3*(k -1)+1+i] = 3*lel + l; (mesh->tria[lel]).tag[l] |= MG_GEO + MG_NOM; pt1->tag[j] |= MG_GEO + MG_NOM; pt->tag[i] |= MG_GEO + MG_NOM; diff --git a/src/mmgs/libmmgs.c b/src/mmgs/libmmgs.c index 4d0a3dbf1..8a06d19b3 100644 --- a/src/mmgs/libmmgs.c +++ b/src/mmgs/libmmgs.c @@ -90,7 +90,7 @@ static inline int MMGS_packMesh(MMG5_pMesh mesh,MMG5_pSol sol,MMG5_pSol met) { MMG5_pTria pt,ptnew; MMG5_pPoint ppt,pptnew; - int np,nc,nr, k,nt,nbl,imet,imetnew,i,na,jel; + int np,nc,nr, k,nt,nbl,imet,imetnew,i,na,jel,kcur,icur,next; int iadr,iadrnew,iadrv,*adjav,*adja,*adjanew,voy; char i1,i2; @@ -127,36 +127,33 @@ int MMGS_packMesh(MMG5_pMesh mesh,MMG5_pSol sol,MMG5_pSol met) { for(i=0 ; i<3 ; i++) { adjanew[i] = adja[i]; if(!adja[i]) continue; - iadrv = 3*(adja[i]/3-1) +1; + + next = adja[i]; + while ( next != 3*k + i ) { + kcur = next/3; + icur = next%3; + next = mesh->adja[ 3*(kcur-1) + icur+1 ]; + } + + iadrv = 3*(kcur-1) +1; adjav = &mesh->adja[iadrv]; voy = i; - adjav[adja[i]%3] = 3*nbl + voy; + adjav[icur] = 3*nbl + voy; } } - nbl++; /* Count the edges */ for(i=0 ; i<3 ; i++) { if ( !MG_EDG(pt->tag[i]) ) continue; - adja = &mesh->adja[3*(k-1)+1]; - jel = adja[i] / 3; - - assert ( jel != k ); - if ( jel ) { - /* If we have an adjacent, treat the edge either from the tria with - * higher ref or, if both tria have the same references, from the tria - * with higher index */ - if ( mesh->tria[jel].ref < mesh->tria[k].ref ) { - continue; - } - else if ( (mesh->tria[jel].ref == mesh->tria[k].ref) && (k < jel) ) { - continue; - } + /* Count the edge only if we come from the lowest tria of the list of tria + * that contains this edge */ + if ( MMGS_is_smallestAdja(mesh, nbl,i) ) { + ++na; } - ++na; } + nbl++; } mesh->nt = nt; @@ -234,18 +231,10 @@ int MMGS_packMesh(MMG5_pMesh mesh,MMG5_pSol sol,MMG5_pSol met) { for (i=0; i<3; i++) { if ( !MG_EDG(pt->tag[i]) ) continue; - adja = &mesh->adja[3*(k-1)+1]; - jel = adja[i] / 3; - if ( jel ) { - /* If we have an adjacent, treat the edge either from the tria with - * higher ref or, if both tria have the same references, from the tria - * with higher index */ - if ( mesh->tria[jel].ref < mesh->tria[k].ref ) { - continue; - } - else if ( (mesh->tria[jel].ref == mesh->tria[k].ref) && (k < jel) ) { - continue; - } + /* Treat the edge only if we come from the lowest tria of the list of + * tria that contains this edge */ + if ( !MMGS_is_smallestAdja(mesh, k,i) ) { + continue; } i1 = MMG5_inxt2[i]; diff --git a/src/mmgs/mmgs.h b/src/mmgs/mmgs.h index 9f9346dfc..4539fcb16 100644 --- a/src/mmgs/mmgs.h +++ b/src/mmgs/mmgs.h @@ -126,6 +126,7 @@ int MMGS_Free_names_var( va_list argptr ); int MMGS_zaldy(MMG5_pMesh mesh); int MMGS_assignEdge(MMG5_pMesh mesh); +int MMGS_is_smallestAdja(MMG5_pMesh mesh, int k, int i); int MMGS_analys(MMG5_pMesh mesh); int MMGS_inqua(MMG5_pMesh,MMG5_pSol); int MMGS_outqua(MMG5_pMesh,MMG5_pSol); @@ -163,6 +164,7 @@ char typelt(MMG5_pPoint p[3],char *ia); int litswp(MMG5_pMesh mesh,int k,char i,double kal); int litcol(MMG5_pMesh mesh,int k,char i,double kal); int MMG5_mmgsChkmsh(MMG5_pMesh,int,int); +int MMGS_chkadja(MMG5_pMesh mesh); int paratmet(double c0[3],double n0[3],double m[6],double c1[3],double n1[3],double mt[6]); int intregmet(MMG5_pMesh mesh,MMG5_pSol met,int k,char i,double s,double mr[6]); int MMG5_intridmet(MMG5_pMesh,MMG5_pSol,int,int,double,double*,double*); diff --git a/src/mmgs/mmgs1.c b/src/mmgs/mmgs1.c index cbb920556..b67d7d528 100755 --- a/src/mmgs/mmgs1.c +++ b/src/mmgs/mmgs1.c @@ -971,6 +971,11 @@ static int colelt(MMG5_pMesh mesh,MMG5_pSol met,char typchk) { for (i=0; i<3; i++) { if ( MS_SIN(pt->tag[i]) ) continue; +#ifndef NDEBUG + if ( mesh->adja[3*(k-1)+1+i] ) + assert ( 3*k+i == mesh->adja[ 3* (mesh->adja[3*(k-1)+1+i]/3-1) + 1 +mesh->adja[3*(k-1)+1+i]%3 ] ); +#endif + i1 = MMG5_inxt2[i]; i2 = MMG5_iprv2[i]; p1 = &mesh->point[pt->v[i1]]; @@ -1154,12 +1159,19 @@ static int adpcol(MMG5_pMesh mesh,MMG5_pSol met) { p1 = &mesh->point[pt->v[i1]]; p2 = &mesh->point[pt->v[i2]]; if ( MS_SIN(p1->tag) ) continue; - else if ( p1->tag > p2->tag || p1->tag > pt->tag[i] ) continue; + +#ifndef NDEBUG + if ( mesh->adja[3*(k-1)+1+i] ) { + assert ( 3*k+i == mesh->adja[ 3* (mesh->adja[3*(k-1)+1+i]/3-1) + 1 +mesh->adja[3*(k-1)+1+i]%3 ] ); + } +#endif + + if ( p1->tag > p2->tag || p1->tag > pt->tag[i] ) continue; /* check if geometry preserved */ ilist = chkcol(mesh,met,k,i,list,2); if ( ilist > 3 ) { - ier = colver(mesh,list,ilist);; + ier = colver(mesh,list,ilist); nc += ier; if ( !ier ) return -1; break; From 88fd583d37a82bf52daefce89fcd909ddc217d4d Mon Sep 17 00:00:00 2001 From: Algiane Date: Mon, 2 Dec 2019 11:00:49 +0100 Subject: [PATCH 05/10] New implementation of split1b for nm edges. --- src/mmgs/split_s.c | 247 ++++++++++++++++++++++++++++----------------- 1 file changed, 155 insertions(+), 92 deletions(-) diff --git a/src/mmgs/split_s.c b/src/mmgs/split_s.c index 22783f5f2..64f8b89e7 100644 --- a/src/mmgs/split_s.c +++ b/src/mmgs/split_s.c @@ -182,74 +182,122 @@ int MMGS_simbulgept(MMG5_pMesh mesh,MMG5_pSol met, int k,int i,int ip) { cal = MMG5_calelt(mesh,met,pt0); if ( cal < MMG5_EPSOK ) return 0; - // Check the validity of the two triangles created from the triangle adjacent - // to k by edge i. - kadja = mesh->adja[3*(k-1)+i+1]/3; - iadja = mesh->adja[3*(k-1)+i+1]%3; + // Check the validity of the triangles (2 for manifold test cases, more for + // non manifold ones) created from the triangle adjacent to k by edge i. + kadja = mesh->adja[3*(k-1)+i+1]; if ( !kadja ) return 1; - pt = &mesh->tria[kadja]; - memcpy(pt0,pt,sizeof(MMG5_Tria)); - is = MMG5_iprv2[iadja]; - pt0->v[is] = 0; - cal = MMG5_calelt(mesh,met,pt0); - if ( cal < MMG5_EPSOK ) return 0; + while ( kadja != 3*k + i ) { + iadja = kadja%3; + kadja /= 3; - pt0->v[is] = pt->v[is]; - is = MMG5_inxt2[iadja]; - pt0->v[is] = 0; - cal = MMG5_calelt(mesh,met,pt0); - if ( cal < MMG5_EPSOK ) return 0; + pt = &mesh->tria[kadja]; + memcpy(pt0,pt,sizeof(MMG5_Tria)); + is = MMG5_iprv2[iadja]; + pt0->v[is] = 0; + cal = MMG5_calelt(mesh,met,pt0); + if ( cal < MMG5_EPSOK ) return 0; + + pt0->v[is] = pt->v[is]; + is = MMG5_inxt2[iadja]; + pt0->v[is] = 0; + cal = MMG5_calelt(mesh,met,pt0); + if ( cal < MMG5_EPSOK ) return 0; + + kadja = mesh->adja[ 3*(kadja-1) + iadja+1 ]; + } return 1; } /** * \param mesh pointer toward the mesh structure. - * \param k index of element to split. - * \param i index of edge to split. + * \param kstart index of element to split. + * \param istart index of edge to split. * \param ip index of the new point. * \return 0 if lack of memory, 1 otherwise. * - * Split element \a k along edge \a i, inserting point \a ip and updating - * the adjacency relations. + * Split element \a kstart along edge \a istart, inserting point \a ip and + * updating the adjacency relations. * * \remark do not call this function in non-manifold case */ -int split1b(MMG5_pMesh mesh,int k,char i,int ip) { +int split1b(MMG5_pMesh mesh,int start,char istart,int ip) { MMG5_pTria pt,pt1; MMG5_pPoint ppt; - MMG5_Bezier b; + MMG5_Bezier b; MMG5_pxPoint go; double uv[2],o[3],no[3],to[3]; - int *adja,iel,jel,kel,mel,ier; - char i1,i2,j,j1,j2,m; + int k,adj,*adja,l,iel,jel,kel,mel,ier,next,*newelt,ilist; + char i,i1,i2,j,m; - iel = MMGS_newElt(mesh); - if ( !iel ) { - MMGS_TRIA_REALLOC(mesh,iel,mesh->gap, - MMG5_INCREASE_MEM_MESSAGE(); - return 0); - } - pt = &mesh->tria[k]; - pt->flag = 0; - pt->base = mesh->base; + /** Count the number of elements in the shell of i */ + ilist = 1; - pt1 = &mesh->tria[iel]; - memcpy(pt1,pt,sizeof(MMG5_Tria)); - memcpy(&mesh->adja[3*(iel-1)+1],&mesh->adja[3*(k-1)+1],3*sizeof(int)); + adja = &mesh->adja[3*(start-1)+1]; + next = adja[istart] ; + + if ( next ) { + while ( next != 3*start + istart ) { + jel = next / 3; + j = next % 3; + + ++ilist; + next = mesh->adja[3*(jel-1)+1+j]; + } + } ppt = &mesh->point[ip]; - if ( pt->edg[i] ) ppt->ref = pt->edg[i]; - if ( pt->tag[i] ) ppt->tag = pt->tag[i]; - adja = &mesh->adja[3*(k-1)+1]; - jel = adja[i] / 3; - j = adja[i] % 3; + MMG5_SAFE_CALLOC(newelt,ilist,int,return 0); + + k = start; + i = istart; + for ( l=0; ltria[k]; + pt->flag = 0; + pt->base = mesh->base; + + iel = MMGS_newElt(mesh); + if ( !iel ) { + MMGS_TRIA_REALLOC(mesh,iel,mesh->gap, + fprintf(stderr,"\n ## Error: %s: unable to allocate" + " a new element.\n",__func__); + MMG5_INCREASE_MEM_MESSAGE(); + l--; + for ( ; l>=0 ; --l ) { + if ( !MMGS_delElt(mesh,newelt[l]) ) return -0; + } + return 0); + } + newelt[l] = iel; + + pt1 = &mesh->tria[iel]; + + memcpy(pt1,pt,sizeof(MMG5_Tria)); + memcpy(&mesh->adja[3*(iel-1)+1],&mesh->adja[3*(k-1)+1],3*sizeof(int)); + + if ( pt->edg[i] ) ppt->ref = pt->edg[i]; + if ( pt->tag[i] ) ppt->tag = pt->tag[i]; + + next = mesh->adja[3*(k-1)+1+i]; + k = next / 3; + i = next % 3; + } + + + /** Update normal n2 if need to be (from the first tria of shell) */ + pt = &mesh->tria[start]; + + adja = &mesh->adja[3*(start-1)+1]; + jel = adja[istart] / 3; + j = adja[istart] % 3; /* update normal n2 if need be */ - if ( jel && pt->tag[i] & MG_GEO ) { + if ( pt->tag[i] & MG_GEO ) { ier = MMG5_bezierCP(mesh,&mesh->tria[jel],&b,1); assert(ier); uv[0] = 0.5; @@ -263,58 +311,73 @@ int split1b(MMG5_pMesh mesh,int k,char i,int ip) { memcpy(go->n2,no,3*sizeof(double)); } - /* update two triangles */ - i1 = MMG5_inxt2[i]; - i2 = MMG5_iprv2[i]; - pt->v[i2] = ip; - pt->tag[i1] = MG_NOTAG; - pt->edg[i1] = 0; - pt1->v[i1] = ip; - pt1->tag[i2] = MG_NOTAG; - pt1->edg[i2] = 0; - - /* update adjacency relations */ - mel = adja[i1] / 3; - m = adja[i1] % 3; - mesh->adja[3*(k-1)+1+i1] = 3*iel+i2; - mesh->adja[3*(iel-1)+1+i2] = 3*k+i1; - mesh->adja[3*(iel-1)+1+i1] = 3*mel+m; - if(mel) - mesh->adja[3*(mel-1)+1+m] = 3*iel+i1; - - if ( jel ) { - kel = MMGS_newElt(mesh); - if ( !kel ) { - MMGS_TRIA_REALLOC(mesh,kel,mesh->gap, - MMG5_INCREASE_MEM_MESSAGE(); - if ( !MMGS_delElt(mesh,iel) ) return 0; - return 0); - } - pt = &mesh->tria[jel]; - pt1 = &mesh->tria[kel]; - pt->flag = 0; - pt->base = mesh->base; - memcpy(pt1,pt,sizeof(MMG5_Tria)); - memcpy(&mesh->adja[3*(kel-1)+1],&mesh->adja[3*(jel-1)+1],3*sizeof(int)); - - j1 = MMG5_inxt2[j]; - j2 = MMG5_iprv2[j]; - pt->v[j1] = ip; - pt->tag[j2] = MG_NOTAG; - pt->edg[j2] = 0; - pt1->v[j2] = ip; - pt1->tag[j1] = MG_NOTAG; - pt1->edg[j1] = 0; - - /* update adjacency */ - adja = &mesh->adja[3*(jel-1)+1]; - mel = adja[j2] / 3; - m = adja[j2] % 3; - mesh->adja[3*(jel-1)+1+j2] = 3*kel+j1; - mesh->adja[3*(kel-1)+1+j1] = 3*jel+j2; - mesh->adja[3*(kel-1)+1+j2] = 3*mel+m; + k = start; + i = istart; + for ( l=0; ltria[k]; + + adja = &mesh->adja[3*(k-1)+1]; + + iel = newelt[l]; + pt1 = &mesh->tria[iel]; + + /** update two triangles created by the splitting of k */ + i1 = MMG5_inxt2[i]; + i2 = MMG5_iprv2[i]; + pt->v[i2] = ip; + pt->tag[i1] = MG_NOTAG; + pt->edg[i1] = 0; + pt1->v[i1] = ip; + pt1->tag[i2] = MG_NOTAG; + pt1->edg[i2] = 0; + + /** update adjacency relations */ + mel = adja[i1] / 3; + m = adja[i1] % 3; + /* Through i1, pt is adja to pt1,i2 */ + mesh->adja[3*(k-1)+1+i1] = 3*iel+i2; + /* Through i2, pt1 is adja to pt,i1 */ + mesh->adja[3*(iel-1)+1+i2] = 3*k+i1; + /* Through i1, pt1 is adja to mel,m */ + mesh->adja[3*(iel-1)+1+i1] = 3*mel+m; + /* Through m, mel is adja to pt1,i1 */ if(mel) - mesh->adja[3*(mel-1)+1+m] = 3*kel+j2; + mesh->adja[3*(mel-1)+1+m] = 3*iel+i1; + + k = adja[i] / 3; + i = adja[i] % 3; + } + + /** Update the adjacents through edge i */ + if ( mesh->adja[3*(start-1)+1+istart] ) { + k = start; + i = istart; + + iel = newelt[0]; + + for ( l=1; ladja[3*(k-1)+1+i]; + j = adj%3; + + kel = newelt[l]; + + assert ( kel ); + + mesh->adja[3*(k-1)+1+i] = 3*kel+j; + k = adj / 3; + + mesh->adja[3*(iel-1)+1+i] = 3*k +j; + + i = j; + iel = kel; + } + + mesh->adja[3*( k-1)+1+i] = 3*newelt[0]+istart; + mesh->adja[3*(iel-1)+1+i] = 3*start +istart; + } mesh->adja[3*(iel-1)+1+i] = 3*kel+j; mesh->adja[3*(kel-1)+1+j] = 3*iel+i; From 04380254ae2888b17c39a913cb2c2ef63776e857 Mon Sep 17 00:00:00 2001 From: Algiane Date: Mon, 2 Dec 2019 13:10:31 +0100 Subject: [PATCH 06/10] Correction of bugs in split1b for nm cases. --- src/mmgs/split_s.c | 64 ++++++++++++++++++++++------------------------ 1 file changed, 30 insertions(+), 34 deletions(-) diff --git a/src/mmgs/split_s.c b/src/mmgs/split_s.c index 64f8b89e7..e4443580b 100644 --- a/src/mmgs/split_s.c +++ b/src/mmgs/split_s.c @@ -113,11 +113,11 @@ int MMGS_split1(MMG5_pMesh mesh,MMG5_pSol met,int k,int i,int *vx) { iel = MMGS_newElt(mesh); if ( !iel ) { MMGS_TRIA_REALLOC(mesh,iel,mesh->gap, - fprintf(stderr,"\n ## Error: %s: unable to allocate" - " a new element.\n",__func__); - MMG5_INCREASE_MEM_MESSAGE(); - fprintf(stderr," Exit program.\n"); - return 0); + fprintf(stderr,"\n ## Error: %s: unable to allocate" + " a new element.\n",__func__); + MMG5_INCREASE_MEM_MESSAGE(); + fprintf(stderr," Exit program.\n"); + return 0); } pt = &mesh->tria[k]; @@ -379,10 +379,6 @@ int split1b(MMG5_pMesh mesh,int start,char istart,int ip) { mesh->adja[3*(iel-1)+1+i] = 3*start +istart; } - mesh->adja[3*(iel-1)+1+i] = 3*kel+j; - mesh->adja[3*(kel-1)+1+j] = 3*iel+i; - } - return 1; } @@ -484,20 +480,20 @@ int MMGS_split2(MMG5_pMesh mesh,MMG5_pSol met,int k,int *vx) { iel = MMGS_newElt(mesh); if ( !iel ) { MMGS_TRIA_REALLOC(mesh,iel,mesh->gap, - fprintf(stderr,"\n ## Error: %s: unable to allocate" - " a new element.\n",__func__); - MMG5_INCREASE_MEM_MESSAGE(); - fprintf(stderr," Exit program.\n"); - return 0); + fprintf(stderr,"\n ## Error: %s: unable to allocate" + " a new element.\n",__func__); + MMG5_INCREASE_MEM_MESSAGE(); + fprintf(stderr," Exit program.\n"); + return 0); } jel = MMGS_newElt(mesh); if ( !jel ) { MMGS_TRIA_REALLOC(mesh,jel,mesh->gap, - fprintf(stderr,"\n ## Error: %s: unable to allocate" - " a new element.\n",__func__); - MMG5_INCREASE_MEM_MESSAGE(); - fprintf(stderr," Exit program.\n"); - return 0); + fprintf(stderr,"\n ## Error: %s: unable to allocate" + " a new element.\n",__func__); + MMG5_INCREASE_MEM_MESSAGE(); + fprintf(stderr," Exit program.\n"); + return 0); } pt = &mesh->tria[k]; @@ -644,29 +640,29 @@ int MMGS_split3(MMG5_pMesh mesh,MMG5_pSol met,int k,int *vx) { iel = MMGS_newElt(mesh); if ( !iel ) { MMGS_TRIA_REALLOC(mesh,iel,mesh->gap, - fprintf(stderr,"\n ## Error: %s: unable to allocate" - " a new element.\n",__func__); - MMG5_INCREASE_MEM_MESSAGE(); - fprintf(stderr," Exit program.\n"); - return 0); + fprintf(stderr,"\n ## Error: %s: unable to allocate" + " a new element.\n",__func__); + MMG5_INCREASE_MEM_MESSAGE(); + fprintf(stderr," Exit program.\n"); + return 0); } jel = MMGS_newElt(mesh); if ( !jel ) { MMGS_TRIA_REALLOC(mesh,jel,mesh->gap, - fprintf(stderr,"\n ## Error: %s: unable to allocate" - " a new element.\n",__func__); - MMG5_INCREASE_MEM_MESSAGE(); - fprintf(stderr," Exit program.\n"); - return 0); + fprintf(stderr,"\n ## Error: %s: unable to allocate" + " a new element.\n",__func__); + MMG5_INCREASE_MEM_MESSAGE(); + fprintf(stderr," Exit program.\n"); + return 0); } kel = MMGS_newElt(mesh); if ( !kel ) { MMGS_TRIA_REALLOC(mesh,kel,mesh->gap, - fprintf(stderr,"\n ## Error: %s: unable to allocate" - " a new element.\n",__func__); - MMG5_INCREASE_MEM_MESSAGE(); - fprintf(stderr," Exit program.\n"); - return 0); + fprintf(stderr,"\n ## Error: %s: unable to allocate" + " a new element.\n",__func__); + MMG5_INCREASE_MEM_MESSAGE(); + fprintf(stderr," Exit program.\n"); + return 0); } pt = &mesh->tria[k]; From 62a47c06f9b59e96f1ebf66d00db7203a21fb19e Mon Sep 17 00:00:00 2001 From: Algiane Date: Tue, 3 Dec 2019 09:51:35 +0100 Subject: [PATCH 07/10] Reset flag at the beginning of the chkedg function. --- src/mmgs/boulep_s.c | 2 ++ src/mmgs/mmgs1.c | 2 ++ 2 files changed, 4 insertions(+) diff --git a/src/mmgs/boulep_s.c b/src/mmgs/boulep_s.c index 8934ea30d..103e16a0a 100644 --- a/src/mmgs/boulep_s.c +++ b/src/mmgs/boulep_s.c @@ -138,6 +138,8 @@ int boulechknm(MMG5_pMesh mesh,int start,int ip,int *list) { /** Stage 2: check if a collapse may lead to a non-convex situation */ /* Flags initialization */ + base = ++mesh->base; + pt = &mesh->tria[start]; ia = MMG5_iprv2[ip]; iq = MMG5_inxt2[ip]; diff --git a/src/mmgs/mmgs1.c b/src/mmgs/mmgs1.c index b67d7d528..677d9a76f 100755 --- a/src/mmgs/mmgs1.c +++ b/src/mmgs/mmgs1.c @@ -222,6 +222,8 @@ int chkedg(MMG5_pMesh mesh,int iel) { p[1] = &mesh->point[pt->v[1]]; p[2] = &mesh->point[pt->v[2]]; + pt->flag = 0; + /* normal recovery */ for (i=0; i<3; i++) { if ( MS_SIN(p[i]->tag) ) { From f52503da48d83687728c26659f542bed48fc6902 Mon Sep 17 00:00:00 2001 From: Algiane Date: Fri, 7 Feb 2020 10:20:48 +0100 Subject: [PATCH 08/10] Corections of list of bugs in colver and chkcol. --- src/mmgs/analys_s.c | 4 +- src/mmgs/anisosiz_s.c | 18 +- src/mmgs/boulep_s.c | 201 ++++++++----- src/mmgs/chkmsh_s.c | 10 +- src/mmgs/colver_s.c | 643 ++++++++++++++++++++++++++++-------------- src/mmgs/mmgs.h | 8 +- src/mmgs/mmgs1.c | 47 +-- 7 files changed, 589 insertions(+), 342 deletions(-) diff --git a/src/mmgs/analys_s.c b/src/mmgs/analys_s.c index 6a50f8e89..d081b7dc5 100644 --- a/src/mmgs/analys_s.c +++ b/src/mmgs/analys_s.c @@ -501,7 +501,7 @@ static int MMG5_singul(MMG5_pMesh mesh) { MMG5_pTria pt; MMG5_pPoint ppt,p1,p2; double ux,uy,uz,vx,vy,vz,dd; - int list[MMGS_LMAX+2],listref[MMGS_LMAX+2],k,nc,xp,nr,ns,nre; + int list[MMGS_LMAX+2],listref[MMGS_LMAX+2],k,nc,xp,nr,ns,nre,ishell; char i; nre = nc = 0; @@ -566,7 +566,7 @@ static int MMG5_singul(MMG5_pMesh mesh) { for (i=0; i<3; i++) { ppt = &mesh->point[pt->v[i]]; if ( !ppt->s ) continue; - nr = boulet(mesh,k,i,list); + nr = MMGS_boulet(mesh,k,i,MMG5_iprv2[i],list,&ishell); if ( nr != ppt->s ) { ppt->tag |= MG_CRN + MG_REQ; ppt->s = 0; diff --git a/src/mmgs/anisosiz_s.c b/src/mmgs/anisosiz_s.c index d5da59c0d..23d9ff7f6 100644 --- a/src/mmgs/anisosiz_s.c +++ b/src/mmgs/anisosiz_s.c @@ -55,7 +55,7 @@ static int MMG5_defmetsin(MMG5_pMesh mesh,MMG5_pSol met,int it,int ip) { double *m,n[3],isqhmin,isqhmax,b0[3],b1[3],ps1,tau[3]; double ntau2,gammasec[3]; double c[3],kappa,maxkappa,alpha,hausd,hausd_v; - int ilist,list[MMGS_LMAX+2],k,i,iel,idp,init_s; + int ilist,list[MMGS_LMAX+2],k,i,iel,idp,init_s,ishell; unsigned char i0,i1,i2; pt = &mesh->tria[it]; @@ -68,8 +68,8 @@ static int MMG5_defmetsin(MMG5_pMesh mesh,MMG5_pSol met,int it,int ip) { isqhmin = mesh->info.hmin; isqhmax = mesh->info.hmax; - ilist = boulet(mesh,it,ip,list); - if ( !ilist ) + ilist = MMGS_boulet(mesh,it,ip,MMG5_iprv2[ip],list,&ishell); + if ( ilist < 1 ) return 0; maxkappa = 0.0; @@ -345,7 +345,7 @@ static int MMG5_defmetref(MMG5_pMesh mesh,MMG5_pSol met,int it,int ip) { MMG5_pPoint p0,p1; MMG5_Bezier b; MMG5_pPar par; - int i,ilist,list[MMGS_LMAX+2],k,iel,ipref[2],idp,isloc; + int i,ilist,list[MMGS_LMAX+2],k,iel,ipref[2],idp,isloc,ishell; double *m,isqhmin,isqhmax,*n,r[3][3],lispoi[3*MMGS_LMAX+1]; double ux,uy,uz,det2d,intm[3],c[3]; double tAA[6],tAb[3],hausd; @@ -357,8 +357,8 @@ static int MMG5_defmetref(MMG5_pMesh mesh,MMG5_pSol met,int it,int ip) { idp = pt->v[ip]; p0 = &mesh->point[idp]; - ilist = boulet(mesh,it,ip,list); - if ( !ilist ) + ilist = MMGS_boulet(mesh,it,ip,MMG5_iprv2[ip],list,&ishell); + if ( ilist < 1 ) return 0; /* Computation of the rotation matrix T_p0 S -> [z = 0] */ @@ -516,7 +516,7 @@ static int MMG5_defmetreg(MMG5_pMesh mesh,MMG5_pSol met,int it,int ip) { MMG5_pPoint p0,p1; MMG5_Bezier b; MMG5_pPar par; - int ilist,list[MMGS_LMAX+2],k,iel,idp,isloc,i; + int ilist,list[MMGS_LMAX+2],k,iel,idp,isloc,i,ishell; double *n,*m,r[3][3],ux,uy,uz,lispoi[3*MMGS_LMAX+1]; double det2d,c[3],isqhmin,isqhmax; double tAA[6],tAb[3],hausd; @@ -526,8 +526,8 @@ static int MMG5_defmetreg(MMG5_pMesh mesh,MMG5_pSol met,int it,int ip) { idp = pt->v[ip]; p0 = &mesh->point[idp]; - ilist = boulet(mesh,it,ip,list); - if ( !ilist ) + ilist = MMGS_boulet(mesh,it,ip,MMG5_iprv2[ip],list,&ishell); + if ( ilist < 1 ) return 0; /* Computation of the rotation matrix T_p0 S -> [z = 0] */ diff --git a/src/mmgs/boulep_s.c b/src/mmgs/boulep_s.c index 103e16a0a..dbd368c29 100644 --- a/src/mmgs/boulep_s.c +++ b/src/mmgs/boulep_s.c @@ -39,19 +39,32 @@ * \param mesh pointer toward the mesh structure. * \param start index of triangle to start. * \param ip index of point for wich we compute the ball. + * \param ia index of edge in tria \a start from which we start to compute + * the ball (so edge \a ia has \a ip as extremity). * \param list pointer toward the computed ball of \a ip. + * \param ishell number of tria in the shell of \a ia (to fill) + * * \return the size of the computed ball (ilist) or -ilist if fail. * - * Find all triangles sharing \a ip, \f$list[0] =\f$ \a start do not stop when - * crossing ridge. + * Find all triangles sharing \a ip, starting from \a start and storing first + * the \a ishell trias in the shell of \a ia. Do not stop when crossing + * ridge. Works for a non-manifold edge and if we cross nm edges within the ball of \a ip. * */ -int boulet(MMG5_pMesh mesh,int start,int ip,int *list) { +int MMGS_boulet(MMG5_pMesh mesh,int start,int ip,int ia,int *list, int *ishell) { MMG5_pTria pt,pt1; MMG5_pPoint ppt; int *adja,cur,k,ilist,nump,k1,j,base; char i,l; + /* Check that ip is at extremity of edge ia */ + assert ( ip != ia ); + + /* Functions related to collapse operators (boulechknm, colver...) implicitely + * state that iq == MMG5_inxt2[ip] which implies ia == MMG5_iprv2[ip]. */ + assert ( MMG5_iprv2[ip] == ia ); + + /* Initialization */ base = ++mesh->base; pt = &mesh->tria[start]; @@ -61,17 +74,47 @@ int boulet(MMG5_pMesh mesh,int start,int ip,int *list) { /* store initial tria */ pt->flag = base; list[0] = 3*start + ip; - ilist = 1; + *ishell = 1; - /* Explore list and travel by adjacency through elements sharing nump */ - cur = 0; + /** First stage: Fill the shell of ia */ + cur = mesh->adja[ 3*(start-1)+ 1 +ia ]; + while ( cur && (cur != 3*start+ia) ) { + k = cur/3; + i = cur%3; + /* Store tria k and look for nump index */ + mesh->tria[k].flag = base; + + j = MMG5_inxt2[i]; + if ( mesh->tria[k].v[j] != nump ) { + j = MMG5_inxt2[j]; + assert ( mesh->tria[k].v[j] == nump ); + } + + list[(*ishell)++] = 3*k + j; + cur = mesh->adja[ 3*(k-1)+ 1 +i ]; + } + +#ifndef NDEBUG + if ( pt->tag[ia] & MG_NOM ) { + assert ( (*ishell) >=3 ); + } + else { + assert ( (*ishell) < 3 ); + } +#endif + + /** Second stage: From each tria of the shell, explore list and travel by + * adjacency through elements sharing nump */ + cur = 0; + ilist = *ishell; while ( cur < ilist ) { k = list[cur] / 3; i = list[cur] % 3; // index of point nump in tria k adja = &mesh->adja[ 3*(k-1)+ 1 ]; + /* Pile up the 2 adjacent of k */ for (l=0; l<2; l++) { i = MMG5_inxt2[i]; k1 = adja[i]; @@ -92,7 +135,7 @@ int boulet(MMG5_pMesh mesh,int start,int ip,int *list) { assert ( pt1->v[j] == nump ); /* Check for overflow */ - if ( ilist > MMGS_LMAX-2 ) return 0; + if ( ilist > MMGS_LMAX-2 ) return -ilist; /* Store adjacent */ list[ilist] = 3*k1 + j; @@ -106,103 +149,115 @@ int boulet(MMG5_pMesh mesh,int start,int ip,int *list) { /** * \param mesh pointer toward the mesh structure. - * \param start index of tetra to start to compute the ball. - * \param ip index of point in tetra \a start for which we want to compute + * \param start index of tria to start to compute the ball. + * \param ip index of point in tria \a start for which we want to compute * the ball. + * \param ie index of edge in tria \a start from which we start to compute + * the ball (so edge \a ie has \a ip as extremity). * \param list pointer toward the computed ball of point. + * \param ishell number of tria in the shell of \a ie (to fill) * * \return -ilist (number of tria in ball) if buffer overflow, ilist if the * collapse is ok, 0 if the collapse lead to a non manifold situation. * - * Find all triangles sharing \a ip, \f$list[0] = start\f$. Do not stop when - * crossing ridge. Check whether resulting configuration is manifold. + * Find all triangles sharing \a ip, starting from \a start and storing first + * the \a ishell tria in the shell of \a ie. Do not stop when crossing + * ridge. Check whether resulting configuration is manifold. * */ -int boulechknm(MMG5_pMesh mesh,int start,int ip,int *list) { +int MMGS_boulechknm(MMG5_pMesh mesh,int start,int ip,int ie,int *list,int *ishell) { MMG5_pTria pt; MMG5_pPoint ppt; - int *adja,k,ilist,base,iel; - char i,i1,i2,ia,iq,voy; + int *adja,k,k1,ilist,base,iel,listq[MMGS_LMAX+2],ilistq,numq,cur; + char i,j,iq,voy; pt = &mesh->tria[start]; if ( !MG_EOK(pt) ) return 0; - -#warning ajeter : on, fait quoi dans ce cas? +#warning ajeter : old test : Ca doit marcher non? if ( mesh->point[pt->v[ip]].tag & MG_NOM ) return 0; - /** Stage 1: fill the ball of point */ - ilist = boulet(mesh,start,ip,list); + /** Stage 1: fill the ball of point and store it in \a list (such as + * list[k]%3 == ip ) */ + ilist = MMGS_boulet(mesh,start,ip,ie,list,ishell); if ( ilist < 1 ) return ilist; /** Stage 2: check if a collapse may lead to a non-convex situation */ - /* Flags initialization */ + /* 1. Flags initialization: set all the flags of the points of the ball to base */ base = ++mesh->base; + for ( k=(*ishell); ktria[iel]; + mesh->point[pt->v[MMG5_inxt2[i]]].s = base; + mesh->point[pt->v[MMG5_iprv2[i]]].s = base; + } - pt = &mesh->tria[start]; - ia = MMG5_iprv2[ip]; - iq = MMG5_inxt2[ip]; - - i1 = ia; - i2 = iq; - ppt = &mesh->point[pt->v[i1]]; - ppt->s = 0; - ppt = &mesh->point[pt->v[i2]]; - ppt->s = 0; - - adja = &mesh->adja[3*(start-1)+1]; - iel = adja[ia] / 3; - voy = adja[ia] % 3; - - if( iel ) { - pt = &mesh->tria[iel]; - ppt = &mesh->point[pt->v[voy]]; - ppt->s = 0; + /* 2. Flags initialization: reset the fkags of the points of the triangles of + * the shell ip-iq to 0. */ + for ( k=0; k<(*ishell); ++k ) { + iel = list[k]/3; + i = list[k]%3; + pt = &mesh->tria[iel]; + mesh->point[pt->v[MMG5_inxt2[i]]].s = 0; + mesh->point[pt->v[MMG5_iprv2[i]]].s = 0; } - /* Unfold the ball af iq and check that it doesn't contains common tria with + /* 3. Unfold the ball af iq and check that it doesn't contains common tria with * the ball of ip (except for the ip-iq shell) */ - k = start; - i = iq; - do { - pt = &mesh->tria[k]; + pt = &mesh->tria[start]; + j = MMG5_iprv2[ip]; + assert ( j == ie && "Input edge index for the collapse don't match with the computed one." ); + iq = MMG5_inxt2[ip]; + numq = pt->v[iq]; - i1 = MMG5_inxt2[i]; - i2 = MMG5_iprv2[i]; - ppt = &mesh->point[pt->v[i1]]; - if ( ppt->s == base ) return 0; - ppt = &mesh->point[pt->v[i2]]; - if ( ppt->s == base ) return 0; + listq[0] = 3*start + iq; + pt->flag = base; - adja = &mesh->adja[3*(k-1)+1]; - k = adja[i1] / 3; - i = adja[i1] % 3; - i = MMG5_inxt2[i]; - } - while ( k && k != start ); - if( k > 0 ) return ilist; + cur = 0; + ilistq = 1; + while ( cur < ilistq ) { + k = listq[cur] / 3; + i = listq[cur] % 3; // index of point iq in tria k - /* Unfold the ball in opposite direction if boundary hit */ - k = start; - i = iq; - do { - adja = &mesh->adja[3*(k-1)+1]; - i1 = MMG5_inxt2[i]; - i2 = MMG5_iprv2[i]; + adja = &mesh->adja[ 3*(k-1)+ 1 ]; + + /* Pile up the 2 adjacent of k */ + for (voy=0; voy<2; voy++) { + i = MMG5_inxt2[i]; + k1 = adja[i]; + if ( !k1 ) continue; + + j = k1%3; + k1 /= 3; + + pt = &mesh->tria[k1]; + if ( pt->flag == base ) continue; + pt->flag = base; + + /* Find the local index of numq */ + j = MMG5_inxt2[j]; + if ( pt->v[MMG5_inxt2[j]] == numq ) { + j = MMG5_inxt2[j]; + } + assert ( pt->v[j] == numq ); - pt = &mesh->tria[k]; - ppt = &mesh->point[pt->v[i1]]; - if ( ppt->s == base ) return 0; - ppt = &mesh->point[pt->v[i2]]; - if ( ppt->s == base ) return 0; + ppt = &mesh->point[pt->v[MMG5_inxt2[j]]]; + if ( ppt->s == base ) return 0; - k = adja[i2] / 3; - if ( k == 0 ) break; - i = adja[i2] % 3; - i = MMG5_iprv2[i]; + ppt = &mesh->point[pt->v[MMG5_iprv2[j]]]; + if ( ppt->s == base ) return 0; + + /* Check for overflow */ + if ( ilistq > MMGS_LMAX-2 ) return -ilist; + + /* Store adjacent */ + listq[ilistq] = 3*k1 + j; + ++ilistq; + } + ++cur; } - while ( k ); return ilist; } diff --git a/src/mmgs/chkmsh_s.c b/src/mmgs/chkmsh_s.c index 9356dc7d1..f7b48d4fc 100644 --- a/src/mmgs/chkmsh_s.c +++ b/src/mmgs/chkmsh_s.c @@ -142,7 +142,7 @@ int MMG5_mmgsChkmsh(MMG5_pMesh mesh,int severe,int base) { MMG5_pPoint ppt; MMG5_pTria pt1,pt2; int adj,adj1,k,kk,l,nk,i,j,ip,lon,len; - int *adja,*adjb,list[MMGS_LMAX+2]; + int *adja,*adjb,list[MMGS_LMAX+2],ishell; char voy,voy1,i1,i2,j1,j2; static char mmgErr0=0,mmgErr1=0,mmgErr2=0,mmgErr3=0,mmgErr4=0; static char mmgErr5=0,mmgErr6=0,mmgErr7=0; @@ -282,8 +282,12 @@ int MMG5_mmgsChkmsh(MMG5_pMesh mesh,int severe,int base) { } else if ( MS_SIN(ppt->tag) ) continue; - lon = boulet(mesh,k,i,list); - if ( lon < 1 ) continue; + lon = MMGS_boulet(mesh,k,i,MMG5_inxt2[i],list,&ishell); + if ( lon < 1 ) { + fprintf(stderr,"\n ## Warning: %s: 6. overflow in the ball of %d (from %d)\n", + __func__,MMGS_indPt(mesh,ip),MMGS_indElt(mesh,k)); + continue; + } for (l=0; ltria[k]; + *fwd = 1; + if ( pt->v[MMG5_inxt2[i]] != ip1 ) { + *fwd = 0; + assert(pt->v[MMG5_iprv2[i]] == ip1); + } + do { + ++ilist; + pt1 = &mesh->tria[jel]; + adja = &mesh->adja[3*(jel-1)+1]; + if ( *fwd ) { + j1 = MMG5_inxt2[j]; + j2 = MMG5_iprv2[j]; + } + else { + j1 = MMG5_iprv2[j]; + j2 = MMG5_inxt2[j]; + } + assert ( pt1->v[j1] == ip1 ); + + if ( pt1->tag[j2] & MG_NOM ) break; + + jel = adja[j2]/3; + j = adja[j2]%3; + } while ( jel && (jel != k) ); + + return ilist; +} + +#define MMGS_UPDATE_PIV(ia,ia1,ia2,fwd) do \ + { \ + if ( fwd==1 ) { \ + ia1 = MMG5_inxt2[ia]; \ + ia2 = MMG5_iprv2[ia]; \ + } \ + else { \ + ia1 = MMG5_iprv2[ia]; \ + ia2 = MMG5_inxt2[ia]; \ + } \ + } while(0) /** * \param mesh pointer toward the mesh * \param met pointer toward the metric - * \param k index of the element in wich we collapse - * \param i index of the edge to collapse - * \param list pointer toward the ball of point + * \param start index of the element in wich we collapse + * \param istart index of the edge to collapse + * \param list pointer toward the ball of point \a ip1 + * \param ishell number of tria in the shell of \a istart * \param typchk type of check to perform * * \return 0 if we can't move of if we fail, 1 if success * - * check if geometry preserved by collapsing edge i + * check if geometry preserved by collapsing non-manifold edge \a istart + * (collapse \a ip1 over \a ip2). * */ -int chkcol(MMG5_pMesh mesh,MMG5_pSol met,int k,char i,int *list,char typchk) { +int MMGS_chkcol(MMG5_pMesh mesh,MMG5_pSol met,int start,char istart,int *list, + int *ishell,char typchk) { MMG5_pTria pt,pt0,pt1,pt2; MMG5_pPoint p1,p2; double len,lon,ps,cosnold,cosnnew,kal,n0old[3],n1old[3],n00old[3]; double n0new[3],n1new[3],n00new[3]; - int *adja,jel,kel,ip1,ip2,l,ll,ilist; - char i1,i2,j,jj,j2,lj,open,voy; + int *adja,jel,kel,ip1,ip2,nump,l,ll,base,kk,k,i; + int ilistfull,ilist,sumilist,fwd; + char i1,i2,j,jj,j1,j2,lj,lj1,lj2,open,voy; pt0 = &mesh->tria[0]; - pt = &mesh->tria[k]; - i1 = MMG5_inxt2[i]; - i2 = MMG5_iprv2[i]; + pt = &mesh->tria[start]; + i1 = MMG5_inxt2[istart]; + i2 = MMG5_iprv2[istart]; ip1 = pt->v[i1]; ip2 = pt->v[i2]; + nump = ip1; + + assert ( !MG_SIN(mesh->point[ip1].tag) ); + if ( typchk == 2 && met->m ) { lon = MMG5_lenSurfEdg(mesh,met,ip1,ip2,0); if ( !lon ) return 0; @@ -72,79 +143,223 @@ int chkcol(MMG5_pMesh mesh,MMG5_pSol met,int k,char i,int *list,char typchk) { } /* collect all triangles around vertex i1 */ - ilist = boulechknm(mesh,k,i1,list); - if ( ilist <= 0 ) return 0; + ilistfull = MMGS_boulechknm(mesh,start,i1,istart,list,ishell); + if ( ilistfull <= 0 ) return 0; /* check for open ball */ - adja = &mesh->adja[3*(k-1)+1]; - open = adja[i] == 0; + adja = &mesh->adja[3*(start-1)+1]; + // open = adja[istart] == 0; + + /** Starting from each tria of the shell of i, travel by adjacency around \a + * ip1 and process the portion of surface until meeting a non-manifold edge */ + base = ++mesh->base; + sumilist = 0; + for ( kk = 0; kk<*ishell; ++ kk ) { + k = list[kk]/3; + i = list[kk]%3; + + adja = &mesh->adja[3*(k-1)+1]; + + pt = &mesh->tria[k]; + if ( pt->base == base ) continue; + pt->base = base; + + /* First tria of the portion of surface */ + jel = k; + j = i; + i1 = i; + i2 = MMG5_inxt2[i1]; + i = MMG5_iprv2[i1]; + if ( pt->v[i2]!= ip2 ) { + i = i2; + i2 = MMG5_inxt2[i2]; + } + assert ( (pt->v[i1] == nump) && (pt->v[i2] == ip2) ); + + fwd = 1; + if ( pt->tag[i] & MG_NOM ) { + //assert ( !open ); + assert ( *ishell > 2 ); + /** Count the number of tria that we can reach without crossing a nm edge */ + ilist = MMGS_count_bouleploc( mesh, k, i, ip1,&fwd); + open = 1; + } + else { + assert ( *ishell < 3 ); + ilist = ilistfull; + open = ( adja[i]==0 ); + if ( !open ) { + /* Closed manifold case: mark the last tria of the ball as treated to + * not unfold the ball a second time */ + mesh->tria[adja[i]/3].base = base; + } + } + sumilist += ilist; - if ( ilist > 3 ) { - /* check references */ + /* check references: */ if ( MG_EDG(pt->tag[i2]) ) { jel = list[1] / 3; pt1 = &mesh->tria[jel]; if ( abs(pt->ref) != abs(pt1->ref) ) return 0; } - /* analyze ball */ - assert ( ilist-1+open > 1 ); - for (l=1; l 3 ) { + /* If only one pt is the only elt of a given ref, refuse the collapse */ + jel = adja[i2] / 3; + assert ( jel ); pt1 = &mesh->tria[jel]; - /* check length */ - if ( typchk == 2 && met->m && !MG_EDG(mesh->point[ip2].tag) ) { - ip1 = pt1->v[j2]; - len = MMG5_lenSurfEdg(mesh,met,ip1,ip2,0); - if ( len > lon || !len ) return 0; + j = adja[i2] % 3; + MMGS_UPDATE_PIV(j,j1,j2,fwd); + assert ( pt1->v[j1] == nump ); + + if ( MG_EDG(pt->tag[i2]) ) { + if ( abs(pt->ref) != abs(pt1->ref) ) return 0; } - /* check normal flipping */ - if ( !MMG5_nortri(mesh,pt1,n1old) ) return 0; - memcpy(pt0,pt1,sizeof(MMG5_Tria)); - pt0->v[j] = ip2; - if ( !MMG5_nortri(mesh,pt0,n1new) ) return 0; + /* analyze ball */ + assert ( ilist-1+open > 1 ); - ps = n1new[0]*n1old[0] + n1new[1]*n1old[1] + n1new[2]*n1old[2]; - if ( ps < 0.0 ) return 0; + for ( l=1; ltria[jel]; + MMGS_UPDATE_PIV(j,j1,j2,fwd); + assert ( pt1->v[j1] == nump ); - /* keep normals at 1st triangles */ - if ( l == 1 && !open ) { - memcpy(n00old,n1old,3*sizeof(double)); - memcpy(n00new,n1new,3*sizeof(double)); - } + pt1->flag = base; - /* check normals deviation */ - if ( !(pt1->tag[j2] & MG_GEO) ) { - if ( l > 1 ) { - cosnold = n0old[0]*n1old[0] + n0old[1]*n1old[1] + n0old[2]*n1old[2]; - cosnnew = n0new[0]*n1new[0] + n0new[1]*n1new[1] + n0new[2]*n1new[2]; - if ( cosnold < MMG5_ANGEDG ) { - if ( cosnnew < cosnold ) return 0; + /* check length */ + if ( typchk == 2 && met->m && !MG_EDG(mesh->point[ip2].tag) ) { + ip1 = pt1->v[j2]; + len = MMG5_lenSurfEdg(mesh,met,ip1,ip2,0); + if ( len > lon || !len ) return 0; + } + + /* check normal flipping */ + if ( !MMG5_nortri(mesh,pt1,n1old) ) return 0; + memcpy(pt0,pt1,sizeof(MMG5_Tria)); + pt0->v[j1] = ip2; + if ( !MMG5_nortri(mesh,pt0,n1new) ) return 0; + + ps = n1new[0]*n1old[0] + n1new[1]*n1old[1] + n1new[2]*n1old[2]; + if ( ps < 0.0 ) return 0; + + /* keep normals at 1st triangles */ + if ( l == 1 && !open ) { + memcpy(n00old,n1old,3*sizeof(double)); + memcpy(n00new,n1new,3*sizeof(double)); + } + + /* check normals deviation */ + if ( !(pt1->tag[j] & MG_GEO) ) { + if ( l > 1 ) { + cosnold = n0old[0]*n1old[0] + n0old[1]*n1old[1] + n0old[2]*n1old[2]; + cosnnew = n0new[0]*n1new[0] + n0new[1]*n1new[1] + n0new[2]*n1new[2]; + if ( cosnold < MMG5_ANGEDG ) { + if ( cosnnew < cosnold ) return 0; + } + else if ( cosnnew < MMG5_ANGEDG ) return 0; } - else if ( cosnnew < MMG5_ANGEDG ) return 0; } - } - /* check geometric support */ - if ( l == 1 ) { - pt0->tag[j2] |= pt->tag[i1]; + /* check geometric support */ + if ( l == 1 ) { + pt0->tag[j] |= pt->tag[i1]; + } + else if ( l == ilist-2+open ) { + if ( !open ) { + ll = mesh->adja[3*(jel-1)+1+j2] / 3; + lj = mesh->adja[3*(jel-1)+1+j2] % 3; + MMGS_UPDATE_PIV(lj,lj1,lj2,fwd); + assert ( mesh->tria[ll].v[lj1]==nump ); + pt0->tag[j2] |= mesh->tria[ll].tag[lj1]; + } + else { + pt0->tag[j2] |= pt->tag[i]; + } + } + if ( chkedg(mesh,0) ) return 0; + + /* check quality */ + if ( typchk == 2 && met->m ) + kal = MMGS_ALPHAD*MMG5_calelt(mesh,met,pt0); + else + kal = MMGS_ALPHAD*MMG5_caltri_iso(mesh,NULL,pt0); + if ( kal < MMGS_NULKAL ) return 0; + + memcpy(n0old,n1old,3*sizeof(double)); + memcpy(n0new,n1new,3*sizeof(double)); + + if ( pt1->tag[j2] & MG_NOM ) break; + + adja = &mesh->adja[3*(jel-1)+1]; + + jel = adja[j2]/3; + j = adja[j2]%3; } - else if ( l == ilist-2+open ) { - if ( !open ) { - ll = list[ilist-1] / 3; - lj = list[ilist-1] % 3; - pt0->tag[jj] |= mesh->tria[ll].tag[lj]; + + /* check angle between 1st and last triangles */ + if ( (!open) && !(pt->tag[i] & MG_GEO) ) { + cosnold = n00old[0]*n1old[0] + n00old[1]*n1old[1] + n00old[2]*n1old[2]; + cosnnew = n00new[0]*n1new[0] + n00new[1]*n1new[1] + n00new[2]*n1new[2]; + if ( cosnold < MMG5_ANGEDG ) { + if ( cosnnew < cosnold ) return 0; } - else { - assert ( list[0]/3 == k ); - pt0->tag[jj] |= pt->tag[i]; + else if ( cosnnew < MMG5_ANGEDG ) return 0; + + /* other checks for reference collapse */ + adja = &mesh->adja[3*(jel-1)+1]; + jel = adja[j]/3; + j = adja[j]%3; + MMGS_UPDATE_PIV(j,j1,j2,-fwd); + + pt = &mesh->tria[jel]; + if ( MG_EDG(pt->tag[j2]) ) { + adja = &mesh->adja[3*(jel-1)+1]; + jel = adja[j2] / 3; + pt1 = &mesh->tria[jel]; + if ( abs(pt->ref) != abs(pt1->ref) ) return 0; } } + } + + /* specific test: no collapse if any interior edge is EDG */ + else if ( ilist == 3 ) { + /* Remark: if ilist==3 and i is an open ridge, we pass in the previous + * test (open+ilist > 3) so here, ip1 is in the middle of the 3 + * triangles */ + p1 = &mesh->point[pt->v[i1]]; + if ( MS_SIN(p1->tag) ) return 0; + else if ( MG_EDG(pt->tag[i2]) && !MG_EDG(pt->tag[i]) ) return 0; + else if ( !MG_EDG(pt->tag[i2]) && MG_EDG(pt->tag[i]) ) return 0; + else if ( MG_EDG(pt->tag[i2]) && MG_EDG(pt->tag[i]) && MG_EDG(pt->tag[i1]) ) return 0; + + /* Check geometric approximation */ + jel = adja[i2] / 3; + assert ( jel ); + pt1 = &mesh->tria[jel]; + pt1->base = base; + + j = adja[i2] % 3; + MMGS_UPDATE_PIV(j,j1,j2,fwd); + assert ( pt1->v[j1] == nump ); + + pt0 = &mesh->tria[0]; + memcpy(pt0,pt1,sizeof(MMG5_Tria)); + pt0->v[j1] = ip2; + + adja = &mesh->adja[3*(jel-1)+1]; + + jel = adja[j2] / 3; + assert ( jel ); + pt1 = &mesh->tria[jel]; + pt1->base = base; + + lj = adja[j2] % 3; + MMGS_UPDATE_PIV(lj,lj1,lj2,fwd); + + pt0->tag[j2] |= pt1->tag[lj1]; + pt0->tag[j ] |= pt1->tag[lj2]; if ( chkedg(mesh,0) ) return 0; /* check quality */ @@ -153,131 +368,169 @@ int chkcol(MMG5_pMesh mesh,MMG5_pSol met,int k,char i,int *list,char typchk) { else kal = MMGS_ALPHAD*MMG5_caltri_iso(mesh,NULL,pt0); if ( kal < MMGS_NULKAL ) return 0; - - memcpy(n0old,n1old,3*sizeof(double)); - memcpy(n0new,n1new,3*sizeof(double)); } - /* check angle between 1st and last triangles */ - if ( !open && !(pt->tag[i] & MG_GEO) ) { - cosnold = n00old[0]*n1old[0] + n00old[1]*n1old[1] + n00old[2]*n1old[2]; - cosnnew = n00new[0]*n1new[0] + n00new[1]*n1new[1] + n00new[2]*n1new[2]; - if ( cosnold < MMG5_ANGEDG ) { - if ( cosnnew < cosnold ) return 0; - } - else if ( cosnnew < MMG5_ANGEDG ) return 0; + /* for specific configurations along open ridge */ + else if ( ilist == 2 ) { + if ( !open ) return 0; - /* other checks for reference collapse */ - jel = list[ilist-1] / 3; - j = list[ilist-1] % 3; - j = MMG5_iprv2[j]; - pt = &mesh->tria[jel]; - if ( MG_EDG(pt->tag[j]) ) { - jel = list[ilist-2] / 3; - pt1 = &mesh->tria[jel]; - if ( abs(pt->ref) != abs(pt1->ref) ) return 0; - } - } - } + jel = adja[i2] / 3; + assert ( jel ); + pt1 = &mesh->tria[jel]; + pt1->base = base; - /* specific test: no collapse if any interior edge is EDG */ - else if ( ilist == 3 ) { + j = adja[i2] % 3; + MMGS_UPDATE_PIV(j,j1,j2,fwd); + assert ( pt1->v[j1] == nump ); - p1 = &mesh->point[pt->v[i1]]; - if ( MS_SIN(p1->tag) ) return 0; - else if ( MG_EDG(pt->tag[i2]) && !MG_EDG(pt->tag[i]) ) return 0; - else if ( !MG_EDG(pt->tag[i2]) && MG_EDG(pt->tag[i]) ) return 0; - else if ( MG_EDG(pt->tag[i2]) && MG_EDG(pt->tag[i]) && MG_EDG(pt->tag[i1]) ) return 0; + /* Topological test */ + adja = &mesh->adja[3*(jel-1)+1]; + kel = adja[j1] / 3; + voy = adja[j1] % 3; + pt2 = &mesh->tria[kel]; - /* Check geometric approximation */ - jel = list[1] / 3; - j = list[1] % 3; - jj = MMG5_inxt2[j]; - j2 = MMG5_iprv2[j]; - pt0 = &mesh->tria[0]; - pt1 = &mesh->tria[jel]; - memcpy(pt0,pt1,sizeof(MMG5_Tria)); - pt0->v[j] = ip2; + if ( pt2->v[voy] == ip2) return 0; - jel = list[2] / 3; - j = list[2] % 3; - pt1 = &mesh->tria[jel]; - pt0->tag[jj] |= pt1->tag[j]; - pt0->tag[j2] |= pt1->tag[MMG5_inxt2[j]]; - if ( chkedg(mesh,0) ) return 0; - - /* check quality */ - if ( typchk == 2 && met->m ) - kal = MMGS_ALPHAD*MMG5_calelt(mesh,met,pt0); - else - kal = MMGS_ALPHAD*MMG5_caltri_iso(mesh,NULL,pt0); - if ( kal < MMGS_NULKAL ) return 0; - } + if ( fwd ) { + jj = MMG5_inxt2[j1]; + j2 = MMG5_iprv2[j1]; + } + else { + jj = MMG5_iprv2[j1]; + j2 = MMG5_inxt2[j1]; + } + if ( abs(pt->ref) != abs(pt1->ref) ) return 0; + else if ( !(pt1->tag[jj] & MG_GEO) ) return 0; - /* for specific configurations along open ridge */ - else if ( ilist == 2 ) { - if ( !open ) return 0; + p1 = &mesh->point[pt->v[i1]]; + p2 = &mesh->point[pt1->v[j]]; + if ( p2->tag > p1->tag || p2->ref != p1->ref ) return 0; - jel = list[1] / 3; - j = list[1] % 3; + /* Check geometric approximation */ + pt0 = &mesh->tria[0]; + memcpy(pt0,pt,sizeof(MMG5_Tria)); + pt0->v[i1] = pt1->v[j2]; - /* Topological test */ - adja = &mesh->adja[3*(jel-1)+1]; - kel = adja[j] / 3; - voy = adja[j] % 3; - pt2 = &mesh->tria[kel]; + if ( chkedg(mesh,0) ) return 0; + + /* check quality */ + if ( typchk == 2 && met->m ) + kal = MMGS_ALPHAD*MMG5_calelt(mesh,met,pt0); + else + kal = MMGS_ALPHAD*MMG5_caltri_iso(mesh,NULL,pt0); + if ( kal < MMGS_NULKAL ) return 0; + } + } - if ( pt2->v[voy] == ip2) return 0; + assert ( sumilist == ilistfull && + "ip1 is at intersection of multiple nm edges but non singular... " ); - jj = MMG5_inxt2[j]; - pt1 = &mesh->tria[jel]; - if ( abs(pt->ref) != abs(pt1->ref) ) return 0; - else if ( !(pt1->tag[jj] & MG_GEO) ) return 0; + return ilistfull; +} - p1 = &mesh->point[pt->v[i1]]; - p2 = &mesh->point[pt1->v[jj]]; - if ( p2->tag > p1->tag || p2->ref != p1->ref ) return 0; +/** + * \param mesh pointer toward the mesh structure + * \param list ball of the point to collapse, \f$list[k]\%3\f$ being the local + * index of the point in tetra \f$list[k]/3\f$. + * \param ishell number of triangles in the shell of the collapsed edge. + * \param ilist number of triangles in the ball of the point to collapse. + * + * \return 0 if fail, 1 if success. + * + * Collapse point \f$list[k]\%3\f$ (i1) over point \f$ MMG5_inxt2[list[k]\%3]\f$ + * (i2) in tria \f$list[k]/3\f$ (iel). + * + */ +int colver(MMG5_pMesh mesh,int *list,int ishell,int ilist) { + MMG5_pTria ptstart,pt,pt1; + int *adja,start,k,iel,jel,kel,ip1,ip2,iadr,fwd; + char istart,istart1,istart2,i,i1,i2,j,jj,j1,j2,kk; + + /** Keep the index in the starting tetra in memory */ + start = list[0] / 3; + istart1 = list[0] % 3; + + /* Index of the edge along which we collapse */ + istart = MMG5_iprv2[istart1]; + istart2 = MMG5_inxt2[istart1]; + ptstart = &mesh->tria[start]; + ip1 = ptstart->v[istart1]; + ip2 = ptstart->v[istart2]; + + /* We are not able to collapse a singular point (which implies that ip1 is not + * a non manifold point or is at the extremity of exactly 2 non-manifold + * edges) */ + assert ( !MG_SIN(mesh->point[ip1].tag) ); + + /* update adjacents */ + for (k=0; ktria[0]; - memcpy(pt0,pt,sizeof(MMG5_Tria)); - pt0->v[i1] = pt1->v[j2]; + pt = &mesh->tria[iel]; + assert ( ip1 == pt->v[i1]); + + /* Look for ip2 to reach the adjacent that need to be updated */ + i2 = MMG5_inxt2[i1]; + i = MMG5_iprv2[i1]; + fwd = 1; + if ( ip2 != pt->v[i2] ) { + i = i2; + i2 = MMG5_inxt2[i2]; + fwd = 0; + } + assert ( ip2 == pt->v[i2] ); - if ( chkedg(mesh,0) ) return 0; + /* Check that we dont collapse a domain with only one element (must be + * refused by chkcol) */ + #warning To remove: is it possible that the edge i2 is nm (collapse of a point over a nm point along a non nm edge) + assert ( !(pt->tag[i2] & MG_NOM) ); - /* check quality */ - if ( typchk == 2 && met->m ) - kal = MMGS_ALPHAD*MMG5_calelt(mesh,met,pt0); - else - kal = MMGS_ALPHAD*MMG5_caltri_iso(mesh,NULL,pt0); - if ( kal < MMGS_NULKAL ) return 0; + adja = &mesh->adja[3*(iel-1)+1]; - } + if ( !adja[i2] ) { + continue; + } - return ilist; -} + jel = adja[i2]/3; + j = adja[i2]%3; -/* collapse edge i of k, i1->i2 */ -int colver(MMG5_pMesh mesh,int *list,int ilist) { - MMG5_pTria pt,pt1,pt2; - int *adja,k,iel,jel,kel,ip1,ip2; - char i,i1,i2,j,jj,open; + pt1 = &mesh->tria[jel]; + if ( ishell!=2 || ilist!=3 || k!=1 ) { + /* If ishell == 2 && ilist == 3 && k==1 pt1->v[j1] has been updated at + * previous iteration */ + if ( fwd ) { + j1 = MMG5_inxt2[j]; + j2 = MMG5_iprv2[j]; + assert ( ip1 == pt1->v[j1] ); + assert ( pt1->v[j2] == pt->v[i] ); + } + else { + j2 = MMG5_inxt2[j]; + j1 = MMG5_iprv2[j]; + assert ( ip1 == pt1->v[j1] ); + assert ( pt1->v[j2] == pt->v[i] ); + } + /* Update pt1 vertex (ip1 becomes ip2) */ + pt1->v[j1] = ip2; + } - iel = list[0] / 3; - i1 = list[0] % 3; - i = MMG5_iprv2[i1]; - i2 = MMG5_inxt2[i1]; - pt = &mesh->tria[iel]; - ip1 = pt->v[i1]; - ip2 = pt->v[i2]; + /* Update tag and ref of edge j of pt1 (merged with edge i1 of pt) */ + pt1->tag[j] |= pt->tag[i1]; + pt1->edg[j] = MG_MAX( pt1->edg[j], pt->edg[i1] ); - /* check for open ball */ - adja = &mesh->adja[3*(iel-1)+1]; - open = adja[i] == 0; + /* Update adjacency relationships */ + iadr = adja[i1]; + kel = iadr/3; + kk = iadr%3; + + mesh->adja[3*(kel-1)+kk+1] = 3*jel + j; + mesh->adja[3*(jel-1)+j +1] = iadr; + } - /* update vertex ip1 -> ip2 */ - for (k=1; k ip2 in the ball of \a ip1 */ + for ( k=ishell; ktria[jel]; @@ -285,55 +538,12 @@ int colver(MMG5_pMesh mesh,int *list,int ilist) { pt1->base = mesh->base; } - /* update adjacent with 1st elt */ - jel = list[1] / 3; - jj = list[1] % 3; - j = MMG5_iprv2[jj]; - pt1 = &mesh->tria[jel]; - pt1->tag[j] |= pt->tag[i1]; - pt1->edg[j] = MG_MAX(pt->edg[i1],pt1->edg[j]); - if ( adja[i1] ) { - kel = adja[i1] / 3; - k = adja[i1] % 3; - mesh->adja[3*(kel-1)+1+k] = 3*jel + j; - mesh->adja[3*(jel-1)+1+j] = 3*kel + k; - pt2 = &mesh->tria[kel]; - pt2->tag[k] |= pt1->tag[j]; - pt2->edg[k] = MG_MAX(pt1->edg[j],pt2->edg[k]); - } - else - mesh->adja[3*(jel-1)+1+j] = 0; - - /* adjacent with last elt */ - if ( !open ) { - iel = list[ilist-1] / 3; - i1 = list[ilist-1] % 3; - pt = &mesh->tria[iel]; - - jel = list[ilist-2] / 3; - jj = list[ilist-2] % 3; - j = MMG5_inxt2[jj]; - pt1 = &mesh->tria[jel]; - pt1->tag[j] |= pt->tag[i1]; - pt1->edg[j] = MG_MAX(pt->edg[i1],pt1->edg[j]); - adja = &mesh->adja[3*(iel-1)+1]; - if ( adja[i1] ) { - kel = adja[i1] / 3; - k = adja[i1] % 3; - mesh->adja[3*(kel-1)+1+k] = 3*jel + j; - mesh->adja[3*(jel-1)+1+j] = 3*kel + k; - pt2 = &mesh->tria[kel]; - pt2->tag[k] |= pt1->tag[j]; - pt2->edg[k] = MG_MAX(pt1->edg[j],pt2->edg[k]); - } - else - mesh->adja[3*(jel-1)+1+j] = 0; - } MMGS_delPt(mesh,ip1); - if ( !MMGS_delElt(mesh,list[0] / 3) ) return 0; - if ( !open ) { - if ( !MMGS_delElt(mesh,list[ilist-1] / 3) ) return 0; + + for ( k=0; ki2 */ +#warning does it still work??? int litcol(MMG5_pMesh mesh,int k,char i,double kali) { MMG5_pTria pt,pt0,pt1; MMG5_pPoint p1,p2; double kal,ps,cosnold,cosnnew,n0old[3],n0new[3],n1old[3],n1new[3],n00old[3],n00new[3]; - int *adja,list[MMGS_LMAX+2],jel,ip2,l,ilist; + int *adja,list[MMGS_LMAX+2],jel,ip2,l,ilist,ishell; char i1,i2,j,jj,j2,open; pt0 = &mesh->tria[0]; @@ -477,7 +688,10 @@ int litcol(MMG5_pMesh mesh,int k,char i,double kali) { /* collect all triangles around vertex i1 */ if ( pt->v[i1] & MG_NOM ) return 0; - ilist = boulet(mesh,k,i1,list); + ilist = MMGS_boulet(mesh,k,i1,MMG5_inxt2[i1],list,&ishell); + if ( ilist <= 1 ) { + return 0; + } /* check for open ball */ adja = &mesh->adja[3*(k-1)+1]; @@ -547,7 +761,7 @@ int litcol(MMG5_pMesh mesh,int k,char i,double kali) { if ( abs(pt->ref) != abs(pt1->ref) ) return 0; } - return colver(mesh,list,ilist); + return colver(mesh,list,ishell,ilist); } /* specific test: no collapse if any interior edge is EDG */ @@ -580,6 +794,3 @@ int litcol(MMG5_pMesh mesh,int k,char i,double kali) { return 0; } - - - diff --git a/src/mmgs/mmgs.h b/src/mmgs/mmgs.h index 4539fcb16..a2fe0918b 100644 --- a/src/mmgs/mmgs.h +++ b/src/mmgs/mmgs.h @@ -135,8 +135,8 @@ int curvpo(MMG5_pMesh ,MMG5_pSol ); int MMG5_mmgs1(MMG5_pMesh ,MMG5_pSol,int* ); int MMGS_mmgs2(MMG5_pMesh ,MMG5_pSol, MMG5_pSol); int MMGS_bdryUpdate(MMG5_pMesh mesh); -int boulet(MMG5_pMesh mesh,int start,int ip,int *list); -int boulechknm(MMG5_pMesh mesh,int start,int ip,int *list); +int MMGS_boulet(MMG5_pMesh mesh,int start,int ip,int ia,int *list,int *ishell); +int MMGS_boulechknm(MMG5_pMesh mesh,int start,int ip,int ia,int *list,int *ishell); int bouletrid(MMG5_pMesh mesh,int start,int ip,int *il1,int *l1,int *il2,int *l2,int *ip0,int *ip1); int MMGS_newPt(MMG5_pMesh mesh,double c[3],double n[3]); void MMGS_delPt(MMG5_pMesh mesh,int ip); @@ -153,8 +153,8 @@ int MMGS_split1(MMG5_pMesh mesh,MMG5_pSol met,int k,int i,int *vx); int MMGS_split2(MMG5_pMesh mesh,MMG5_pSol met,int k,int *vx); int MMGS_split3(MMG5_pMesh mesh,MMG5_pSol met,int k,int *vx); int split1b(MMG5_pMesh mesh,int k,char i,int ip); -int chkcol(MMG5_pMesh mesh,MMG5_pSol met,int k,char i,int *list,char typchk); -int colver(MMG5_pMesh mesh,int *list,int ilist); +int MMGS_chkcol(MMG5_pMesh mesh,MMG5_pSol met,int k,char i,int *list,int *ishell,char typchk); +int colver(MMG5_pMesh mesh,int *list,int ishell,int ilist); int colver3(MMG5_pMesh mesh,int*list); int colver2(MMG5_pMesh mesh,int *ilist); int swapar(MMG5_pMesh mesh,int k,int i); diff --git a/src/mmgs/mmgs1.c b/src/mmgs/mmgs1.c index 677d9a76f..92b7c14f6 100755 --- a/src/mmgs/mmgs1.c +++ b/src/mmgs/mmgs1.c @@ -423,7 +423,7 @@ static int swpmsh(MMG5_pMesh mesh,MMG5_pSol met,char typchk) { static int movtri(MMG5_pMesh mesh,MMG5_pSol met,int maxit) { MMG5_pTria pt; MMG5_pPoint ppt; - int it,k,ier,base,nm,ns,nnm,list[MMGS_LMAX+2],ilist; + int it,k,ier,base,nm,ns,nnm,list[MMGS_LMAX+2],ilist,ishell; char i; if ( abs(mesh->info.imprim) > 5 || mesh->info.ddebug ) @@ -445,7 +445,7 @@ static int movtri(MMG5_pMesh mesh,MMG5_pSol met,int maxit) { if ( ppt->flag == base || MS_SIN(ppt->tag) || ppt->tag & MG_NOM ) continue; - ilist = boulet(mesh,k,i,list); + ilist = MMGS_boulet(mesh,k,i,MMG5_iprv2[i],list,&ishell); if ( !ilist ) continue; @@ -959,7 +959,7 @@ static int colelt(MMG5_pMesh mesh,MMG5_pSol met,char typchk) { MMG5_pPoint p1,p2; MMG5_pPar par; double ll,ux,uy,uz,hmin; - int list[MMGS_LMAX+2],ilist,k,nc,l,isloc,ier; + int list[MMGS_LMAX+2],ilist,k,nc,l,isloc,ier,ishell; char i,i1,i2; nc = 0; @@ -1019,23 +1019,12 @@ static int colelt(MMG5_pMesh mesh,MMG5_pSol met,char typchk) { } /* check if geometry preserved */ - ilist = chkcol(mesh,met,k,i,list,typchk); - if ( ilist > 3 ) { - ier = colver(mesh,list,ilist); - if ( !ier ) return -1; - nc += ier; - break; - } - else if ( ilist == 3 ) { - ier = colver3(mesh,list); - if ( !ier ) return -1; - nc += ier; - break; - } - else if ( ilist == 2 ) { - ier = colver2(mesh,list); + ilist = MMGS_chkcol(mesh,met,k,i,list,&ishell,typchk); + + if ( ilist > 1 ) { + ier = colver(mesh,list,ishell,ilist); + nc += ier; if ( !ier ) return -1; - nc += ier; break; } } @@ -1134,7 +1123,7 @@ static int adpcol(MMG5_pMesh mesh,MMG5_pSol met) { MMG5_pTria pt; MMG5_pPoint p1,p2; double len; - int k,list[MMGS_LMAX+2],ilist,nc,ier; + int k,list[MMGS_LMAX+2],ilist,nc,ier,ishell; char i,i1,i2; nc = 0; @@ -1171,25 +1160,13 @@ static int adpcol(MMG5_pMesh mesh,MMG5_pSol met) { if ( p1->tag > p2->tag || p1->tag > pt->tag[i] ) continue; /* check if geometry preserved */ - ilist = chkcol(mesh,met,k,i,list,2); - if ( ilist > 3 ) { - ier = colver(mesh,list,ilist); + ilist = MMGS_chkcol(mesh,met,k,i,list,&ishell,2); + if ( ilist > 1 ) { + ier = colver(mesh,list,ishell,ilist); nc += ier; if ( !ier ) return -1; break; } - else if ( ilist == 3 ) { - ier = colver3(mesh,list); - nc += ier; - if ( !ier ) return -1; - break; - } - else if ( ilist == 2 ) { - ier = colver2(mesh,list); - nc += ier; - if ( !ier ) return -1; - break; - } } } return nc; From 5d055b1facd41af82174a0570ad8ca3bf1affcdd Mon Sep 17 00:00:00 2001 From: Algiane Date: Fri, 7 Feb 2020 21:36:01 +0100 Subject: [PATCH 09/10] Debug mmgs collapse. --- src/mmgs/colver_s.c | 32 +++++++++++++++++++------------- src/mmgs/mmgs1.c | 3 ++- 2 files changed, 21 insertions(+), 14 deletions(-) diff --git a/src/mmgs/colver_s.c b/src/mmgs/colver_s.c index b329180b2..f3a4b87ee 100644 --- a/src/mmgs/colver_s.c +++ b/src/mmgs/colver_s.c @@ -148,12 +148,14 @@ int MMGS_chkcol(MMG5_pMesh mesh,MMG5_pSol met,int start,char istart,int *list, /* check for open ball */ adja = &mesh->adja[3*(start-1)+1]; - // open = adja[istart] == 0; /** Starting from each tria of the shell of i, travel by adjacency around \a * ip1 and process the portion of surface until meeting a non-manifold edge */ base = ++mesh->base; + + sumilist = 0; + for ( kk = 0; kk<*ishell; ++ kk ) { k = list[kk]/3; i = list[kk]%3; @@ -178,7 +180,6 @@ int MMGS_chkcol(MMG5_pMesh mesh,MMG5_pSol met,int start,char istart,int *list, fwd = 1; if ( pt->tag[i] & MG_NOM ) { - //assert ( !open ); assert ( *ishell > 2 ); /** Count the number of tria that we can reach without crossing a nm edge */ ilist = MMGS_count_bouleploc( mesh, k, i, ip1,&fwd); @@ -198,7 +199,7 @@ int MMGS_chkcol(MMG5_pMesh mesh,MMG5_pSol met,int start,char istart,int *list, /* check references: */ if ( MG_EDG(pt->tag[i2]) ) { - jel = list[1] / 3; + jel = adja[i2] / 3; pt1 = &mesh->tria[jel]; if ( abs(pt->ref) != abs(pt1->ref) ) return 0; } @@ -230,7 +231,7 @@ int MMGS_chkcol(MMG5_pMesh mesh,MMG5_pSol met,int start,char istart,int *list, /* check length */ if ( typchk == 2 && met->m && !MG_EDG(mesh->point[ip2].tag) ) { - ip1 = pt1->v[j2]; + ip1 = pt1->v[j]; len = MMG5_lenSurfEdg(mesh,met,ip1,ip2,0); if ( len > lon || !len ) return 0; } @@ -299,24 +300,29 @@ int MMGS_chkcol(MMG5_pMesh mesh,MMG5_pSol met,int start,char istart,int *list, } /* check angle between 1st and last triangles */ - if ( (!open) && !(pt->tag[i] & MG_GEO) ) { - cosnold = n00old[0]*n1old[0] + n00old[1]*n1old[1] + n00old[2]*n1old[2]; - cosnnew = n00new[0]*n1new[0] + n00new[1]*n1new[1] + n00new[2]*n1new[2]; - if ( cosnold < MMG5_ANGEDG ) { - if ( cosnnew < cosnold ) return 0; + if ( !open ) { + + if ( !(pt->tag[i] & MG_GEO) ) { + cosnold = n00old[0]*n1old[0] + n00old[1]*n1old[1] + n00old[2]*n1old[2]; + cosnnew = n00new[0]*n1new[0] + n00new[1]*n1new[1] + n00new[2]*n1new[2]; + if ( cosnold < MMG5_ANGEDG ) { + if ( cosnnew < cosnold ) return 0; + } + else if ( cosnnew < MMG5_ANGEDG ) return 0; } - else if ( cosnnew < MMG5_ANGEDG ) return 0; /* other checks for reference collapse */ + MMGS_UPDATE_PIV(j,j1,j2,fwd); + adja = &mesh->adja[3*(jel-1)+1]; jel = adja[j]/3; j = adja[j]%3; - MMGS_UPDATE_PIV(j,j1,j2,-fwd); + MMGS_UPDATE_PIV(j,j1,j2,fwd); pt = &mesh->tria[jel]; - if ( MG_EDG(pt->tag[j2]) ) { + if ( MG_EDG(pt->tag[j]) ) { adja = &mesh->adja[3*(jel-1)+1]; - jel = adja[j2] / 3; + jel = adja[j] / 3; pt1 = &mesh->tria[jel]; if ( abs(pt->ref) != abs(pt1->ref) ) return 0; } diff --git a/src/mmgs/mmgs1.c b/src/mmgs/mmgs1.c index 92b7c14f6..208ae2b14 100755 --- a/src/mmgs/mmgs1.c +++ b/src/mmgs/mmgs1.c @@ -1022,13 +1022,14 @@ static int colelt(MMG5_pMesh mesh,MMG5_pSol met,char typchk) { ilist = MMGS_chkcol(mesh,met,k,i,list,&ishell,typchk); if ( ilist > 1 ) { - ier = colver(mesh,list,ishell,ilist); + ier = colver(mesh,list,ishell,ilist); nc += ier; if ( !ier ) return -1; break; } } } + if ( nc > 0 && (abs(mesh->info.imprim) > 5 || mesh->info.ddebug) ) fprintf(stdout," %8d vertices removed\n",nc); From 128e45f61bc21275af8f9af4e4ef763db1dd268e Mon Sep 17 00:00:00 2001 From: Algiane Froehly Date: Mon, 12 Dec 2022 15:02:54 +0100 Subject: [PATCH 10/10] Add last local modifications (debug of partial implementation in progress ). --- src/mmgs/anisomovpt_s.c | 8 +++++- src/mmgs/boulep_s.c | 62 +++++++++++++++++++++++++++++++++++++++++ src/mmgs/mmgs.h | 1 + src/mmgs/mmgs1.c | 14 ++++++++-- src/mmgs/movpt_s.c | 20 +++++++++++-- 5 files changed, 99 insertions(+), 6 deletions(-) diff --git a/src/mmgs/anisomovpt_s.c b/src/mmgs/anisomovpt_s.c index b89119cbc..322ebdc44 100644 --- a/src/mmgs/anisomovpt_s.c +++ b/src/mmgs/anisomovpt_s.c @@ -52,7 +52,7 @@ int movintpt_ani(MMG5_pMesh mesh,MMG5_pSol met,int *list,int ilist) { double r[3][3],ux,uy,uz,*n,area,lispoi[3*MMGS_LMAX+1],*m0;//,m[6],mo[6]; double gv[2],detloc,step,lambda[3],o[3],no[3],to[3],uv[2]; double calold,calnew,caltmp; - int k,iel,kel,nump,nbeg,nend; + int k,iel,kel,nump,nbeg,nend,*adja; char i0,i1,i2,ier; static int warn=0; step = 0.1; @@ -69,6 +69,12 @@ int movintpt_ani(MMG5_pMesh mesh,MMG5_pSol met,int *list,int ilist) { m0 = &met->m[6*nump]; assert( !p0->tag ); + /* Check that we have a closed manifold ball */ + assert ( !(p0->tag & MG_NOM) ); + assert ( !(p0->tag & MG_GEO) ); + assert ( mesh->adja[ 3*(iel-1) + 1 + i1] ); + assert ( mesh->adja[ 3*(iel-1) + 1 + MMG5_iprv2[i0]] ); + iel = list[ilist-1] / 3; i0 = list[ilist-1] % 3; i2 = MMG5_iprv2[i0]; diff --git a/src/mmgs/boulep_s.c b/src/mmgs/boulep_s.c index dbd368c29..d6563a410 100644 --- a/src/mmgs/boulep_s.c +++ b/src/mmgs/boulep_s.c @@ -147,6 +147,68 @@ int MMGS_boulet(MMG5_pMesh mesh,int start,int ip,int ia,int *list, int *ishell) return ilist; } +/** + * \param mesh pointer toward the mesh structure. + * \param start index of triangle to start. + * \param ip index of point for wich we compute the ball. + * \param list pointer toward the computed ball of \a ip. + * \return the size of the computed ball or 0 if fail. + * + * Find all triangles sharing \a ip with ip a manifold point, \f$list[0] =\f$ \a + * start do not stop when crossing ridge. + * + */ +int MMGS_bouletmani(MMG5_pMesh mesh,int start,int ip,int *list) { + MMG5_pTria pt; + MMG5_pPoint ppt; + int *adja,k,ilist; + char i,i1,i2; + + pt = &mesh->tria[start]; + + ppt = &mesh->point[pt->v[ip]]; + ilist = 0; + + /* store neighbors */ + k = start; + i = ip; + do { + if ( ilist > MMGS_LMAX-2 ) return 0; + list[ilist] = 3*k + i; + ++ilist; + + adja = &mesh->adja[3*(k-1)+1]; + i1 = MMG5_inxt2[i]; + k = adja[i1] / 3; + i = adja[i1] % 3; + i = MMG5_inxt2[i]; + } + while ( k && k != start ); + if ( k > 0 ) return ilist; + + if ( ppt->tag & MG_NOM ) + return 0; + + /* check if boundary hit */ + k = start; + i = ip; + do { + adja = &mesh->adja[3*(k-1)+1]; + i2 = MMG5_iprv2[i]; + k = adja[i2] / 3; + if ( k == 0 ) break; + i = adja[i2] % 3; + i = MMG5_iprv2[i]; + + if ( ilist > MMGS_LMAX-2 ) return 0; + list[ilist] = 3*k + i; + ilist++; + } + while ( k ); + + return ilist; +} + /** * \param mesh pointer toward the mesh structure. * \param start index of tria to start to compute the ball. diff --git a/src/mmgs/mmgs.h b/src/mmgs/mmgs.h index a2fe0918b..79baa29fe 100644 --- a/src/mmgs/mmgs.h +++ b/src/mmgs/mmgs.h @@ -136,6 +136,7 @@ int MMG5_mmgs1(MMG5_pMesh ,MMG5_pSol,int* ); int MMGS_mmgs2(MMG5_pMesh ,MMG5_pSol, MMG5_pSol); int MMGS_bdryUpdate(MMG5_pMesh mesh); int MMGS_boulet(MMG5_pMesh mesh,int start,int ip,int ia,int *list,int *ishell); +int MMGS_bouletmani(MMG5_pMesh mesh,int start,int ip,int *list); int MMGS_boulechknm(MMG5_pMesh mesh,int start,int ip,int ia,int *list,int *ishell); int bouletrid(MMG5_pMesh mesh,int start,int ip,int *il1,int *l1,int *il2,int *l2,int *ip0,int *ip1); int MMGS_newPt(MMG5_pMesh mesh,double c[3],double n[3]); diff --git a/src/mmgs/mmgs1.c b/src/mmgs/mmgs1.c index 208ae2b14..868eb24f3 100755 --- a/src/mmgs/mmgs1.c +++ b/src/mmgs/mmgs1.c @@ -443,9 +443,16 @@ static int movtri(MMG5_pMesh mesh,MMG5_pSol met,int maxit) { for (i=0; i<3; i++) { ppt = &mesh->point[pt->v[i]]; - if ( ppt->flag == base || MS_SIN(ppt->tag) || ppt->tag & MG_NOM ) + if ( ppt->flag == base || MS_SIN(ppt->tag) ) { continue; - ilist = MMGS_boulet(mesh,k,i,MMG5_iprv2[i],list,&ishell); + } + else if ( ppt->tag & MG_NOM ) { + ilist = MMGS_boulet(mesh,k,i,MMG5_iprv2[i],list,&ishell); + continue; + } + else { + ilist = MMGS_bouletmani(mesh,k,i,list); + } if ( !ilist ) continue; @@ -453,8 +460,9 @@ static int movtri(MMG5_pMesh mesh,MMG5_pSol met,int maxit) { ier = movridpt(mesh,met,list,ilist); if ( ier ) ns++; } - else + else { ier = movintpt(mesh,met,list,ilist); + } if ( ier ) { nm++; ppt->flag = base; diff --git a/src/mmgs/movpt_s.c b/src/mmgs/movpt_s.c index 1e5eb5d95..9a25fa636 100644 --- a/src/mmgs/movpt_s.c +++ b/src/mmgs/movpt_s.c @@ -37,14 +37,24 @@ #include -/* compute movement of an internal point whose ball is passed */ +/** + * \param mesh pointer toward the mesh structure. + * \param met pointer toward the metric structure. + * \param list pointer toward the ball of point. + * \param ilist number of elements in the ball of point. + * + * \return 0 if impossible to move the point, 1 otherwise. + * + * Compute movement of an internal point whose ball is passed. + * + */ int movintpt_iso(MMG5_pMesh mesh,MMG5_pSol met,int *list,int ilist) { MMG5_pPoint p0,p1,ppt0; MMG5_pTria pt,pt0; MMG5_Bezier b; double aa,bb,ab,ll,l,mlon,devmean,GV[3],gv[2],cosalpha,sinalpha,r[3][3],*n,lispoi[3*MMGS_LMAX+1]; double ux,uy,uz,det2d,detloc,step,lambda[3],uv[2],o[3],no[3],to[3],Vold,Vnew,calold,calnew,caltmp; - int ier,iel,ipp,k,kel,npt,ibeg,iend; + int ier,iel,ipp,k,kel,npt,ibeg,iend,*adja; char i0,i1,i2; step = 0.1; @@ -59,6 +69,12 @@ int movintpt_iso(MMG5_pMesh mesh,MMG5_pSol met,int *list,int ilist) { ipp = pt->v[i0]; p0 = &mesh->point[ipp]; /* point to move */ + /* Check that we have a closed manifold ball */ + assert ( !(p0->tag & MG_NOM) ); + assert ( !(p0->tag & MG_GEO) ); + assert ( mesh->adja[ 3*(k-1) + 1 + i1] ); + assert ( mesh->adja[ 3*(k-1) + 1 + MMG5_iprv2[i0]] ); + k = list[ilist-1] / 3; i0 = list[ilist-1] % 3; i1 = MMG5_inxt2[i0];