mshptg.cpp

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

		      ap = 0;
		      a = reft[ss1];
		    L90:
		      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;
/*             recherche si l' element coupe l''arete a */
			   is1 = is;
			   s3t = nu[p3[p3[is - 1] - 1] + t * 6];
			   det2 = (LONG8) (c[(s2t << 1) + 1] - c[(s1 << 1) + 1]) * (LONG8) (c[( s2 << 1) + 2] - c[(s1 << 1) + 2]) 
			        - (LONG8) (c[(s2t << 1) + 2] - c[(s1 << 1) + 2]) * (LONG8) (c[( s2 << 1) + 1] - c[(s1 << 1) + 1]);
			   det3 = (LONG8) (c[(s3t << 1) + 1] - c[(s1 << 1) + 1]) * (LONG8) (c[( s2 << 1) + 2] - c[(s1 << 1) + 2]) 
			        - (LONG8) (c[(s3t << 1) + 2] - c[(s1 << 1) + 2]) * (LONG8) (c[( s2 << 1) + 1] - c[(s1 << 1) + 1]);
			   if (det2 > 0 && det3 < 0)
			      {
				mshfr1_ (&c[3], &nu[7], &t, &ta, &is1, &s2, err);
				if (*err != 0)
				   {
				     return 0;
				   }
				goto L100;
			      }
			   else if (det2 == 0 && det3 < 0 && reft[s2t] == 0)
			      {
				err1 = 10;
				printf(" det = %ld %ld %ld %ld %ld == %ld \n ",
				   (c[(s2t << 1) + 1] - c[(s1 << 1) + 1]), (c[( s2 << 1) + 2] - c[(s1 << 1) + 2]) ,
				   (c[(s2t << 1) + 2] - c[(s1 << 1) + 2]), (c[( s2 << 1) + 1] - c[(s1 << 1) + 1]) ,
				   (c[(s2t << 1) + 1] - c[(s1 << 1) + 1]) * (c[( s2 << 1) + 2] - c[(s1 << 1) + 2]) ,
				   (c[(s2t << 1) + 2] - c[(s1 << 1) + 2]) * (c[( s2 << 1) + 1] - c[(s1 << 1) + 1])
				   );
				   
				 
				printf("bug 2, mshptg: point %ld is on boundary edge %ld %ld  \n",s2t,i_2,i_3);
			      }
			   else if (det2 > 0 && det3 == 0 && reft[s3t] == 0)
			      {
				err1 = 10;
				printf(" det = %ld %ld %ld %ld  %ld %ld \n ",
				    (c[(s3t << 1) + 1] - c[(s1 << 1) + 1]),  (c[( s2 << 1) + 2] - c[(s1 << 1) + 2]) ,
				    (c[(s3t << 1) + 2] - c[(s1 << 1) + 2]),  (c[( s2 << 1) + 1] - c[(s1 << 1) + 1]) ,
				    (c[(s3t << 1) + 1] - c[(s1 << 1) + 1]) * (c[( s2 << 1) + 2] - c[(s1 << 1) + 2]) ,
				    (c[(s3t << 1) + 2] - c[(s1 << 1) + 2]) *  (c[( s2 << 1) + 1] - c[(s1 << 1) + 1])
				    );
				printf("bug 2, mshptg: point %ld is on  boundary %ld %ld\n",s3t,i_2,i_3);
			      }
			   ap = a;
			   a = w[a];
			   goto L90;
			 }
		      goto L110;
		    L100:
		      ++nbacpp;
		      if (ap == 0)
			 {
			   reft[ss1] = w[a];
			 }
		      else
			 {
			   w[ap] = w[a];
			 }
		      if (nbac + nbacpp == nba)
			 {
			   goto L130;
			 }
		    L110:
		      ;
		    }
	       }
	  }
       nbac += nbacpp;
       goto L50;
     }
L130:
/* ----------------------------------------------------------------------- */
/*     prise en compte des sous domaines */
/* ----------------------------------------------------------------------- */
/*  add FH   si pas de nbsd ---  jan 2004 */
    
    if (nbsd<0) {
     i_1 = nbt;
    for (ie = 1; ie <= i_1; ++ie)
     { 
       if (nu[ie * 6 + 1] > 0)
	     {
	    nu[ie * 6 + 1] = -nu[ie * 6 + 1];
	     }
	 reft[ie]=1;
     }
     
      goto L205; //  pas triangle retire
      }    
     if (nbsd == 0) {
	 long i__1 = nbt,nbssd,exter,i__,headt,nst,j;
	for (t = 1; t <= i__1; ++t) {
	    reft[t] = -1073741824;
	}
	nbssd = 0;
/*         if(impre.gt.1) print *,'nbsd .eq. 0 =>  recherche de ssd' */
	i__1 = nbt;
	for (t = 1; t <= i__1; ++t) {
	    if (nu[t * 6 + 1] > 0 && nu[t * 6 + 5] != 0) {
		++nbssd;
		exter = 0;
		i__ = 2;
		w[i__ - 1] = t;
		w[i__] = 3;
		reft[t] = 0;
		headt = t;
		nu[t * 6 + 1] = -nu[t * 6 + 1];
		nst = 1;
/*              print *,' ssd ',nbssd */
L131:
		if (i__ > 0) {
		    ++w[i__];
		    if (w[i__] <= 6) {
			ta = nu[w[i__] + w[i__ - 1] * 6];
			if (ta <= 0) {
			    if (ta != -1073741824) {
				exter = 1;
			    }
			} else if (ta > 0) {
			    ta /= 8;
			    if (nu[ta * 6 + 1] > 0) {
/*                           print *,ta */
				nu[ta * 6 + 1] = -nu[ta * 6 + 1];
				++nst;
				reft[ta] = headt;
				headt = ta;
				w[i__ + 1] = ta;
				w[i__ + 2] = 3;
				i__ += 2;
			    }
			}
		    } else {
			i__ += -2;
		    }
		    goto L131;
		}

		if (exter) {
		    --nbssd;
		    i__ = headt;
L133:
		    if (i__ > 0) {
			j = reft[i__];
/*                     print *,i */
			reft[i__] = 0;
			nu[i__ * 6 + 1] = 0;
			i__ = j;
			goto L133;
		    }
		} else {
		    i__ = headt;
L136:
		    if (i__ > 0) {
			j = reft[i__];
			reft[i__] = nbssd;
			i__ = j;
			goto L136;
		    }
		}
	    }
	}
	goto L205;
    }
/*  fin ajoute FH jan 2004  */


  i_1 = *nbs + nbsd + nbsd;
  for (i = 1; i <= i_1; ++i)
     {
       w[i] = 0;
     }
  i_1 = nbsd;
  for (i = 1; i <= i_1; ++i)
     {
       a = (i_2 = sd[(i << 1) + 1], aabs (i_2));
/* Computing MIN */
       i_2 = arete[(a << 1) + 1], i_3 = arete[(a << 1) + 2];
       s1 = amin (i_2, i_3);
       w[i + i] = w[s1 + nbsd + nbsd];
       w[s1 + nbsd + nbsd] = i;
     }
  i_1 = nbt;
  for (t = 1; t <= i_1; ++t)
     {
       reft[t] = nothing;
       if (nu[t * 6 + 6] != 0)
	  {
	    for (i = 1; i <= 3; ++i)
	       {
/* Computing MIN */
		 i_2 = nu[i + t * 6], i_3 = nu[p3[i - 1] + t * 6];
		 ss1 = amin (i_2, i_3);
		 jsd = nbsd + nbsd + ss1;
	       L160:
		 isd = w[jsd];
		 if (isd > 0)
		    {
		      a = sd[(isd << 1) + 1];
		      if (a > 0)
			 {
			   if (nu[i + t * 6] == arete[(a << 1) + 1] && nu[p3[i -
					 1] + t * 6] == arete[(a << 1) + 2])
			      {
				reft[t] = sd[(isd << 1) + 2];
				w[isd + isd - 1] = t;
				w[jsd] = w[isd + isd];
				goto L170;
			      }
			 }
		      else if (a < 0)
			 {
			   if (nu[i + t * 6] == arete[(-a << 1) + 2] && nu[p3[i
				      - 1] + t * 6] == arete[(-a << 1) + 1])
			      {
				reft[t] = sd[(isd << 1) + 2];
				w[isd + isd - 1] = t;
				w[jsd] = w[isd + isd];
				goto L170;
			      }
			 }
		      else
			 {
			    cout << " Err sous domaine " << isd << "ref par  d'arete +/-  = " << a  << " n'est dans aucun triangle " << endl;
			   *err = 11;
			 }
		      jsd = isd + isd;
		      goto L160;
		    }
	       L170:
		 ;
	       }
	  }
     }
  i_1 = nbsd;
  for (isd = 1; isd <= i_1; ++isd)
     {
       if (w[isd + isd - 1] == 0)
	  {
	   cout << " Err sous domaine " << isd << "ref par  d'arete +/-  = " << sd[(isd << 1) + 1]  << " n'est dans aucun triangle " << endl;
	    *err = 11;
	  }
       else
	  {
	    w[isd + isd] = 3;
	  }
     }
  if (*err != 0)
     {
       return 0;
     }
  i = nbsd + nbsd;
L200:
  if (i > 0)
     {
       ++w[i];
       if (w[i] <= 6)
	  {
	    ta = nu[w[i] + w[i - 1] * 6];
	    if (ta > 0)
	       {
		 ta /= 8;
		 if (nu[ta * 6 + 1] > 0)
		    {
		      nu[ta * 6 + 1] = -nu[ta * 6 + 1];
		      if (reft[ta] != reft[w[i - 1]])
			 {
			   if (reft[ta] != nothing)
			      {
			      }
			   else
			      {
				reft[ta] = reft[w[i - 1]];
			      }
			   w[i + 1] = ta;
			   w[i + 2] = 3;
			   i += 2;
			 }
		    }
	       }
	  }
       else
	  {
	    ta = w[i-1] ; 
	    if( nu[ta * 6 + 1]>=0)    nu[ta * 6 + 1]= -nu[ta * 6 + 1]; // pour etre sur que le traingle est bien marque 

	    i += -2;
	  }
       goto L200;
     }
L205:     
  i_1 = nbt;
  for (ie = 1; ie <= i_1; ++ie)
     {
       if (nu[ie * 6 + 1] < 0)
	  {
	    nu[ie * 6 + 1] = -nu[ie * 6 + 1];
	  }
       else
	  {
	    for (i = 1; i <= 6; ++i)
	       {
		 nu[i + ie * 6] = 0;
	       }
	  }
     }
  return 1;
}				/* mshfrt_ */

/* ********************************************************************** */
int 
mshfr1_ (long *c, long *nu, long *it1, long *ita, long *is1, long *s2, long *err)
{
  /* Initialized data */

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

  static long     nbac, t, x, y, l1, l2, l3, s1, s3;

//    extern int mshfr2_();
  static long     la, ta;
  static long     s2t, s3t,  lst[768] /* was [3][256] */ ;
  LONG8 det;

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

  /* Function Body */
  t = *it1;
  s1 = nu[*is1 + t * 6];
  x = c[(*s2 << 1) + 1] - c[(s1 << 1) + 1];
  y = c[(*s2 << 1) + 2] - c[(s1 << 1) + 2];
  nbac = 0;
  l1 = *is1;
  l2 = p3[l1 - 1];
  l3 = p3[l2 - 1];
  s2t = nu[l2 + t * 6];
  s3t = nu[l3 + t * 6];
  la = l2 + 3;
L20:
  ++nbac;
  if (nbac > 256)
     {
       *err = 8;
       return 0;
     }
  lst[nbac * 3 - 2] = t;
  lst[nbac * 3 - 1] = la;
  ta = nu[la + t * 6];
  if (ta <= 0)
     {
       *err = 9;
       return 0;
     }
  t = ta / 8;
  la = ta - (t << 3);
  s3 = nu[p3[la - 3] + t * 6];
  if (s3 != *s2)
     {
       det =   (LONG8) x * (LONG8) (c[(s3 << 1) + 2] - c[(s1 << 1) + 2]) 
             - (LONG8) y * (LONG8) (c[(s3 << 1) + 1] - c[(s1 << 1) + 1]);
       if (det > 0)
	  {
	    la = p3[la - 4] + 3;
	  }
       else if (det < 0)
	  {
	    la = p3[la - 3] + 3;
	  }
       else
	  {
	    printf("mshptg: bug the point %ld  is on boundary \n", s3 ); 
	    *err = 10+s3*10;
	    return 0;
	  }
       goto L20;
     }
  mshfr2_ (&c[3], &nu[7], lst, &nbac, it1, ita, &s1, s2, err);
  return 0;
}				/* mshfr1_ */

/* ********************************************************************** */
int 
mshfr2_ (long *c, long *nu, long *lst, long *nbac, long *t, long *ta,
	 long *ss1, long *ss2, long *err)
{
  /* Initialized data */

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

  /* System generated locals */
  long            i_1;

  /* Local variables */
  static long     i, x, y, a1, a2, pplst, s1, pslst, ptlst, s2, s3, s4,
                  ttlst, t1, t2, aa, i11, i12, i13, i21, i22, i23, x41,
                  y41, tt;

//    extern int mshopt_();
  static long     tt1, aas;
LONG8  det1, det2, det3, det4;

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

  /* Function Body */
  x = c[(*ss1 << 1) + 1] - c[(*ss2 << 1) + 1];
  y = c[(*ss1 << 1) + 2] - c[(*ss2 << 1) + 2];
  i_1 = *nbac - 1;
  for (i = 1; i <= i_1; ++i)
     {
       lst[i * 3 + 1] = i + 1;
     }
  lst[*nbac * 3 + 1] = 0;
  ttlst = 1;
L20:
  ptlst = ttlst;
  pplst = 0;
L30:
  if (ptlst > 0)
     {
       t1 = lst[ptlst * 3 + 2];
       a1 = lst[ptlst * 3 + 3];
       tt1 = nu[a1 + t1 * 6];
       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];
       x41 = c[(s4 << 1) + 1] - c[(s1 << 1) + 1];
       y41 = c[(s4 << 1) + 2] - c[(s1 << 1) + 2];
       det2 = (LONG8) (c[(s2 << 1) + 1] - c[(s1 << 1) + 1]) * (LONG8) y41 
            - (LONG8) (c[(s2 << 1) + 2] - c[(s1 << 1) + 2]) * (LONG8) x41;
            
       det3 = (LONG8) (c[(s3 << 1) + 1] - c[(s1 << 1) + 1]) * (LONG8)  y41 
            - (LONG8) (c[(s3 << 1) + 2] - c[(s1 << 1) + 2]) * (LONG8)  x41;
       if (det2 > 0 && det3 < 0)
	  {
/*         le quadrilataire est convexe on le retourne */
/*         update des sommets */
/* ------------------------- */
	    nu[i12 + t1 * 6] = s4;
	    nu[i22 + t2 * 6] = s1;
/*         update du pointeur suivant */
/* ----------------------------------- */
	    pslst = lst[ptlst * 3 + 1];
	    if (pslst > 0)
	       {
		 aas = lst[pslst * 3 + 3];
		 if (aas == i22 + 3)
		    {
		      lst[pslst * 3 + 2] = t1;
		      lst[pslst * 3 + 3] = i11 + 3;
		    }
	       }
/*         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);
	    
	    det1 = (LONG8) (c[(s1 << 1) + 1] - c[(*ss1 << 1) + 1]) * (LONG8) y 
	         - (LONG8) (c[(s1 << 1) + 2] - c[(*ss1 << 1) + 2]) * (LONG8) x;
	         
	    det4 = (LONG8) (c[(s4 << 1) + 1] - c[(*ss1 << 1) + 1]) * (LONG8) y 
	         - (LONG8) (c[(s4 << 1) + 2] - c[(*ss1 << 1) + 2]) * (LONG8) x;
	    if (det1 < 0 && det4 > 0)
	       {
/*           le sommets s4 est dans omega */
		 lst[ptlst * 3 + 2] = t2;
		 lst[ptlst * 3 + 3] = i22 + 3;
	       }
	    else if (det1 > 0 && det4 < 0)
	       {
/*           le sommets s1 est dans omega */
		 lst[ptlst * 3 + 2] = t1;
		 lst[ptlst * 3 + 3] = i12 + 3;
	       }
	    else
	       {
		 if (pplst == 0)
		    {
		      ttlst = lst[ptlst * 3 + 1];
		      ptlst = ttlst;
		    }
		 else
		    {
		      ptlst = lst[ptlst * 3 + 1];
		      lst[pplst * 3 + 1] = ptlst;
		    }
		 goto L30;
	       }
	  }
       pplst = ptlst;
       ptlst = lst[ptlst * 3 + 1];
       goto L30;
     }
  if (ttlst != 0)
     {
       goto L20;
     }
  nu[i12 + 3 + t1 * 6] = nothing;
  nu[i22 + 3 + t2 * 6] = nothing;
  *t = t2;
  *ta = t1;
  i_1 = *nbac;
  for (i = 1; i <= i_1; ++i)
     {
       mshopt_ (&c[3], &nu[7], &lst[i * 3 + 2], 4L, err);
       mshopt_ (&c[3], &nu[7], &lst[i * 3 + 2], 5L, err);
       mshopt_ (&c[3], &nu[7], &lst[i * 3 + 2], 6L, err);
     }
  return 1;
}				/* mshfr2_ */


}