Feature/mmgs nm

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 {

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) ) { */
@@ -462,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;
@@ -520,14 +559,14 @@ 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;
     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;

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];

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] */