Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Feature/mmgs nm #190

Draft
wants to merge 10 commits into
base: develop
Choose a base branch
from
2 changes: 1 addition & 1 deletion src/common/libmmgtypes.h
Original file line number Diff line number Diff line change
Expand Up @@ -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 {
Expand Down
51 changes: 45 additions & 6 deletions src/mmgs/analys_s.c
Original file line number Diff line number Diff line change
Expand Up @@ -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
*
Expand Down Expand Up @@ -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 */
Expand Down Expand Up @@ -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++;
Expand Down Expand Up @@ -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;

Expand Down Expand Up @@ -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) ) { */
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand Down
8 changes: 7 additions & 1 deletion src/mmgs/anisomovpt_s.c
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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];
Expand Down
18 changes: 9 additions & 9 deletions src/mmgs/anisosiz_s.c
Original file line number Diff line number Diff line change
Expand Up @@ -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];
Expand All @@ -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;
Expand Down Expand Up @@ -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;
Expand All @@ -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] */
Expand Down Expand Up @@ -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;
Expand All @@ -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] */
Expand Down
Loading