genericmesh.hpp

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

	  R lb[Rd::d+1];//{1.-PHat.sum(),PHat}; 
	  R lbb[Rd::d+1];//{1.-PHat.sum(),PHat}; 
	  PHat.toBary(lb); // R1 R2 R3 
	  if(Abs(lb[j])>1e-10)
	   assert(Abs(lb[j])<1e-10);
	int sigma[T::nva];	  
	const void * nvkj[T::nva], *nvkkjj[T::nva];
	int jj=p%nea;
	int kk=p/nea;

	Element & K(elements[CheckT(k)]);
	Element & KK(elements[CheckT(kk)]);
	Rd Pin=K(PHat);  
	for (int l=0;l<T::nva;++l)
	    nvkj[l] =&K[T::nvadj[j][l]];
	for (int l=0;l<T::nva;++l)
	    nvkkjj[l] = &KK[T::nvadj[jj][l]];
	//  il faut permute ll.
	PermI2J<nva>(nvkj,nvkkjj,sigma);
	for (int l=0;l<T::nva;++l)
	    lbb[T::nvadj[jj][l]]=lb[T::nvadj[j][sigma[l]]];
	lbb[jj]=0;
#ifdef DEBUG	  
	Rd PH=PHat;  
#endif
	PHat=Rd(lbb+1);
#ifdef DEBUG	  
	Rd Pout=KK(PHat);
	if( (Pin-Pout).norme2() > 1e-10 )
	    {
		for (int l=0;l<=T::nva;++l)
		    cout << lbb[l] <<" < -- " << lb[l] << endl;
		for (int l=0;l<T::nva;++l)
		    cout <<l << " :    o=  " << nvkkjj[l]  << "   i= " << nvkj[l] << " " <<  sigma[l] 
		         << " -- " << &KK[T::nvadj[jj][l]]  << " == " << &K[T::nvadj[j][sigma[l]]] 
		    << " -- " << &K[T::nvadj[j][l]]  << " == " << &KK[T::nvadj[jj][sigma[l]]] 
		    << " -- " << lbb[T::nvadj[jj][l]] << " == " << lb[T::nvadj[j][sigma[l]]]
		    << " ++ " << T::nvadj[jj][l] << " <-- " << T::nvadj[j][sigma[l]] 
		    << endl;
		cout << "Adj:  j= " << j << " ," << Pin << " != " << Pout << " , " << PH << " -> " << PHat << "  jj = " << jj <<  endl;
		assert(0);
	    }
#endif	  
	j=jj;
	return kk;
      }
    return -1;//  on border 
  }
    
  int GetAllElementAdj(int it,int *tabk) const
  { //  get the tab of all adj element (max ne)
    //  and return the size of the tab 
    int i=0;
    for(int j=0;j<nea;++j)
      {
	tabk[i]=TheAdjacencesLink[3*it+j]/3;
	if(tabk[i] >=0 && tabk[i]!=it) i++; 
      }
    return i;
  }

  int BoundaryElement(int be,int & ItemInK) const {
    int i= BoundaryElementHeadLink[be]; 
    ItemInK = i%nea; 
    return i/nea;}
  
  // Add J. Morice 
  template<int N,int M>
  SortArray<int,N> itemadjs(const int (* const  nu )[N],int k,int i, int *sens) 
  {
    int nv[N];
    B & K(borderelements[CheckBE(k)]);
    ASSERTION(i>=0 && i <M);
    for (int j=0;j<N;++j){
      nv[j] = operator()(K[nu[i][j]]);
    }
    if(nv[0] > nv[1] )
      *sens = 1;
    else
      *sens =-1;
    return SortArray<int,N>(nv);
  }

  SortArray<int,B::nva> items(int k,int i,int *sens) 
  {
    return itemadjs<B::nva,B::nv>(B::nvadj,k,i,sens);
  }

  
  template<int N,int M>
  SortArray<int,N> iteme(const int (* const  nu )[N],int k,int i) 
  {
    int nv[N];
    Element & K(elements[CheckT(k)]);
    ASSERTION(i>=0 && i <M);
    for (int j=0;j<N;++j){
      nv[j] = operator()(K[nu[i][j]]);
    }

    return SortArray<int,N>(nv);
  }

  SortArray<int,B::nv> itemadj(int k,int i) 
  {
    return iteme<B::nv,T::nea>(T::nvadj,k,i);
  }
  
  SortArray<int,B::nv> itembe(int k) 
  {
    int nv[B::nv];
    B & K(borderelements[CheckBE(k)]);
    
    for (int j=0;j<B::nv;++j){
      nv[j] = operator()(K[j]);
    }

    return SortArray<int,B::nv>(nv);
  }

  //  const Element * Find(const Rd & P) const ;
  const Element * Find(Rd P, RdHat & Phat,bool & outside,const Element * tstart=0) const  
  {return EF23::Find<GMesh>(*this,this->gtree,P,Phat,outside,tstart);}
  
  R mesure(){ return mes;}
  R bordermesure(){ return mesb;}
  virtual ~GenericMesh() { 
    cout << "~GenericMesh\n";
   
    delete [] ElementConteningVertex;
    delete [] TheAdjacencesLink;
    delete [] BoundaryElementHeadLink;
    delete [] borderelements;
    if(nt>0) delete [] elements;
    delete [] vertices;
    delete [] bnormalv;
    if(gtree) delete gtree;
    
  }

  Serialize serialize() const;
  
private:
  GenericMesh(const GenericMesh &); // pas de construction par copie
   void operator=(const GenericMesh &);// pas affectation par copy 
};

template<typename T,typename B,typename V>
void GenericMesh<T,B,V>::BuildjElementConteningVertex()
{
  const int nkv= T::nv;
  if(!ElementConteningVertex) ElementConteningVertex = new int[nv];

    for(int i=0;i<nv;++i)
	ElementConteningVertex[i]=-1; 
    
    for (int k=0;k<nt;++k)
	for (int i=0;i<nkv;++i)
	    ElementConteningVertex[operator()(elements[k][i])]=k ;
    int kerr=0;
    for(int i=0;i<nv;++i)
	if (ElementConteningVertex[i]<0) 
	    kerr++; 
    assert(kerr==0);

} 
template<typename T,typename B,typename V>
void GenericMesh<T,B,V>::BuildAdj()
{
  const int nva   = T::nva;
  const int nea   = T::nea;
  assert(TheAdjacencesLink==0);
  TheAdjacencesLink = new int[nea*nt];
  BoundaryElementHeadLink = new int[nbe];
  HashTable<SortArray<int,nva>,int> h(nea*nt,nv);
  int nk=0,nba=0;
  int err=0;

  cout << "nva=// nea=" << nva << " " << nea << " "<< nbe << endl;
  for (int k=0;k<nt;++k)
    for (int i=0;i<nea;++i)
      {
        SortArray<int,nva> a(itemadj(k,i));
	//cout << " ### "   << " item(k,i)= " << itemadj(k,i) << " a= " << a << " k " << k << " i " << i << endl;
	typename HashTable<SortArray<int,nva>,int>::iterator p= h.find(a);
	if(!p) 
	  { 
	    h.add(a,nk);
	    TheAdjacencesLink[nk]=-1;
	    nba++;
	  }
	else 
	  {	  
	    ASSERTION(p->v>=0);
	    TheAdjacencesLink[nk]=p->v;
	    TheAdjacencesLink[p->v]=nk;
	    p->v=-1-nk;
	    nba--;
	  }
	++nk;
      }
    
  for (int k=0;k<nbe;++k)
     {
	SortArray<int,nva> a(itembe(k));

	typename HashTable<SortArray<int,nva>,int>::iterator p= h.find(a);
	//cout << k << " ### "   << " item(k,i)= " << itembe(k) << " a= " << a << endl;
	if(!p) { err++;
	if(err==1) cerr << "Err  Border element not in mesh \n";
	if (err<10)  cerr << " \t " << k << " " << a << endl;
	}
	 else
	   {
	     BoundaryElementHeadLink[k] = p->v <0 ? -p->v-1 : p->v;
	     #ifndef NDEBUG
	     int t=BoundaryElementHeadLink[k]/nea;
	     int e=BoundaryElementHeadLink[k]%nea;
	     //cout << k << " ### "   << a << " = " << itemadj(t,e) << " t " << t << " e " << e << endl;
	     assert(itemadj(t,e)==a);
	     #endif
	   }
     }

    
  assert(err==0);
  int na= h.n;
  cout << "  -- Nb adj  = "<< na << " on border " << nba << " nea = " << nea << " nva = " << nva ;
  if(nea==2)
    cout << " Const d'Euler: " << nt - na + nv << endl;
  else
    cout << endl;
}

template<typename T,typename B,typename V>
void GenericMesh<T,B,V>::BuildSurfaceAdj()
{
 
  // assert(TheSurfaceAdjacencesLink==0); plus tard
  int *TheSurfaceAdjacencesLink = new int[B::nea*nbe];
  HashTable<SortArray<int,B::nva>,int> h(B::nea*nbe,nv);
  int nk=0;
  int err=0;
  int sens;

  cout << "nea/nva" << B::nea << " "  << B::nva << endl;
  for (int k=0;k<nbe;++k)
    for (int i=0;i<B::nea;++i)
      {
        SortArray<int,B::nva> a(items(k,i,&sens));

	typename HashTable<SortArray<int,B::nva>,int>::iterator p= h.find(a);
	if(!p) 
	  { 
	    h.add(a,nk);
	    TheSurfaceAdjacencesLink[nk]=sens;
	  }
	else 
	  {	  
	    ASSERTION(p->v>=0);
	    if( sens*TheSurfaceAdjacencesLink[p->v] != -1){
	 
	      B & K(borderelements[CheckBE(k)]);
	      int firstVertex  =  operator()(K[B::nvadj[i][0]])+1;
	      int secondVertex =  operator()(K[B::nvadj[i][1]])+1;
	      cout << " The edges defined by vertex is " << firstVertex << "-" << secondVertex << " is oriented twicd in the same direction "<< endl;
	      err++;
	    }
	    TheSurfaceAdjacencesLink[nk]=p->v;
	    TheSurfaceAdjacencesLink[p->v]=nk;    
	  }
	if( err > 10 ) 
	  exit(1); 
	nk++;
      }
    
  assert(err==0);
  delete [ ] TheSurfaceAdjacencesLink; 
  cout << "number of adjacents edges " << nk << endl; 
}



template<typename T,typename B,typename V>
DataFENodeDF GenericMesh<T,B,V>::BuildDFNumbering(int ndfon[NbTypeItemElement]) const
{
  const GenericMesh & Th(*this);
  int nnodeK = T::NbNodes(ndfon);
  int *p = 0, *pp=0; 
  unsigned int tinfty=-1;
  const int nkv= T::nv;
  const int nkf= T::nf;
  const int nke= T::ne;
  const int nkt= T::nt;
  const int nk[]={nkv,nke,nkf,nkt};
  int MaxNbNodePerElement=0;
  int MaxNbDFPerElement=0;
  int nbNodes=0;
  int NbOfDF=0;
  int n=0;
  int minndf=100000000;
  int maxndf=0;
  int nbnzero=0;
  for (int dd=0;dd<NbTypeItemElement;++dd)
    if(ndfon[dd])
      {
	
        nbnzero++;
	minndf=Min(minndf,ndfon[dd]);
	maxndf=Max(maxndf,ndfon[dd]);
	MaxNbDFPerElement   += ndfon[dd]*nk[dd];
	MaxNbNodePerElement += nk[dd];
      }
  bool constndfpernode = minndf == maxndf;
  bool nodearevertices = nbnzero ==1  && ndfon[0];
  assert(maxndf>0);
  const int nkeys=4+6+4+1;
  assert(nnodeK<= nkeys);
  if(nodearevertices)
    {
      nbNodes=nv;
      NbOfDF=nbNodes*ndfon[0];
    }
  else
    {      
      //  construction of nodes numbering 
      p =  new int[nnodeK*nt];
      typedef SortArray<unsigned int,2> Key;
      // ------------
      Key keys[nkeys];
      int keysdim[nkeys];
      int of = Th.nv+1;
      int ndim[NbTypeItemElement]={0,0,0,0};
      NbOfDF=0;
      {
	HashTable<Key,int> h(nnodeK*nt,of+nt);
	for(int k=0;k<nt;++k)
	  {
	    const T & K(Th[k]);
	    int m=0;
	    if( ndfon[0] )//  node on vertex
	      for(int i=0;i<nkv;++i)
		keysdim[m]=0,keys[m++]=Key(Th(K[i]),tinfty);
	    if(  ndfon[1] )//  node on Edge
	      for(int i=0;i<nke;++i)
		keysdim[m]=1,keys[m++]=Key(Th(K[T::nvedge[i][0]]),Th(K[T::nvedge[i][1]]));
	    if(  ndfon[2])//  node on Face
	      if (nkf==1) 
		keysdim[m]=2,keys[m++]=Key(k+of,tinfty);
	      else
		for(int ii,i=0;i<nkf;++i)
		  keysdim[m]=2,keys[m++]=Key(k+of,ElementAdj(k,ii=i) +of);
	    if(  ndfon[3] )//  node on Tet
	      if(nkt==1)
		keysdim[m]=3,keys[m++]=Key(k+of,tinfty);
	    
	    if(k<0)
	    {
	      for(int i=0;i<nke;++i)
		cout << " e= " << T::nvedge[i][0] << " " << T::nvedge[i][1] << endl;
	      cout << ndfon[0] << " " << ndfon[1] << " " << ndfon[2] << " " << ndfon[3] << ": "
		   <<  " m = "<< m << "  " <<nnodeK
		   << " " << T::nv
		   << " " << T::ne 
		   << " " << T::nf 
		   << " " << T::nt  
		   << endl;
	    }
	    assert(m==nnodeK);
	    for(int i=0;i<m;i++)
	      {
		typename HashTable<Key,int>::iterator pk= h.find(keys[i]);
		if(!pk)
		 {
		  pk = h.add(keys[i],nbNodes++);
		  NbOfDF += ndfon[keysdim[i]];
		 }
		p[n++] = pk->v;
		ndim[keysdim[i]]++;
			  
	      }
	  }
      }
      cout << " nb of Nodes " << nbNodes << endl;
      cout << " nb of DoF   " << NbOfDF << endl ;
      if( ! constndfpernode) 
	{ 
	  pp=new int[nbNodes+1];
	  int kk=0,nn=0;
	  for(int k=0;k<nt;++k)
	    for(int i=0;i<nnodeK;i++)
		pp[p[nn++]]=ndfon[keysdim[i]];;
	  for(int n=0;n<nbNodes;++n)	
	    {  
	      int ndfn=pp[n];
	      pp[n]=kk;
	      kk+=ndfn;
	    }
	  pp[nbNodes]=NbOfDF;//  add last 
	  assert(kk==NbOfDF);
	}
    }
  return DataFENodeDF(ndfon,nt,nbNodes,NbOfDF,p,pp,MaxNbNodePerElement,MaxNbDFPerElement);
}
template<typename T,typename B,typename V>
void GenericMesh<T,B,V>::BuildBound()
{
    if(vertices && (nv>0))
    {
	Pmin=vertices[0];
	Pmax=vertices[0];
	 for(int i=1;i<nv;++i)
	 {
	     Pmin=Minc(Pmin,vertices[i]);
	     Pmax=Maxc(Pmax,vertices[i]);

	 }
    }
	
}

template<typename T,typename B,typename V>
void GenericMesh<T,B,V>::Buildbnormalv()
{
    const int nkv= T::nv;
   // const int nkf= T::nf;
   // const int nke= T::ne;
   // const int nkt= T::nt;
    
    if (bnormalv) 
      {assert(0);return;}
    int nb=0;
    for (int k=0;k<nt;k++)
	for (int i=0;i<nkv;i++)
	{  
	    int ii(i),kk;
	    kk=ElementAdj(k,ii);
	    if (kk<0 || kk==k) nb++;
	}
    if(verbosity>2)
	cout << " number of real boundary  " << nb << endl;
    bnormalv= new Rd[nb];
    Rd *n=bnormalv;
    for (int k=0;k<nt;k++)
	for (int i=0;i<nea;i++)
	{  
	    int ii,kk=ElementAdj(k,ii=i);
	    if (kk<0 || kk==k) {
		Element & K(elements[k]);
		Rd N;//=K.n(i);
		for(int j=0;j<nva;++j)
		{
		    K[Element::nvadj[i][j]].SetNormal(n,N);
		}
		    
	    }
	}
    // cout << nb << " == " << n-bnormalv << endl;
    assert(n - bnormalv <= nb );
}

}
#endif