mshptg.cpp

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

		 for (i = i1; i <= i_3; ++i)
		    {
		      if (w[i] % 3 == 0)
			 {
			   k = w[i] - 2;
			 }
		      else
			 {
			   k = w[i] + 1;
			 }
		      bx += c[(nu[k] << 1) + 1];
		      by += c[(nu[k] << 1) + 2];
		    }
		 bx /= i2 - i1 + 1;
		 by /= i2 - i1 + 1;
		 depx = omega * (c[(is << 1) + 1] - bx);
		 depy = omega * (c[(is << 1) + 2] - by);
		 c[(is << 1) + 1] -= depx;
		 c[(is << 1) + 2] -= depy;
/* Computing MAX */
		 r_1 = err, r_2 = aabs (depx), r_1 = amax (r_1, r_2), r_2 = aabs (depy);
		 err = amax (r_1, r_2);
	       }
	  }
/* -------------------------------- */
       if (err <= eps * dx)
	  {
	    return 0;
	  }
     }
  return 1;
}				/* mshrgl_ */

/* ********************************************************************** */
int 
mshgpt_ (long *c, double *cr, long *nu, double *h, long *reft, long *nbs,
     long nbsmx, long *nbt, double coef, double puis, double *trfri, long *err)
{
  /* System generated locals */
  long            i_1;
  double           r_1, r_2, r_3;
  double          d_1, d_2, d_3, d_4, d_5, d_6, d_7, d_8;

  /* Local variables */
  static double    aire;
  static long     tete;
  static long     t;
  static double    x, y;
  static long     itera;
  static double    h1, h2, h3;
  static long     s1, s2, s3;
  static double    hs;
  static long     ix, iy, nbsold;
  static double    det, pui;

  /* Parameter adjustments */
  --trfri;
  --reft;
  --h;
  nu -= 7;
  cr -= 3;
  c -= 3;

  /* Function Body */
  pui = puis;
  *nbt = (*nbs << 1) - 2;
  if (*nbs >= nbsmx)
     { /*  ADD FH  avril 2007 */
      i_1 = *nbt;
      for (t = 1; t <= i_1; ++t)
        if (nu[t * 6 + 6] != 0)
        { 
	  mshopt_ (&c[3], &nu[7], &t, 4L, err);
	  mshopt_ (&c[3], &nu[7], &t, 5L, err);
	  mshopt_ (&c[3], &nu[7], &t, 6L, err);
        }	
	/* FIN  ADD FH  vril 2007 */    
       return 0;
     }
  tete = 0;
/*     initialisation de la liste des triangles libre */
  i_1 = *nbt;
  for (t = 1; t <= i_1; ++t)
     {
       if (nu[t * 6 + 6] == 0)
	  {
	    nu[t * 6 + 1] = tete;
	    tete = t;
	  }
     }
  itera = 0;
L20:
  ++itera;
  nbsold = *nbs;
  i_1 = *nbt;
  for (t = 1; t <= i_1; ++t)
     {
       if (nu[t * 6 + 6] != 0)
	  {
	    s1 = nu[t * 6 + 1];
	    s2 = nu[t * 6 + 2];
	    s3 = nu[t * 6 + 3];
/*        calcul de 2 fois l'aire du triangle */
	    det = (cr[(s2 << 1) + 1] - cr[(s1 << 1) + 1]) 
	        * (cr[(s3 << 1) + 2] - cr[(s1 << 1) + 2]) 
	        - (cr[(s2 << 1) + 2] - cr[(s1 << 1) + 2]) 
	        * (cr[(s3 << 1) + 1] - cr[(s1 << 1) + 1]);
	    aire = det * coef;
	    if (puis > (double) 0.)
	       {
		 d_2 = (double) h[s1];
		 d_3 = (double) pui;
		 d_4 = (double) h[s2];
		 d_5 = (double) pui;
		 d_6 = (double) h[s3];
		 d_7 = (double) pui;
		 d_1 = (pow (d_2, d_3) + pow (d_4, d_5) + pow (d_6, d_7)) / (double) 3.;
		 d_8 = (double) (1.F / pui);
		 hs = (double)pow (d_1, d_8);
	       }
	    else if (puis > (double) -1.)
	       {
		 d_1 = (double) (h[s1] * h[s2] * h[s3]);
		 hs = (double)pow (d_1, 1. / 3);
	       }
	    else if (puis > (double) -2.)
	       {
		 hs = h[s1] * (double) 3. *h[s2] * h[s3] / (h[s1] * h[s2] + h[
					       s1] * h[s3] + h[s2] * h[s3]);
	       }
	    else
	       {
/* Computing 2nd power */
		 r_1 = (double)sqrt (h[s1] * h[s2]);
/* Computing 2nd power */
		 r_2 = h[s1] * h[s3];
/* Computing 2nd power */
		 r_3 = h[s2] * h[s3];
		 hs = (double)sqrt (3.0) * (h[s1] * h[s2] * h[s3]
			 / (r_1 * r_1) + r_2 * r_2 + r_3 * r_3);
	       }
	    if (aire > hs * hs)
	       {
		 h1 = (double) 1.;
		 h2 = (double) 1.;
		 h3 = (double) 1.;
		 x = (cr[(s1 << 1) + 1] * h1 + cr[(s2 << 1) + 1] * h2 + cr[(s3
					  << 1) + 1] * h3) / (h1 + h2 + h3);
		 y = (cr[(s1 << 1) + 2] * h1 + cr[(s2 << 1) + 2] * h2 + cr[(s3
					  << 1) + 2] * h3) / (h1 + h2 + h3);
		 r_1 = trfri[1] * (x - trfri[2]);
		 ix = (long) i_nint (r_1);
		 r_1 = trfri[1] * (y - trfri[3]) - trfri[4];
		 iy = (long) i_nint (r_1);
		 if ((c[(s2 << 1) + 1] - ix) * (c[(s3 << 1) + 2] - iy) - (c[(
		       s2 << 1) + 2] - iy) * (c[(s3 << 1) + 1] - ix) <= 0 ||
		     (ix - c[(s1 << 1) + 1]) * (c[(s3 << 1) + 2] - c[(s1 <<
			 1) + 2]) - (iy - c[(s1 << 1) + 2]) * (c[(s3 << 1) +
		      1] - c[(s1 << 1) + 1]) <= 0 || (c[(s2 << 1) + 1] - c[(
			s1 << 1) + 1]) * (iy - c[(s1 << 1) + 2]) - (c[(s2 <<
		       1) + 2] - c[(s1 << 1) + 2]) * (ix - c[(s1 << 1) + 1])
		     <= 0)
		    {
		    }
		 else
		    {
		      if (*nbs >= nbsmx)
			 {
			   return 0;
			 }
		      ++(*nbs);
		      c[(*nbs << 1) + 1] = ix;
		      c[(*nbs << 1) + 2] = iy;
		      cr[(*nbs << 1) + 1] = ix / trfri[1] + trfri[2];
		      cr[(*nbs << 1) + 2] = (iy + trfri[4]) / trfri[1] + trfri[
									 3];
		      h[*nbs] = hs;
		      msha1p_ (&t, nbs, &c[3], &nu[7], &reft[1], &tete, nbt, err);
		      if (*err != 0)
			 {
			   return 0;
			 }
		    }
	       }
	  }
     }
  if (nbsold != *nbs)
     {
       goto L20;
     }
  return 1;
}				/* mshgpt_ */

/* ********************************************************************** */
int 
msha1p_ (long *t, long *s, long *c, long *nu, long *reft, long *tete, long *nbt,
	 long *err)
{
  static long     t1, t2, t3, ia2, ia3;
  static long     ta2, ta3;

//    extern int mshopt_();
  static long     tta;

  /* Parameter adjustments */
  --reft;
  nu -= 7;
  c -= 3;

  /* Function Body */
  t1 = *t;
  if (*tete == 0)
     {
       ++(*nbt);
       t2 = *nbt;
     }
  else
     {
       t2 = *tete;
       *tete = nu[*tete * 6 + 1];
     }
  if (*tete == 0)
     {
       ++(*nbt);
       t3 = *nbt;
     }
  else
     {
       t3 = *tete;
       *tete = nu[*tete * 6 + 1];
     }
  nu[t2 * 6 + 1] = *s;
  nu[t2 * 6 + 2] = nu[*t * 6 + 2];
  nu[t2 * 6 + 3] = nu[*t * 6 + 3];
  nu[t2 * 6 + 4] = (t1 << 3) + 5;
  nu[t2 * 6 + 5] = nu[*t * 6 + 5];
  nu[t2 * 6 + 6] = (t3 << 3) + 5;
  nu[t3 * 6 + 1] = nu[*t * 6 + 1];
  nu[t3 * 6 + 2] = *s;
  nu[t3 * 6 + 3] = nu[*t * 6 + 3];
  nu[t3 * 6 + 4] = (t1 << 3) + 6;
  nu[t3 * 6 + 5] = (t2 << 3) + 6;
  nu[t3 * 6 + 6] = nu[*t * 6 + 6];
  tta = nu[*t * 6 + 5];
  if (tta > 0)
     {
       ta2 = tta / 8;
       ia2 = tta - (ta2 << 3);
       nu[ia2 + ta2 * 6] = (t2 << 3) + 5;
     }
  tta = nu[*t * 6 + 6];
  if (tta > 0)
     {
       ta3 = tta / 8;
       ia3 = tta - (ta3 << 3);
       nu[ia3 + ta3 * 6] = (t3 << 3) + 6;
     }
  nu[t1 * 6 + 3] = *s;
  nu[t1 * 6 + 5] = (t2 << 3) + 4;
  nu[t1 * 6 + 6] = (t3 << 3) + 4;
  reft[t2] = reft[*t];
  reft[t3] = reft[*t];
  mshopt_ (&c[3], &nu[7], &t1, 4L, err);
  if (*err != 0)
     {
       return 0;
     }
  mshopt_ (&c[3], &nu[7], &t2, 5L, err);
  if (*err != 0)
     {
       return 0;
     }
  mshopt_ (&c[3], &nu[7], &t3, 6L, err);
  if (*err != 0)
     {
       return 0;
     }
  return 1;
}				/* msha1p_ */

/* ********************************************************************** */
long 
mshlcl_ (long *c, long *nu, long *tete, long *s)
{
  /* System generated locals */
  long            ret_val;

  /* Local variables */
  static long     init;
  static long     x, y, pt, ppt;
  LONG8 det;

  /* Parameter adjustments */
  nu -= 7;
  c -= 3;

  /* Function Body */
  x = c[(*s << 1) + 1];
  y = c[(*s << 1) + 2];
  init = 1;
  pt = *tete;
L10:
  ppt = pt;
  pt = nu[pt * 6 + 4];
  if (pt != *tete)
     {
       det = (LONG8) x * (LONG8) c[(nu[pt * 6 + 1] << 1) + 2] 
           - (LONG8) y * (LONG8) c[(nu[pt * 6 + 1] << 1) + 1];
       if (det < 0)
	  {
	    init = 0;
	    goto L10;
	  }
       else if (init && det == 0)
	  {
	    goto L10;
	  }
     }
  ret_val = ppt;
  return ret_val;
}				/* mshlcl_ */

/* ********************************************************************** */
int 
mshtri_ (double *cr, long *c, long *nbs, long *tri, LONG8  *nu, double *trfri, long *err)
{
  /* System generated locals */
  long            i_1, i_2, i_3;
  double           r_1;

  /* Local variables */
  static long     ierr, trik;
  static double    xmin, ymin, xmax, ymax;
  static long     i, j, k, ic, jc;

//    extern int mshtr1_();
  static long     ip, xx;
  static double    aa1, aa2;
  static long     iii, tri3;
  LONG8 det;

  /* Parameter adjustments */
  --trfri;
  --nu;
  --tri;
  c -= 3;
  cr -= 3;

  /* Function Body */
  ierr = 0;
  iii = 1;
  xmin = cr[3];
  ymin = cr[4];
  xmax = cr[3];
  ymax = cr[4];
  i_1 = *nbs;
  for (ic = 1; ic <= i_1; ++ic)
     {
/* Computing MIN */
       r_1 = cr[(ic << 1) + 1];
       xmin = amin (r_1, xmin);
/* Computing MIN */
       r_1 = cr[(ic << 1) + 2];
       ymin = amin (r_1, ymin);
/* Computing MAX */
       r_1 = cr[(ic << 1) + 1];
       xmax = amax (r_1, xmax);
/* Computing MAX */
       r_1 = cr[(ic << 1) + 2];
       ymax = amax (r_1, ymax);
       tri[ic] = ic;
       if (cr[(ic << 1) + 1] < cr[(iii << 1) + 1])
	  {
	    iii = ic;
	  }
     }
  aa1 = MAXCOOR  / (xmax - xmin);
  aa2 = MAXCOOR / (ymax - ymin);
  aa1 = amin (aa1, aa2);
  aa2 = aa1 * (cr[(iii << 1) + 2] - ymin);
  trfri[1] = aa1;
  trfri[2] = cr[(iii << 1) + 1];
  trfri[3] = ymin;
  trfri[4] = aa2;
  i_1 = *nbs;
  for (ic = 1; ic <= i_1; ++ic)
     {
       r_1 = aa1 * (cr[(ic << 1) + 1] - cr[(iii << 1) + 1]);
       c[(ic << 1) + 1] = (long) i_nint (r_1);
       r_1 = aa1 * (cr[(ic << 1) + 2] - ymin) - aa2;
       c[(ic << 1) + 2] = (long) i_nint (r_1);
/* Computing 2nd power */
       i_2 = c[(ic << 1) + 1];
/* Computing 2nd power */
       i_3 = c[(ic << 1) + 2];
       nu[ic] = (LONG8) i_2 * (LONG8) i_2 + (LONG8) i_3 * (LONG8) i_3;
     }
/* ---------------------------------------------------------- */
  mshtr1_ (&nu[1], &tri[1], nbs);
  ip = 1;
  xx = nu[ip];
  i_1 = *nbs;
  for (jc = 1; jc <= i_1; ++jc)
     {
       if (nu[jc] > xx)
	  {
	    i_2 = jc - ip;
	    mshtr1_ (&nu[ip], &tri[ip], &i_2);
	    i_2 = jc - 2;
	    for (i = ip; i <= i_2; ++i)
	       {
		 if (nu[i] == nu[i + 1])
		    {
		      ++ierr;
	             if (ierr <10) 
	             printf(" The points %ld and %ld are too close \n",tri[i],tri[i+1]);
		    }
	       }
	    xx = nu[jc];
	    ip = jc;
	  }
       ic = tri[jc];
       nu[jc] = c[(ic << 1) + 2];
     }
  i_1 = *nbs - ip;
  mshtr1_ (&nu[ip], &tri[ip], &i_1);
  i_1 = jc - 2;
  for (i = ip; i <= i_1; ++i)
     {
       if (nu[i] == nu[i + 1])
	  {
	    if (ierr <10) 
	      printf(" The points %ld and %ld are to close \n",tri[i],tri[i+1]);
	    ++ierr;
	  }
     }
  if (ierr != 0)
     {
       *err = 2;
       printf("mshptg bug 2 \n");

       return 0;
     }
  k = 2;
L50:
  if (k <= *nbs)
     {
       ++k;
       det = (LONG8) c[(tri[2] << 1) + 1] * (LONG8) c[(tri[k] << 1) + 2] 
           - (LONG8) c[(tri[2] << 1) + 2] * (LONG8) c[(tri[k] << 1) + 1];
       if (det == 0)
	  {
	    goto L50;
	  }
     }
  else
     {
       *err = 3;
       return 0;
     }
/*     k est le premier point non aligne */
  trik = tri[k];
  for (j = k - 1; j >= 3; --j)
     {
       tri[j + 1] = tri[j];
     }
  tri[3] = trik;
  if (det < 0)
     {
/*       on inverse les  points 2 3 tries */
       tri3 = tri[3];
       tri[3] = tri[2];
       tri[2] = tri3;
     }
  return 1;
}				/* mshtri_ */

/* ********************************************************************** */
int 
mshtr1_ (LONG8 *criter, long *record, long *n)
{
  /* System generated locals */
  long            i_1;

  /* Local variables */
  long     i, j, l, r, rec;
  LONG8 crit;
/*     trie selon les valeurs de criter croissantes */
/*     record suit le reordonnancement */


  /* Parameter adjustments */
  --record;
  --criter;

  /* Function Body */
  if (*n <= 1)
     {
       return 0;
     }
  l = *n / 2 + 1;
  r = *n;
L2:
  if (l <= 1)
     {
       goto L20;
     }
  --l;
  rec = record[l];
  crit = criter[l];
  goto L3;
L20:
  rec = record[r];
  crit = criter[r];
  record[r] = record[1];
  criter[r] = criter[1];
  --r;
  if (r == 1)
     {
       goto L999;
     }
L3:
  j = l;
L4:
  i = j;
  j <<= 1;
  if ((i_1 = j - r) < 0)
     {
       goto L5;
     }
  else if (i_1 == 0)
     {
       goto L6;
     }
  else
     {
       goto L8;
     }
L5:
  if (criter[j] < criter[j + 1])
     {
       ++j;
     }
L6:
  if (crit >= criter[j])
     {
       goto L8;
     }
  record[i] = record[j];
  criter[i] = criter[j];
  goto L4;
L8:
  record[i] = rec;
  criter[i] = crit;
  goto L2;
L999:
  record[1] = rec;
  criter[1] = crit;
  return 0;
}				/* mshtr1_ */