femgibbs.cpp

FreeFEM is an implementation of the GFEM language dedicated to the finite element method. It provid

/*       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 
femMesh::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 
femMesh::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 
femMesh::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 * nfout)
{
  /* 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 (abs ((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;
       i2 = nv[i + 1];
       i__3 = i2 - i1 + 1;
       gibbs2_ (&i__3, &nv[i1], &m[1]);
       i__3 = i2;
       for (j = i1; j <= i__3; ++j)
	  {
	    s = nv[j];
	    i__4 = ptvois[s + 1] - 1;
	    for (k = ptvois[s]; k <= i__4; ++k)
	       {
/* Computing MIN */
		 i__5 = m[vois[k]];
		 m[vois[k]] = mmin (i__5, j);
/* L20: */
	       }
/* L30: */
	  }
/* L40: */
     }
  if (*option > 0)
     {
       knew = *new_;
       plus = 1;
     }
  else
     {
       knew = *new_ + nbsc + 1;
       plus = -1;
     }
  *new_ += nbsc;
/*      if(option.gt.0) then */
/*        do 60 k = debut , fin , step */
/*          do 60 j = nv(k+1),nv(k)+1,-1 */
/*            knew = knew + plus */
/*            r(nv(j)) = knew */
/* 60      continue */
/*      else */
  i__2 = fin;
  i__1 = step;
  for (k = debut; i__1 < 0 ? k >= i__2 : k <= i__2; k += i__1)
     {
       i__3 = nv[k + 1];
       for (j = nv[k] + 1; j <= i__3; ++j)
	  {
	    knew += plus;
	    r[nv[j]] = knew;
/* L70: */
	  }
     }
/*      endif */
  *pfnew = 0;
  bnew = 0;
  i__3 = *n;
  for (i = 1; i <= i__3; ++i)
     {
       k = r[i];
       if (k > 0)
	  {
	    i__1 = ptvois[i + 1] - 1;
	    for (j = ptvois[i]; j <= i__1; ++j)
	       {
		 if (r[vois[j]] > 0)
		    {
/* Computing MIN */
		      i__2 = k, i__4 = r[vois[j]];
		      k = mmin (i__2, i__4);
		    }
/* L100: */
	       }
	    *pfnew = *pfnew + r[i] - k + 1;
/* Computing MAX */
	    i__1 = bnew, i__2 = r[i] - k + 1;
	    bnew = mmax (i__1, i__2);
	  }
/* L110: */
     }
/*      if(impre.lt.0.or.impre.gt.2) then */
/*        write(nfout,*) '      option =',option,', profile =',pfnew */
/*     +       ,', 1/2 bande =',bnew,', new=',new,', nbss composante=',nbsc
 */
/*      endif */
  return 0;
}				/* gibbst_ */

/* function */
int 
femMesh::gibbsv (integer * ptvoi,integer * vois, integer * lvois, integer * w, integer * v)
{
  /* System generated locals */
  integer         i__2;

  /* Local variables */
  integer         i, j, k, T, ss, iii, ptv, ptv1;
  integer         nbss = np, nbt = nt;

/*--- Prepare les donees pour gibbsa en construisant ptvoi, vois, lvois -
------------*/
/*     in */
/*     ---   nbnt =3 pour des triangles 2D, */
/*      nbt =  nb de triangle */
/*        nbss = nb de sommets */
/*           nsea = numeros de 3 sommets de chaque triangle (me) */
/*     out */
/*     ---  ptvoi, vois, lvois, err */
/*      tableaux de travail w, v */
/*-----------------------------------------------------------------------
----------*/
  /* Parameter adjustments */
  --v;
  --w;
  --vois;
  --ptvoi;

  /* Function Body */
  for (i = 1; i <= nbss; ++i)
     {
       w[i] = -1;
       ptvoi[i] = 0;
     }
  ptvoi[nbss + 1] = 0;
  for (i = 0; i < nbt; ++i)
     {
       for (j = 0; j < 3; ++j)
	  {
	    ss = tr[i][j] + 1;
	    ++ptvoi[ss + 1];
	    w[ss] = 0;
	  }
     }

  for (i = 1; i <= nbss; ++i)
    ptvoi[i + 1] += ptvoi[i];

  for (i = 0; i < nbt; ++i)
     {
       for (j = 0; j < 3; ++j)
	  {
	    ss = tr[i][j] + 1;
	    ++ptvoi[ss];
	    v[ptvoi[ss]] = i;
	  }
     }
  ptv1 = 0;
  iii = 1;
  for (i = 1; i <= nbss; ++i)
     {
       ptv = ptv1 + 1;
       ptv1 = ptvoi[i];
       ptvoi[i] = iii;
       i__2 = ptv1;
       for (j = ptv; j <= i__2; ++j)
	  {
	    T = v[j];
	    for (k = 0; k < 3; ++k)
	       {
		 ss = tr[T][k] + 1;	/*  nsea[k + T * nsea_dim1]; */
		 if (w[ss] != i)
		    {
		      w[ss] = i;
		      if (iii > *lvois)
			return 2;
		      /* print*,'pas assez de place memoire' */

		      vois[iii] = ss;
		      ++iii;
		    }
	       }
	  }
     }
  ptvoi[nbss + 1] = iii;
  *lvois = iii - 1;
  return 0;			/* OK */
}				/* gibbsv_ */

int
femMesh::renumerotate()
/* -------- 
   renumber vertices by gibbs method
   in:   ns,nt,q,me,ngg  
   out:   q,me,ngg 
   auxiliary arrays: ptvois,vois,r,m,nv,nx,ny,nn,w1,w2,f 
   all of size ns+1 except vois (10(ns+1)) and nv (2(ns+1))
   err = -3: fatal erreur  gibbs 2 : pb racine
 */
{
  long            i, j, pfold, pfnew;
  long           *ptvois = NULL;
  long           *vois = NULL;
  long           *nn = NULL;
  long           *r = NULL;
  long           *m = NULL;
  long           *nv = NULL;
  long           *nx = NULL;
  long           *ny = NULL;
  long           *w1 = NULL;
  long           *w2 = NULL;
  femPoint         *f = NULL;
  long            ns = np;
  long            nbvoisin = 10 * ns;
  long            printint = 0, iodev = 6;
  int             err = 0;
  int            *ngg = ng;

  ptvois = new long[ns+1];
  nn = new long[3*nt];
  vois = new long[nbvoisin+10];
  r = new long[ns+1];

  err = gibbsv (ptvois, vois, &nbvoisin, r, nn);
  delete [] nn;

  if (err == 0)
     {
       m = new long[ns+1];
       nn = new long[ns+1];
       nv = new long[(ns+1)<<1];
       nx = new long[ns+1];
       ny = new long[ns+1];
       w1 = new long[ns+1];
       w2 = new long[ns+1];

       err = gibbsa_ (&ns, ptvois, vois, r, m, nv, nx, ny, nn, w1, w2, &pfold, &pfnew,
		      &printint, &iodev);
       delete [] m;
       delete [] nv;
       delete [] nn;
       delete [] nx;
       delete [] ny;
       delete [] w1;
       delete [] w2;
     }

  delete [] vois;
  if (err == 0 && (pfnew <= pfold))
     {
       for (i = 0; i < ns; ++i)
	 ptvois[i] = ngg[i];
       for (i = 0; i < ns; ++i)
	 ngg[r[i] - 1] = ptvois[i];
     }
  delete [] ptvois;
  if (err == 0 && (pfnew <= pfold))
     {
       f = new femPoint[ns];
       for (i = 0; i < ns; ++i)
	 {
	   f[i][0] = rp[i][0];
	   f[i][1] = rp[i][1];
	 }
       for (i = 0; i < ns; ++i)
	 {
	   rp[r[i] - 1][0] = f[i][0];
	   rp[r[i] - 1][1] = f[i][1];
	 }

       for (j = 0; j < nt; ++j)
	 for (i = 0; i < 3; i++)
	   tr[j][i] = r[tr[j][i]] - 1;

       delete [] f;
     }
  delete [] r;
  return err;
}
/*  message d'erreur:         *err = 2;    print*,'pas assez de place memoire'   */

}