gibbs.cpp

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

	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 Mesh::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 = nv , 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 = number(triangles[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 = number(triangles[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 = number(triangles[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 */
     return 0;} /* gibbsv_ */

int Mesh::renum()
/* -------- 
	renumber vertices by gibbs method; updates triangle and edge array
	in:   mesh  
	out:   mesh
 	auxiliary arrays: ptvois,vois,r,m,nv,nx,ny,nn,w1,w2,f 
 	all of size nv+1 except vois (10(nv+1)) and nv (2(nv+1))
 	err = -1 : memory alloc pb; err = -3: fatal erreur  gibbs 2 : pb racine
*/
{
    long   pfold, pfnew;
    long* ptvois=NULL;
    long* vois=NULL;
    long* nn =NULL;
    long* r =NULL;
    long* m =NULL;
    long* nnv =NULL;
    long* nx =NULL;
    long* ny =NULL;
    long* w1 =NULL;
    long* w2=NULL;
    long nbvoisin =  10*nv;
    long printint=0, iodev=6;
    int err=0;
  	ptvois = new long[nv+1]; 		//(long*)calloc((long)(nv + 1) , sizeof(long));
	nn = 	 new long[3*nt]; 			//(long*)calloc(3 * nt ,sizeof(long));
	vois = 	 new long[nbvoisin+10];	//(long*)calloc((long)(nbvoisin + 10) , sizeof(long)); 
	r = 	 new long[nv+1];				//(long*)calloc((long)(nv + 1) , sizeof(long));
	if((!ptvois)||(!nn)||(!vois)||(!r)) return -1;
	err = gibbsv(ptvois,vois,&nbvoisin,r,nn) ;
	delete [] nn;					// free(nn);
	if(err==0)
	{
       m = new long[nv+1];
       nn = new long[nv+1];
       nnv = new long[(nv+1)<<1];
       nx = new long[nv+1];
       ny = new long[nv+1];
       w1 = new long[nv+1];
       w2 = new long[nv+1];
	   long lnv = nv;
       err = gibbsa_ (&lnv, ptvois, vois, r, m, nnv, nx, ny, nn, w1, w2, &pfold, &pfnew,
		      &printint, &iodev);
       delete [] m;
       delete [] nnv;
       delete [] nn;
       delete [] nx;
       delete [] ny;
       delete [] w1;
       delete [] w2;
     }

  delete [] vois;
  delete [] ptvois;
     if(verbosity>1)
  cout << " -- Mesh: Gibbs: old skyline = " << pfold << "  new skyline = " << pfnew << endl;
  if (err == 0 && (pfnew <= pfold))
     {
       int i,j;
       for ( i=0;i<nv;i++)
	 r[i] -= 1; 
       Vertex  * f= new Vertex[nv];
       for (i = 0; i < nv; ++i)
	 f[i] = vertices[i];
       for (i = 0; i < nv; ++i)
	 vertices[r[i]]  = f[i];
       delete [] f;       
       
       for (j = 0; j < nt; ++j)  // updates triangle array
	  triangles[j].Renum(vertices,r);
       
       for (j = 0; j < neb; ++j)	// updates edge array
	    bedges[j].Renum(vertices,r);
	  
	  if (BoundaryAdjacencesHead )
	   {
          for (int i=0;i<nv;i++) 
            BoundaryAdjacencesHead[i]=-1;
          int j2=0;
          for (int j=0;j<neb;j++)
           for (int k=0;k<2;k++,j2++)
             {  
              int v = number(bedges[j][k]);
              ffassert(v >=0 && v < nv);
              BoundaryAdjacencesLink[j2]=BoundaryAdjacencesHead[v];
              BoundaryAdjacencesHead[v]=j2;
            }
        }
      for (int it=0;it<nt;it++)
       for (int j=0;j<3;j++)
         TriangleConteningVertex[(*this)(it,j)]=it;
        
       if (quadtree) {
         delete quadtree;quadtree=0;
         MakeQuadTree(); }

      
     }
  delete [] r;
  return err;
} 
int FESpace::renum()
/* -------- 
	renumber vertices by gibbs method; updates triangle and edge array
	in:   mesh  
	out:   mesh
 	auxiliary arrays: ptvois,vois,r,m,nv,nx,ny,nn,w1,w2,f 
 	all of size nv+1 except vois (10(nv+1)) and nv (2(nv+1))
 	err = -1 : memory alloc pb; err = -3: fatal erreur  gibbs 2 : pb racine
*/
{
  if (cdef==0) return -2;
  if (cdef->NodesOfElement==0) return -2;
    int nv = NbOfNodes;
    int nt = NbOfElements;// ;Th.nt;
    long  pfold =0, pfnew=0;
    long* ptvois=NULL;
    long* vois=NULL;
    long* nn =NULL;
    long* r =NULL;
    long* m =NULL;
    long* nnv =NULL;
    long* nx =NULL;
    long* ny =NULL;
    long* w1 =NULL;
    long* w2=NULL;
    long nbvperelem = 20; // pour le P1
  //  cout << "gibbs: nbvperelem =" << nbvperelem << endl;
    long nbvoisin =  (nbvperelem)*nv; 
    long printint=0, iodev=6;
    int err=0;
    int nnx= SizeToStoreAllNodeofElement();
  	ptvois = new long[nv+1]; 		
	nn = 	 new long[nnx]; 		
	vois = 	 new long[nbvoisin+10];	
	r = 	 new long[nv+1];			
	if((!ptvois)||(!nn)||(!vois)||(!r)) return -1;
	err = gibbsv(ptvois,vois,&nbvoisin,r,nn) ;
	delete [] nn;					
	if(err==0)
	{
       m = new long[nv+1];
       nn = new long[nv+1];
       nnv = new long[(nv+1)<<1];
       nx = new long[nv+1];
       ny = new long[nv+1];
       w1 = new long[nv+1];
       w2 = new long[nv+1];
	   long lnv = nv;
       err = gibbsa_ (&lnv, ptvois, vois, r, m, nnv, nx, ny, nn, w1, w2, &pfold, &pfnew,
		      &printint, &iodev);
       delete [] m;
       delete [] nnv;
       delete [] nn;
       delete [] nx;
       delete [] ny;
       delete [] w1;
       delete [] w2;
     }
    else 
      cerr << " Not enought memory (bug)" <<  nbvoisin << " " << nbvoisin/NbOfElements <<  endl;
  delete [] vois;
  delete [] ptvois;
  if (err==0 && verbosity>1) 
  cout << " FESpace:Gibbs: old skyline = " << pfold << "  new skyline = " << pfnew << endl;
  if (err == 0 && (pfnew <= pfold))
     {
       int i;
      // cout << *this << endl;
       for ( i=0;i<nv;i++)
	     r[i] -= 1; 
	    cdef->renum(r,nnx);  
        }
     
  delete [] r;
  return err;
} 
/* function */ int FESpace::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 = NbOfNodes , nbt = NbOfElements;
     integer jj,kk;
/*--- 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 < NbOfNodesInElement(i) ; ++j) {
	     ss =  (*this)(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,jj= NbOfNodesInElement(i); j <jj  ; ++j) {
	     ss = (*this)(i,j) +1;//number(triangles[i][j])+1;
	     if ( ! ( ss>0 && ss <= nbss)) 
	       {
	          cout << "bug " << ss << " " <<  i << " " << j << endl;
	          exit(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,kk=NbOfNodesInElement(T); k < kk; ++k) {
	   ss = (*this)(T,k) +1;//number(triangles[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 */
     return 0;} /* gibbsv_ */						
/*  message d'erreur:         *err = 2;    print*,'pas assez de place memoire'   */
			
/*  message d'erreur:         *err = 2;    print*,'pas assez de place memoire'   */