mshptg.cpp

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

/* ********************************************************************** */
int 
mshcvx_ (long direct, long *c, long *nu, long *pfold, long *err)
{
  static long     t, a4, a5, i1, i2, i3, i4, i5, i6, s1, s2, s3, t4, t5,
                  pf, pp, ps;

//    extern int mshopt_();
  static long     tt4, tt5,  ppf, psf;
  LONG8 det;

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

  /* Function Body */
  if (direct)
     {
       pp = 3;
       ps = 4;
       i1 = 1;
       i2 = 3;
       i3 = 2;
       i4 = 6;
       i5 = 5;
       i6 = 4;
     }
  else
     {
       pp = 4;
       ps = 3;
       i1 = 1;
       i2 = 2;
       i3 = 3;
       i4 = 4;
       i5 = 5;
       i6 = 6;
     }
L10:
  ppf = *pfold;
  pf = nu[ps + *pfold * 6];
  psf = nu[ps + pf * 6];
  s1 = nu[ppf * 6 + 1];
  s2 = nu[pf * 6 + 1];
  s3 = nu[psf * 6 + 1];
  det =    (LONG8) (c[(s2 << 1) + 1] - c[(s1 << 1) + 1]) 
         * (LONG8) (c[(s3 << 1) + 2] - c[(s1 << 1) + 2]) 
         - (LONG8) (c[(s2 << 1) + 2] - c[(s1 << 1) + 2]) 
         * (LONG8) (c[(s3 << 1) + 1] - c[(s1 << 1) + 1]);
         
  if (!(direct) && det > 0 || direct && det < 0)
     {
/*       on ajoute un triangle t et on detruit une arete */
/*       ----------------------------------------------- */
       if (direct)
	  {
	    tt4 = nu[ppf * 6 + 2];
	    tt5 = nu[pf * 6 + 2];
	  }
       else
	  {
	    tt4 = nu[pf * 6 + 2];
	    tt5 = nu[psf * 6 + 2];
	  }
       t4 = tt4 / 8;
       t5 = tt5 / 8;
       a4 = tt4 - (t4 << 3);
       a5 = tt5 - (t5 << 3);
/*       destruction de l'arete frontiere en pf */
/*       -------------------------------------- */
       nu[ps + ppf * 6] = psf;
       nu[pp + psf * 6] = ppf;
/*       on remplace l'arete frontiere par l'element genere */
/*       --------------------------------------------------- */
       t = pf;
/*       update de l'arete non detruite */
/*       ------------------------------ */
       if (direct)
	  {
	    nu[ppf * 6 + 2] = (t << 3) + i6;
	  }
       else
	  {
	    nu[psf * 6 + 2] = (t << 3) + i6;
	  }
/*       on cree l'element */
/*       ----------------- */
       nu[i1 + t * 6] = s1;
       nu[i2 + t * 6] = s2;
       nu[i3 + t * 6] = s3;
       nu[i4 + t * 6] = (t4 << 3) + a4;
       nu[i5 + t * 6] = (t5 << 3) + a5;
       if (direct)
	  {
	    nu[i6 + t * 6] = -ppf;
	  }
       else
	  {
	    nu[i6 + t * 6] = -psf;
	  }
       nu[a4 + t4 * 6] = (t << 3) + i4;
       nu[a5 + t5 * 6] = (t << 3) + i5;
       mshopt_ (&c[3], &nu[7], &t5, a5, err);
       if (*err != 0)
	  {
	    return 0;
	  }
       goto L10;
     }
  return 1;
}				/* mshcvx_ */

/* ********************************************************************** */
int 
mshcxi_ (long *c, long *nu, long *tri, long *nbs, long *tete, long *err)
{
  /* System generated locals */
  long            i_1;

  /* Local variables */
  static long     sfree, ttaf, i, j, s, t, pf;

//    extern long mshlcl_();
  //    extern int mshcvx_(), mshopt_();
  static long     iaf, taf, npf, ppf, psf;

/*     initialisation de la sfree liste dans nu */
  /* Parameter adjustments */
  --tri;
  nu -= 7;
  c -= 3;

  /* Function Body */
  i_1 = *nbs + *nbs - 2;
  for (i = 1; i <= i_1; ++i)
     {
       nu[i * 6 + 1] = i + 1;
       for (j = 2; j <= 6; ++j)
	  {
	    nu[j + i * 6] = 0;
	  }
     }
  nu[(*nbs + *nbs - 2) * 6 + 1] = 0;
  sfree = 1;
/*     initialisation du premier triangle */
  t = sfree;
  sfree = nu[sfree * 6 + 1];
/*     initialisation de la liste frontiere */
  *tete = sfree;
  pf = sfree;
  for (i = 1; i <= 3; ++i)
     {
       nu[i + t * 6] = tri[i];
       nu[i + 3 + t * 6] = -pf;
       ppf = pf;
       sfree = nu[pf * 6 + 1];
       pf = sfree;
       if (i == 3)
	  {
	    pf = *tete;
	  }
       nu[ppf * 6 + 1] = tri[i];
       nu[ppf * 6 + 2] = i + 3 + (t << 3);
       nu[ppf * 6 + 4] = pf;
       nu[pf * 6 + 3] = ppf;
     }
  i_1 = *nbs;
  for (i = 4; i <= i_1; ++i)
     {
       s = tri[i];
       pf = mshlcl_ (&c[3], &nu[7], tete, &s);
/*      creation d'un nouveau triangle et modification de la
   frontiere */
/*      
   --------------------------------------------------------------
 */
       t = sfree;
       sfree = nu[sfree * 6 + 1];
       npf = sfree;
       sfree = nu[sfree * 6 + 1];
       ppf = nu[pf * 6 + 3];
       psf = nu[pf * 6 + 4];
       ttaf = nu[pf * 6 + 2];
       taf = ttaf / 8;
       iaf = ttaf - (taf << 3);

/*                  npf */
/*               1  x s               --- */
/*                 / \                --- */
/*              4 /   \ 6        ---  vide --- */
/*               /  t  \              --- */
/*            2 /   5   \ 3           --- */
/* ------ --<---x---------x---------x- frontiere--<--- */
/*          psf \  iaf  /  pf         --- */
/*               \ taf /         --- omega --- */
/*                \   /               --- */
/*                 \ /                --- */
/*                  x                 --- */
/*                                    --- */
/*     generation  de l'element t */
       nu[t * 6 + 1] = s;
       nu[t * 6 + 2] = nu[psf * 6 + 1];
       nu[t * 6 + 3] = nu[pf * 6 + 1];
       nu[t * 6 + 4] = -npf;
       nu[t * 6 + 5] = (taf << 3) + iaf;
       nu[t * 6 + 6] = -pf;
       nu[iaf + taf * 6] = (t << 3) + 5;
/*      update de la liste frontiere */
       nu[npf * 6 + 4] = psf;
       nu[pf * 6 + 4] = npf;
       nu[npf * 6 + 3] = pf;
       nu[psf * 6 + 3] = npf;
       nu[npf * 6 + 1] = s;
       nu[npf * 6 + 2] = (t << 3) + 4;
       nu[pf * 6 + 2] = (t << 3) + 6;
       mshopt_ (&c[3], &nu[7], &t, 5L, err);
       if (*err != 0)
	  {
	    return 0;
	  }
       mshcvx_ (1, &c[3], &nu[7], &npf, err);
       if (*err != 0)
	  {
	    return 0;
	  }
       mshcvx_ (0, &c[3], &nu[7], &npf, err);
       if (*err != 0)
	  {
	    return 0;
	  }
     }
  return 1;
}				/* mshcxi_ */

/* ********************************************************************** */
int 
mshopt_ (long *c, long *nu, long *t, long a, long *err)
{
  /* Initialized data */

  static long     mod3[3] =
  {2, 3, 1};

  /* System generated locals */
  long            i_1;
  double          d_1;

  /* Local variables */
  static long     pile[4096] /* was [2][256] */ ;
  static double    reel1, reel2;
  static double   reel8;
  static long     i, a1, a2, s1, t1, t2, s2, s3, s4, aa, i11, i12, i13,
                  i21, i22, i23, tt;
  static long     tt1;
  LONG8  cos1, cos2, sin1, sin2, sgn;
  int kerr21=0;
  /* Parameter adjustments */
  nu -= 7;
  c -= 3;

  /* Function Body */
  i = 1;
  pile[(i << 1) - 2] = *t;
  pile[(i << 1) - 1] = a;
L10:
  if (i > 0)
     {
       t1 = pile[(i << 1) - 2];
       a1 = pile[(i << 1) - 1];
       --i;
       if (t1 <= 0)
	  {
	    goto L10;
	  }
       tt1 = nu[a1 + t1 * 6];
       if (tt1 <= 0)
	  {
	    goto L10;
	  }
       t2 = tt1 / 8;
       a2 = tt1 - (t2 << 3);
       i11 = a1 - 3;
       i12 = mod3[i11 - 1];
       i13 = mod3[i12 - 1];
       i21 = a2 - 3;
       i22 = mod3[i21 - 1];
       i23 = mod3[i22 - 1];
       s1 = nu[i13 + t1 * 6];
       s2 = nu[i11 + t1 * 6];
       s3 = nu[i12 + t1 * 6];
       s4 = nu[i23 + t2 * 6];
       sin1 =     (LONG8) (c[(s3 << 1) + 2] - c[(s1 << 1) + 2]) 
                * (LONG8) (c[(s2 << 1) + 1] - c[(s1 << 1) + 1]) 
              -   (LONG8) (c[(s3 << 1) + 1] - c[(s1 << 1) + 1]) 
                * (LONG8) (c[(s2 << 1) + 2] - c[(s1 << 1) + 2]);
                
       cos1 =     (LONG8) (c[(s3 << 1) + 1] - c[(s1 << 1) + 1]) 
                * (LONG8) (c[(s3 << 1) + 1] - c[(s2 << 1) + 1]) 
              +   (LONG8) (c[(s3 << 1) + 2] - c[(s1 << 1) + 2]) 
                * (LONG8) (c[(s3 << 1) + 2] - c[(s2 << 1) + 2]);
       if (sin1 == 0 && cos1 == 0)
	  {
	    *err = 20;
	    return 0;
	  }
/*       b est la cotangente de angle (s1,s3,s2) */

       sin2 =   (LONG8) (c[(s4 << 1) + 1] - c[(s1 << 1) + 1]) 
              * (LONG8) (c[(s2 << 1) + 2] - c[(s1 << 1) + 2]) 
              - (LONG8) (c[(s4 << 1) + 2] - c[(s1 << 1) + 2]) 
              * (LONG8) (c[(s2 << 1) + 1] - c[(s1 << 1) + 1]);
     // FH correct 6/11/2005 forgotted cast          
       cos2 =   (LONG8) (c[(s4 << 1) + 1] - c[(s2 << 1) + 1]) 
              * (LONG8) (c[(s4 << 1) + 1] - c[(s1 << 1) + 1]) 
              + (LONG8) (c[(s4 << 1) + 2] - c[(s2 << 1) + 2]) 
              * (LONG8) (c[(s4 << 1) + 2] - c[(s1 << 1) + 2]);
              
       reel1 = (double) cos2 *(double) sin1;
       reel2 = (double) cos1 *(double) sin2;

       if (aabs (reel1) + aabs (reel2) >= (double ) DLONG8LONG8)
	  {
	    reel8 = (double) cos2 *(double) sin1 + (double) cos1 *(double) sin2;

/* Computing MIN */
	    d_1 = amax (reel8, -1.);
	    reel8 = amin (d_1, 1.);
	    sgn = (LONG8) reel8;
	  }
       else
	  {
	    sgn = cos2 * sin1 + cos1 * sin2;
	  }
/* Computing MIN */
       i_1 = amin(amax (sgn, -1),1);
      // cout << sgn << " " << i_1 << endl;
       if ( i_1 * sin1 >= 0)
	  {
	    goto L10;
	  }
/*       on inverse le quadrilatere */
/*       update des sommets */
/* ------------------------- */
       nu[i12 + t1 * 6] = s4;
       nu[i22 + t2 * 6] = s1;
/*       update des aretes a1,a2 */
/* ------------------------------- */
       tt1 = nu[i22 + 3 + t2 * 6];
       nu[a1 + t1 * 6] = tt1;
       if (tt1 > 0)
	  {
	    tt = tt1 / 8;
	    aa = tt1 - (tt << 3);
	    nu[aa + tt * 6] = a1 + (t1 << 3);
	  }
       else if (tt1 != nothing)
	  {
	    nu[-tt1 * 6 + 2] = a1 + (t1 << 3);
	  }
       tt1 = nu[i12 + 3 + t1 * 6];
       nu[a2 + t2 * 6] = tt1;
       if (tt1 > 0)
	  {
	    tt = tt1 / 8;
	    aa = tt1 - (tt << 3);
	    nu[aa + tt * 6] = a2 + (t2 << 3);
	  }
       else if (tt1 != nothing)
	  {
	    nu[-tt1 * 6 + 2] = a2 + (t2 << 3);
	  }
       nu[i12 + 3 + t1 * 6] = i22 + 3 + (t2 << 3);
       nu[i22 + 3 + t2 * 6] = i12 + 3 + (t1 << 3);
       if (i + 4 > 1024)
	  {
	    if(kerr21)
	      cout << " Bizarre mshptg err 21 pile too small (continue)  "<< endl;
	    if(kerr21++<10000)
		    goto L10; 
	    *err = 21;
	    return 0;
	  }
       ++i;
       pile[(i << 1) - 2] = t1;
       pile[(i << 1) - 1] = a1;
       ++i;
       pile[(i << 1) - 2] = t2;
       pile[(i << 1) - 1] = a2;
       ++i;
       pile[(i << 1) - 2] = t1;
       pile[(i << 1) - 1] = i13 + 3;
       ++i;
       pile[(i << 1) - 2] = t2;
       pile[(i << 1) - 1] = i23 + 3;
       goto L10;
     }
  return 1;
}				/* mshopt_ */

/* ********************************************************************** */
int 
mshfrt_ (long *c, long *nu, long *nbs, long *arete, long nba, long *sd,
	 long nbsd, long *reft, long *w, long *err)
{
  /* Initialized data */

  static long     p3[3] =
  {2, 3, 1};

  /* System generated locals */
  long            i_1, i_2, i_3;

  /* Local variables */
  static long     nbac, ifrt, a, i, t, itera, s1, s2;

//    extern int mshfr1_();
  static long     ie, ap, ta, is, nbacpp;
  static long     is1, ss1, s2t, s3t, isd, jsd, nbt, err1;
  LONG8 det2, det3;

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

  /* Function Body */
  if (nba == 0)
     {
       return 0;
     }
  ifrt = 0;
  nbt = *nbs + *nbs - 2;
  i_1 = *nbs;
  for (i = 1; i <= i_1; ++i)
     {
       reft[i] = 0;
     }
  i_1 = nba;
  for (i = 1; i <= i_1; ++i)
     {
       reft[arete[(i << 1) + 1]] = nothing;
       reft[arete[(i << 1) + 2]] = nothing;
     }
  nbac = 0;
  i_1 = nba;
  for (a = 1; a <= i_1; ++a)
     {
/* Computing MIN */
       i_2 = arete[(a << 1) + 1], i_3 = arete[(a << 1) + 2];
       s1 = amin (i_2, i_3);
/* Computing MAX */
       i_2 = arete[(a << 1) + 1], i_3 = arete[(a << 1) + 2];
       s2 = amax (i_2, i_3);
       if (s1 == s2)
	  {
	    ++nbac;
	  }
       else
	  {
	    i = reft[s1];
	  L25:
	    if (i != nothing)
	       {
/* Computing MAX */
		 i_2 = arete[(i << 1) + 1], i_3 = arete[(i << 1) + 2];
		 if (s2 == amax (i_2, i_3))
		    {
		      ++nbac;
		    }
		 else
		    {
		      i = w[i];
		      goto L25;
		    }
	       }
	    else
	       {
		 w[a] = reft[s1];
		 reft[s1] = a;
	       }
	  }
     }
  nbacpp = 1;
  itera = 0;
  err1 = 0;
L50:
  ++itera;
  if (err1 != 0)
     {
       *err = err1;
       return 0;
     }
  if (nbac < nba)
     {
       if (nbacpp == 0)
	  {
	    i_1 = *nbs;
	    for (i = 1; i <= i_1; ++i)
	       {
		 a = reft[i];
	       L60:
		 if (a > 0)
		    {
		      s1 = arete[(i << 1) + 1];
		      s2 = arete[(i << 1) + 2];
		      a = w[a];
		      goto L60;
		    }
	       }
	    *err = 7;
	    return 0;
	  }
/* --------------------------------------------------------------------- */
/*     on s'occupe des aretes a forcer */
/* --------------------------------------------------------------------- */
       nbacpp = 0;
       i_1 = nbt;
       for (ie = 1; ie <= i_1; ++ie)
	  {
	    if (nu[ie * 6 + 5] != 0)
	       {
		 for (is = 1; is <= 3; ++is)
		    {
		      s1 = nu[is + ie * 6];
		      s2t = nu[p3[is - 1] + ie * 6];
		      ss1 = amin (s1, s2t);
		      ap = 0;
		      a = reft[ss1];
		    L80:
		      if (a > 0)
			 {
/* Computing MAX */
			   i_2 = arete[(a << 1) + 1], i_3 = arete[(a << 1) + 2];
			   s2 = amax (i_2, i_3);
			   t = ie;
			   ta = 0;
			   if (s2 == amax (s1, s2t))
			      {
				if (nu[is + 3 + ie * 6] > 0)
				   {
				     ta = nu[is + 3 + ie * 6] / 8;
				     i = nu[is + 3 + ie * 6] - (ta << 3);
				     nu[i + ta * 6] = nothing;
				   }
				nu[is + 3 + ie * 6] = nothing;
				goto L100;
			      }
			   ap = a;
			   a = w[a];
			   goto L80;
			 }
		      if (itera == 1)
			 {
			   goto L110;
			 }
		      ss1 = s1;