gquadtree.cpp

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

    bc =b;
    be = b +ll;
    assert(b);
  }
  
  template<class Vertex>    GTree<Vertex>::StorageQuadTreeBox::~StorageQuadTreeBox()
  { //cout <<  "~StorageQuadTreeBox " << this << " n " << n << " b " << b << endl;
    if(n) delete n;
    delete [] b;
  }
  
  template<class Vertex> GTree<Vertex>::~GTree()
  {
    delete sb; 
  }
  
template<class Vertex> ostream& operator <<(ostream& f, const  GTree<Vertex> & qt)
{ 
  f << " the tree "  << endl;
  f << " NbTreeBox = " << qt.NbQuadTreeBox 
    << " Nb Vertices = " <<  qt.NbVertices << endl;
  f << " NbTreeBoxSearch " << qt.NbQuadTreeBoxSearch  
    << " NbVerticesSearch " << qt.NbVerticesSearch << endl;
  f << " SizeOf QuadTree" << qt.SizeOf() << endl;
  //     return  dump(f,*qt.root);
  return  f;
}
  
  template<class Vertex> Vertex *  GTree<Vertex>::NearestVertexWithNormal(const Rd &P)//(long xi,long yj)
  {
    QuadTreeBox * pb[ MaxDeep ];
    int  pi[ MaxDeep  ];
    Zd pp[ MaxDeep];
    int l; // level
    QuadTreeBox * b;
    IntQuad  h=MaxISize,h0;
    IntQuad hb =  MaxISize;
    Zd   p0;
    Zd  plus(RdtoZd(P) );//xi<MaxISize?(xi<0?0:xi):MaxISize-1,yj<MaxISize?(yj<0?0:yj):MaxISize-1);
    
    Vertex *vn=0;
    // init for optimisation ---
    b = root;
    long  n0;
    if (!root->n)
      return vn; // empty tree 
    
    while( (n0 = b->n) < 0) 
      {
	// search the non empty 
	// QuadTreeBox containing  the point (i,j)
	long hb2 = hb >> 1 ;
	int k = plus.Case(hb2);//(iplus,jplus,hb2);// QuadTreeBox number of size hb2 contening i;j
	QuadTreeBox * b0= b->b[k];
	if ( ( b0 == 0) || (b0->n == 0) ) 
	  break; // null box or empty   => break 	    
	NbQuadTreeBoxSearch++;
	b=b0;	
	p0.Add(k,hb2);	
	hb = hb2; 
      }
    
    
    if ( n0 > 0) 
    {  
      for(int k=0;k<n0;k++)
	{
	  Vertex * v=b->v[k];
	  if (v->ninside(P)) {
	    Zd i2 =  VtoZd(v);
	    //   try if is in the right sens -- 
	    h0 = Zd(i2,plus).norm();// h0 = NORM(iplus,i2.x,jplus,i2.y);
	    if (h0 <h) {
	      h = h0;
	      vn = v;}
	    NbVerticesSearch++;}
	}
      if (vn) return vn; 
    }
    // general case -----
    // INITIALISATION OF THE STACK 
    l =0; // level 
    pb[0]= b;
    pi[0]= b->n>0 ?(int)  b->n : N  ;
    pp[0]=p0;
    h=hb;
  L1: 
    do {   // walk on the tree  
      b= pb[l];
      while (pi[l]--) // loop on 4 element of the box
	{ 	      
	  int k = pi[l];
	  
	  if (b->n>0) // Vertex QuadTreeBox none empty
	    { 
	      Vertex * v=b->v[k];
	      if (v->ninside(P) ) {	    
		NbVerticesSearch++;
		Zd i2 =  VtoZd(v);
		// if good sens when try -- 
		h0 = Zd(i2,plus).norm();//  NORM(iplus,i2.x,jplus,i2.y);
		if (h0 <h) 
		{
		  h = h0;
		  vn =v;
		}}
	    }
	  else // Pointer QuadTreeBox 
	    { 
	      QuadTreeBox *b0=b;
	      NbQuadTreeBoxSearch++;
	      if ((b=b->b[k])) 
		{
		  hb >>=1 ; // div by 2
		  Zd ppp(pp[l],k,hb);
		
		  if  ( ppp.interseg(plus,hb,h) )//(INTER_SEG(iii,iii+hb,iplus-h,iplus+h) && INTER_SEG(jjj,jjj+hb,jplus-h,jplus+h)) 
		    {
		      pb[++l]=  b;
		      pi[l]= b->n>0 ?(int)  b->n : N  ;
		      pp[l]=ppp;		    
		    }
		  else
		    b=b0, hb <<=1 ;
		}
	      else
		b=b0;
	    }
	}
      hb <<= 1; // mul by 2 
    } while (l--);
    if (!vn && b != root )
   {// cas particulier on repart du sommet on avais rien trouver 
     b=root;
     hb =  MaxISize;
     p0=Zd();
     l=0;
     pb[0]= b;
     pi[0]= b->n>0 ?(int)  b->n : N  ;
     pp[0]=Zd();
     
     goto L1;
   }
    return vn;
}




  //  static int kfind=0;
  // static int kthrough=0;
  
  inline void CoorBary(const Triangle2 & K,const  R2  & P, R *l)
  {
    R detK = 2.*K.mesure() ;
    l[1]=det(K[0],P,K[2])/detK;
    l[2]=det(K[0],K[1],P)/detK;
    l[0]=1-l[1]-l[2];
  }
  

  inline void CoorBary(const Tet & K,const  R3  & P, R *l)
  {
    R detK = 6.*K.mesure() ;
    l[1]=det(K[0],P,K[2],K[3])/detK;
    l[2]=det(K[0],K[1],P,K[3])/detK;
    l[3]=det(K[0],K[1],K[2],P)/detK;
    l[0]=1-l[1]-l[2]-l[3];
  }
  
  
  inline    int nRand(int n) {
    return  rand()%n; //avant random()%n;
  }
  
  inline int find5(int i,int *k,int l)
  {
    if(l<5)
      {
	for(int j=0;j<l;++j)
	  if(i==k[j]) return j;
      }
    else  if(i>=k[0] && i<=k[l-1])
      {
	int i0=0,i1=l-1;
	while(i0<=i1)
	  { int im=(i0+i1)/2;
	    if(i<k[im]) i1=im-1;
	    else 
	      if(i>k[im]) i0=im+1;
	      else
		if(i==k[im]){
		  return im;
		}
	  }
      }
    return -1;
  } 
  
  template<class Mesh>
  const typename  Mesh::Element * Find(const Mesh & Th,
				       GTree<typename Mesh::Vertex> *quadtree,
				       typename Mesh::Rd P,
				       typename Mesh::RdHat & Phat,
				       bool & outside,
				       const typename  Mesh::Element * tstart)
  {
    typedef  typename Mesh::Element Element;
    typedef  typename Mesh::Vertex Vertex;
    typedef  typename Mesh::Rd Rd;
    const int nkv=Element::nv;
    const int d=Rd::d;
    int it,j;
    const int mxbord=100;
    int kbord[mxbord];
    int nbord=0;
    if ( tstart )
      it =  Th(tstart);
    else  
      {  
	const Vertex * v=quadtree->NearestVertexWithNormal(P);
	if (!v) 
	  { 
	  v=quadtree->NearestVertex(P);
	  assert(v);
	  }
	it=Th.Contening(v);
      }
    if(verbosity>200)
    cout << "tstart=" << tstart << " "<< "it=" << it << " P="<< P << endl; 

    //     int itdeb=it;     
    //     int count=0;
    //     L1: 
  outside=true; 
  //int its=it;
  //dPdP	int iib=-1;//,iit=-1;
  R dP=DBL_MAX;
  Rd PPhat;
  int k=0;    
  Mesh::kfind++;
  while (1)
    { 
      //if(verbosity>199) cout << "it " << it <<endl;
      const Element & K(Th[it]);
      Mesh::kthrough++;
      assert(k++<1000);
      int kk,n=0,nl[nkv];
      R l[nkv];
      for(int iii=0; iii<nkv; iii++)
	l[iii]=0.;

      CoorBary(K,P,l);

      // CoorBary :: donner entre 0 et 1
      // Pb si K.mesure*1.e-10 precision machine ==> bug
      
      // avant:
      // R eps =  -K.mesure()*1e-10;
      R eps = -1e-10;
      for(int i=0;i<nkv;++i)
	if( l[i] < eps){
	  nl[n++]=i;
	}
	if(verbosity>200){
      cout << "tet it=" << it << endl;
      cout << "K.mesure=" << K.mesure() ;
      cout << " eps=" << eps << endl;
      for(int i=0;i<nkv;++i)
	cout<< " l["<< i <<"]=" <<  l[i] << endl;
      cout << " n=" << n << endl;
	}
      if (n==0)
	{  // interior => return
	  outside=false; 
	  Phat=Rd(l +1);
	  // cout << Phat <<endl;
#ifndef NDEBUG
	  Rd pp=K(Phat)-P;
	  if( pp.norme()>1e-5)
	    {
	      cout << " Bug find P " << P << " ==" << K(Phat)  << endl;
	      cout << "Phat== " << Phat << " diff= " << pp << endl;
	      assert(0);
	    }

#endif	  
	  return &K;
	}
      
      

      kk=n==1 ? 0 : nRand(n);
      j= nl[ kk ];
      int itt =  Th.ElementAdj(it,j);
      if(itt!=it && itt >=0)  
	{
	  dP=DBL_MAX;
	  it=itt;
	  continue;
	}  
      int  inkbord=find5(it,kbord,nbord);
      if(inkbord<0 )  
	{
	  assert(nbord < mxbord);
	  kbord[nbord++]=it;
	  if(nbord>=5)  HeapSort(kbord,nbord);
	}
      if(verbosity>1)
	{
	cout << " bord "<< it<< "   nbf < 0 : " <<n << " (inb) " << inkbord << " nfb" << nbord<<endl;
	    R ss=0; 
	    for(int i=0;i<nkv;++i)
	      {  ss += l[i];
	      cout << l[i] << " ";}
	    cout << " s=" << ss << endl;;
	    
	}
	if(verbosity>200)
      cout << "GQuadTree::value of n " << n << endl;

      if ( n!=1 )  // on est sur le bord, mais plusieurs face <0 => on test les autre
	{  // 1) existe t'il un adj interne
	  int nn=0,ii;
	  int nadj[d+1],ni[d+1];
	  for(int i=0;i<nkv;++i)
	    if (l[i] < eps && (itt=Th.ElementAdj(it,ii=i)) != it && (itt>=0) && find5(itt,kbord,nbord) < 0 ) 
	      ni[nn++]=i,nadj[i]=itt;
	  if(verbosity>100)
	    cout << " nn : "<< nn << endl;
	  if (nn>0)
	    {
	      //j=nadj[nRand(nn)];
	      j=ni[nRand(nn)];
	      it=nadj[j];
	      dP=DBL_MAX;
	      cout << "new it= " << it << endl;
	      continue;
	    }
	}
      // toutes les faces <0  sont sur le bord.
      //   pour l'instant on s'ar阾e
      // le point est externe. (mais trop cher pour faire mieux)
      // on projet le points sur le bord via le coordonne 
      //  barycentrique 
      { // a ridge on border  (to hard to do the correct stuff) 
	//  or a corner    just do the projection on lambda  
	R s=0.;
	for(int i=0;i<nkv;++i)
	  s += (l[i]= max(l[i],0.));
	for(int i=0;i<nkv;++i)
	  l[i]/=s;
	Phat=Rd(l +1);
	if(verbosity>100)
	  {
	    cout << P << " " << n << " l: ";
	      R ss=0; 
	  for(int i=0;i<nkv;++i)
	    {  ss += l[i];
	    cout << l[i] << " ";}
	  cout << " s=" << ss <<" " << s <<" exit by bord " << it << " "  << Phat << endl;;
          }
	outside=true;
	return &Th[it] ;
      }		    
    }
}
// Instantiation du manuel des templates

  template class GTree<Vertex2>;
  template class GTree<Vertex3>;
  template class GTree<Vertex1>;
  typedef Mesh3::GMesh GMesh3;
  typedef Mesh2::GMesh GMesh2;
  typedef Mesh1::GMesh GMesh1;
  template
  const   GMesh3::Element * Find<GMesh3>(const GMesh3 & Th,
				       GTree< GMesh3::Vertex> *quadtree,
				       GMesh3::Rd P,
				       GMesh3::RdHat & Phat,
				       bool & outside,
				       const   GMesh3::Element * tstart);
  template
  const   GMesh2::Element * Find<GMesh2>(const GMesh2 & Th,
				       GTree< GMesh2::Vertex> *quadtree,
				       GMesh2::Rd P,
				       GMesh2::RdHat & Phat,
				       bool & outside,
				       const   GMesh2::Element * tstart);

/*
  template
  const   GMesh1::Element * Find<GMesh1>(const GMesh1 & Th,
				       GTree< GMesh1::Vertex> *quadtree,
				       GMesh1::Rd P,
				       GMesh1::RdHat & Phat,
				       bool & outside,
				       const   GMesh1::Element * tstart);
*/
}