gibbs.cpp

FreeFem++可以生成高质量的有限元网格。可以用于流体力学

	bnew = mmax(i__2,i__3);
/* L110: */
    }
/*      if(impre.ne.0) then */
/*        write(nfout,*)'profile  old  = ',pfold,', profile  new = ',pfnew
 */
/*        write(nfout,*)'1/2 bande old = ',bold ,', 1/2 band new = ',bnew 
*/
/*      endif */
return 0;
} /* gibbsa_ */

/* Subroutine */ int gibbsb_(integer* x,integer*  y,integer*  n,integer*  ptvois,
integer*  vois,integer*  nx,integer*  ny,integer*  nv,integer*  nn,integer*  m,
integer*  wh,integer*  wl,integer* r, integer* , integer* )
{
    /* System generated locals */
    integer i__1, i__2;

    /* Local variables */
    static int flag_;
    static integer i, j, k, p, s, h0, i1, l0, i2;
/*    extern  Subroutine  int gibbs1_(); */
    static integer lg;
/*    extern  Subroutine  int gibbsd_(), gibbsc_();*/
    static integer niveau, mxcanx, mxcany, nbc;

/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 */
/* ...................................................................... 
*/
/*     attention on met la descente optimiser dans r <0 ou nulle */
/* .......................................................................
 */
    /* Parameter adjustments */
    --r;
    --m;
    --ny;
    --nx;
    --vois;
    --ptvois;

    /* Function Body */
    gibbsd_(y, n, &ptvois[1], &vois[1], nv, &r[1], &niveau);
    gibbsc_(&ny[1], nv, &niveau, n, &mxcany);
    gibbsd_(x, n, &ptvois[1], &vois[1], nv, &r[1], &niveau);
    p = niveau;
    gibbsc_(&nx[1], nv, &niveau, n, &mxcanx);
    flag_ = ffalse;
    i__1 = *n;
    for (i = 1; i <= i__1; ++i) {
	if (nx[i] + ny[i] == p) {
	    r[i] = -nx[i];
	} else if (nx[i] >= 0) {
	    flag_ = ttrue;
	    r[i] = -1073741824;
	} else {
	    if (r[i] <= 0) {
		r[i] = -1073741822;
	    }
	}
/* L20: */
    }
    if (flag_) {
/*       calcul des composantes connexe du graphe sans les sommets de 
nn */
/*       ------------------------------------------------------------
--- */
	j = *n;
	k = 0;
	nbc = 0;
	nv[nbc] = j;
L30:
	++k;
	if (k <= *n) {
	    if (r[k] == -1073741824) {
/*           recherche de la fermeture transitive partant de k
 */
		++nbc;
		i = -1;
		s = k;
L40:
		++i;
		wl[i] = ptvois[s];
		wh[i] = ptvois[s + 1];
		++j;
		nv[j] = s;
		r[s] = -1073741823;
L50:
		if (i >= 0) {
		    if (wl[i] < wh[i]) {
			s = vois[wl[i]];
			++wl[i];
			if (r[s] == -1073741824) {
			    goto L40;
			}
			goto L50;
		    }
		    --i;
		    goto L50;
		}
		nv[nbc] = j;
		m[nbc] = nbc;
	    }
	    goto L30;
	}
/*        if(impre.lt.0) write(nfout,*) */
/*     +         ' nb de composante connexe du graphe reduit =',nbc */

/* --------------- fin de construction des composantes connexes------
--- */
/*        nv(0)=n */
/*        if(impre.le.-10) write(nfout,5555)'nv(0:n+n) = ',(nv(i),i=0,
n+n) */
	gibbs1_(&nbc, &m[1], nv);
/*        if(impre.le.-10)write(nfout,5555)'trie m =',(m(i),i=1,nbc) 
*/
	i__1 = p;
	for (i = 0; i <= i__1; ++i) {
	    nn[i] = 0;
/* L60: */
	}
	i__1 = *n;
	for (i = 1; i <= i__1; ++i) {
	    j = -r[i];
	    if (j >= 0 && j <= p) {
		++nn[j];
	    }
/* L70: */
	}

/*       boucle sur les composante connexes par ordre croissantes */
/*       -------------------------------------------------------- */
	for (k = nbc; k >= 1; --k) {
	    i = m[k];
	    i1 = nv[i - 1] + 1;
	    i2 = nv[i];
	    lg = i2 - i1 + 1;
/*         if(impre.le.-7) */
/*     +       write(nfout,*) k,' composante ',i,',lg=',lg,',i1,i2
=',i1,i2 */
/*         if(impre.le.-8) */
/*     +       write (nfout,5555)' ',(nv(i),i=i1,i2) */
	    h0 = 0;
	    l0 = 0;
	    i__1 = p;
	    for (j = 0; j <= i__1; ++j) {
		wh[j] = nn[j];
		wl[j] = nn[j];
/* L90: */
	    }
	    i__1 = i2;
	    for (i = i1; i <= i__1; ++i) {
		s = nv[i];
		++wh[nx[s]];
		++wl[p - ny[s]];
/* L100: */
	    }
	    i__1 = p;
	    for (j = 0; j <= i__1; ++j) {
		if (wh[j] != nn[j]) {
/* Computing MAX */
		    i__2 = wh[j];
		    h0 = mmax(i__2,h0);
		}
		if (wl[j] != nn[j]) {
/* Computing MAX */
		    i__2 = wl[j];
		    l0 = mmax(i__2,l0);
		}
/* L110: */
	    }
	    if (h0 < l0 || h0 == l0 && mxcanx <= mxcany) {
/*           if(impre.le.-2) write(nfout,*) */
/*     +       '         h0 = ',h0,',l0 = ',l0,'  ------- XXXX
 --------' */
		i__1 = i2;
		for (i = i1; i <= i__1; ++i) {
		    s = nv[i];
		    r[s] = -nx[s];
		    ++nn[-r[s]];
/* L120: */
		}
	    } else {
/*           if (impre.le.-2) write(nfout,*) */
/*     +       '         h0 = ',h0,',l0 = ',l0,'  ------- YYYY
 --------' */
		i__1 = i2;
		for (i = i1; i <= i__1; ++i) {
		    s = nv[i];
		    r[s] = -p + ny[s];
		    ++nn[-r[s]];
/* L130: */
		}
	    }
/* L140: */
	}
    }
/*     on met les nouveaux niveaux de la descendance optimiser dans nn */
/*     ----------------------------------------------------------------- 
*/
    i__1 = *n;
    for (i = 1; i <= i__1; ++i) {
	if (r[i] > 0) {
	    nn[i] = -1;
	} else if (r[i] == -1073741822) {
	    nn[i] = -2;
	} else {
	    nn[i] = -r[i];
	}
/* L150: */
    }
/*      if(impre.le.-10)write (nfout,5555)' nn(i)=',(nn(i),i=1,n) */
/* 5555  format('            --------   ',a,/,5(15x,10(i5)/)) */
return 0;} /* gibbsb_ */

/* Subroutine */ int gibbsc_(integer* nz,integer*  nv,integer*  niveau,integer*  n,integer*  mxz)
{
    /* System generated locals */
    integer i__1, i__2, i__3;

    /* Local variables */
    static integer i, j;

/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 */
    /* Parameter adjustments */
    --nz;

    /* Function Body */
    i__1 = *n;
    for (i = 1; i <= i__1; ++i) {
	nz[i] = -1;
/* L10: */
    }
    *mxz = 0;
    i__1 = *niveau;
    for (i = 0; i <= i__1; ++i) {
/* Computing MAX */
	i__2 = *mxz, i__3 = nv[i + 1] - nv[i];
	*mxz = mmax(i__2,i__3);
	i__2 = nv[i + 1];
	for (j = nv[i] + 1; j <= i__2; ++j) {
	    nz[nv[j]] = i;
/* L20: */
	}
    }
return 0;} /* gibbsc_ */

/* Subroutine */ int gibbsd_(integer* racine,integer*  n,integer*  ptvois,integer*  vois,integer*  nv,integer*  r,integer*  niveau)
{
    /* System generated locals */
    integer i__1, i__2;

    /* Local variables */
    static integer i, k, s, sv, stk, stk1, stk2;

/* -----------------------------------------------------------------------
 */
/*     but construire la structure des descendant de racine  du graphe */
/* -----------------------------------------------------------------------
 */
/*     sortie : */
/*     -------- */
/*     nv est la structure des niveaux */
/*     les sommets du niveau (i =0,niveau_ sont defini par : */
/*        (nv(j),j=nv(i),nv(i+1)-1) */

/*     le tableau r(i) n'est modifier que sur les sommets */
/*       de la composante connexe du graphe contenant la racine */

/* -----------------------------------------------------------------------
 */

/*     on demark tout les sommets non remuneroter */
/* -------------------------------------------------- */
    /* Parameter adjustments */
    --r;
    --vois;
    --ptvois;

    /* Function Body */
    i__1 = *n;
    for (i = 1; i <= i__1; ++i) {
	if (r[i] < 0) {
	    r[i] = 0;
	}
/* L10: */
    }

/*    initialisation */

    stk = *n - 1;
    nv[0] = stk;
    stk2 = stk;
    *niveau = 0;
    ++stk;
    nv[stk] = *racine;
    r[*racine] = -1;
L20:
    if (stk2 < stk) {
	++(*niveau);
	stk1 = stk2 + 1;
	nv[*niveau] = stk;
	stk2 = stk;
/*        print *,' ------- niveau =',niveau,' stk=',stk1,stk2 */
	i__1 = stk2;
	for (k = stk1; k <= i__1; ++k) {
	    s = nv[k];
/*         print *,'----------------- s=',s */
	    i__2 = ptvois[s + 1] - 1;
	    for (i = ptvois[s]; i <= i__2; ++i) {
/*               pour tout les sommets (sv) voisin */
/*                d'un sommet (s) du niveau precedent */
		sv = vois[i];
/*          print *,' voisin =',sv */
/*               si le sommet n'est pas marque on le marque et
 on l'ajout */
		if (r[sv] == 0) {
		    ++stk;
		    nv[stk] = sv;
		    r[sv] = -1;
		}
/* L30: */
	    }
/* L40: */
	}
	goto L20;
    }
    --(*niveau);
/*      call pnv(' gibbsd ',n,nv,niveau) */
return 0;} /* gibbsd_ */


/* Subroutine */ int gibbst_(integer* n,integer*  p,integer*  nv,integer*  nn,integer*  ptvois,integer*  vois,integer*  m,integer*  r,integer*  new_,integer*  option, 
	integer* pfnew,integer*  impre,integer*  )
{
    /* System generated locals */
    integer i__1, i__2, i__3, i__4, i__5;

    /* Local variables */
    static integer nbsc, bnew, knew, step, plus, i, j, k, s, debut, i1, i2;
/*    extern Subroutine int gibbs2_();*/
    static integer fin;


/*     construction de la stucture de niveau dans nv a partir de nn */
/*     ------------------------------------------------------------ */
    /* Parameter adjustments */
    --r;
    --m;
    --vois;
    --ptvois;

    /* Function Body */
    nv[0] = *n;
    i__1 = *p + 1;
    for (i = 1; i <= i__1; ++i) {
	nv[i] = 0;
/* L150: */
    }
    i__1 = *n;
    for (i = 1; i <= i__1; ++i) {
	if (nn[i] >= 0) {
	    ++nv[nn[i] + 1];
	}
/* L160: */
    }
    i__1 = *p;
    for (i = 0; i <= i__1; ++i) {
	nv[i + 1] += nv[i];
/* L170: */
    }
    i__1 = *n;
    for (i = 1; i <= i__1; ++i) {
	if (nn[i] >= 0) {
	    j = nn[i];
	    ++nv[j];
	    nv[nv[j]] = i;
	}
/* L180: */
    }
    for (i = *p; i >= 0; --i) {
	nv[i + 1] = nv[i];
/* L190: */
    }
    nv[0] = *n;
    nbsc = nv[*p + 1] - nv[0];
/*     --- fin de la construction ------------------------------------ */
    if (*option == -2) {
	i__1 = *impre - 1;
    }
    i__1 = *n;
    for (i = 1; i <= i__1; ++i) {
	m[i] = *n * 3 + ptvois[i + 1] - ptvois[i];
/* L10: */
    }
    if ((((int)*option) == 1)||(((int)*option) == -1)) {
	debut = 0;
	fin = *p;
	step = 1;
    } else {
	debut = *p;
	fin = 0;
	step = -1;
    }
    i__1 = fin;
    i__2 = step;
    for (i = debut; i__2 < 0 ? i >= i__1 : i <= i__1; i += i__2) {
	i1 = nv[i] + 1;