fem.cpp

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

	nv =0;
	quadtree = new FQuadTree(this,Pmin,Pmax,nv); //  put all the old vertices in the quadtree 
						     //  copy of the old vertices 
	for (int i=0;i<Th.nt;i++)
	    if (split[i]) 
		for (int j=0;j<3;j++)
		{
		    
		    Vertex * pV=quadtree->ToClose(Th[i][j],seuil);
		    if (pV ==0) { 
			// cout << Th(i,j) << " " ;
			vertices[nv] = Th[i][j];
			vertices[nv].normal=0;
			quadtree->Add(vertices[nv]);
			nv++;}
		}
		    // nv = Th.nv;
		    if(verbosity>3) 
		    {
			cout << "  -- number of old vertices use: " << nv << endl;      
			cout << "  -- number of  neb : " << nebmax << endl;
		    }
	bedges = new BoundaryEdge[nebmax];
	BoundaryEdge * bedgesi= new BoundaryEdge[nebimax];  // jan 2007 FH
	int * sdd= new int[Th.nt]; 
	for (int i=0;i<Th.nt;++i)
		sdd[i]= 0; 
	/*   
	    for (int i=0;i<Th.neb;i++)
	{
		cout << i << " " << " " << Th(Th.bedges[i][0]) << " " << Th(Th.bedges[i][1]) << endl;
	} */
	assert(bedges && bedgesi);
	//  generation of the boundary edges
	neb =0;
	nebi=0;   //   generaztion des arete interne pour les force ...   
	nbsdd=0;
	//   for (int ieb=0,jj;ieb<Th.neb;ieb++)
	for (int it=0;it<Th.nt;it++)
	    if (  split[it] ) 
	    {
		// sdd[it]=-1;
		for (int jt=0;jt<3;jt++)
		{
		    int jtt=jt,itt=Th.ElementAdj(it,jtt);
		    //  cout << it <<  " " << jt << " " << jt << " " << itt << !split[itt] << endl;
		    int ie0,ie1;
		    Label  re(label); 
		    Th.VerticesNumberOfEdge(Th[it],jt,ie0,ie1);
		    //  cout << "++ ";
		    BoundaryEdge * pbe = Th.TheBoundaryEdge(ie0,ie1);
		    bool bbe= ( itt == it || itt <0 || !split[itt] || (pbe && &(*pbe)[0] == &Th(ie0))) ;
		    
		    BoundaryEdge * bbedges = bbe ? bedges : bedgesi;
		    int &  nneb = bbe ? neb : nebi;
		    int nnebmax = bbe ? nebmax : nebimax;
		    int offset= bbe ? 0 : nebmax;
		    if (bbe ||  (!pbe && (ie0 < ie1) ) ) // arete interne ou frontiere 
		    {
			int sens = 1; // par defaul le bon sens 
			int kold = it;   //Th.BoundaryElement(ieb,jj);
			int n=split[kold];
			if( itt>=0) n = max(n,split[itt]); //  pour les aretes internes (FH juillet 2005)
			if (!n) continue; 
			if (pbe ) {
			    re = *pbe;
			    if( & (*pbe)[0] == &Th(ie1) ) sens = -sens; // pour les aretes non decoupe avril 2007
			   // cout << " ### " << ie0 << " " << ie1 << " " <<  sens << endl;
			}
			// cout << " lab = " <<  re.lab << endl;
			Vertex *pva= quadtree->NearestVertex(Th(ie0));
			Vertex *pvb= quadtree->NearestVertex(Th(ie1));
			Vertex *pv0=pva;
			R2 A(*pva),B(*pvb);
			R la = 1,lb=0, delta=1.0/n;
			
			for (int j=1;j<n;j++) 
			{ 
			    sens = 1; //  arete decoupe => le sens change avril 2007
			    la-=delta;
			    lb+=delta;
			    assert(nv<nvmax);
			    Vertex *pv1= vertices + nv;
			    
			    (R2 &) *pv1 = A*la + B*lb;
			    (Label &) *pv1 = re ; //= (Label) be;
			    quadtree->Add(*pv1);
			    nv++;
			    
			    assert(nneb<nnebmax);
			    bbedges[nneb].vertices[0]=pv0;
			    bbedges[nneb].vertices[1]=pv1;
			    (Label &)  bbedges[nneb] = re;
			    nneb++;
			    
			    pv0=pv1;
			}
			//if (bbe)
			 //   cout << nneb+1 << " yyyy  t" << it << " a=" << jt << " taj=" << itt << " sa =" <<    ie0 << " " << ie1 << " + = " << n << endl;

			assert(nneb<nnebmax);
			bbedges[nneb].vertices[0]=pv0;
			bbedges[nneb].vertices[1]=pvb;
			(Label &) bbedges[nneb]= re ;		         			    
			sdd[it]= (1+ (nneb++ + offset))*sens; // numero de la derniere arete attention au sens si pas decoupe avril 2007
			//cout << " ### " << pv0 - vertices << " " << pvb - vertices << " " <<  sens << " it = " << it << " " <<  sdd[it] << " itt " << itt  << " -- " << sdd[406] << endl;

			if(  ( itt !=it || itt <0)  ) // interne 
			   sdd[itt]= -sdd[it]; 
		    }
		} 
		nbsdd++; 
	    }

		//   cout << "          (debug)  nb vertices on egdes " << nv << endl;  
		    //  cout << " " <<  nebmax << " " << neb << endl;
	assert(neb==nebmax);
	assert(nebi==nebimax);
	assert(nbsdd==nbsddmax);
	
	//   create the new vertices and the new triangle 
	int   kt=0;
	for (int K=0;K<Th.nt;K++)
	{ // cout << K << endl;
	    Triangle &T(Th[K]);
	    R2 A(T[0]);
	    R2 B(T[1]);
	    R2 C(T[2]);
	    long N=split[K];
	    if (!N) continue;
	    long N2=N*N;
	    int vt[3];
	    for (int n=0;n<N2;n++,kt++) //  loop on all sub triangle
	    {
		//(Label&) triangles[kt] = (Label&) T; 
		for(int j=0;j<3;j++) //  Loop on the 3 vertices
		{ 
		    R2 PTj=SubTriangle(N,n,j);
		    R la=1-PTj.x-PTj.y,lb=PTj.x,lc=PTj.y;
		    R lmin = Min(la,lb,lc);
		    R2 Pj= A*la+B*lb+C*lc; 
		    Vertex *pV;
		    pV=quadtree->ToClose(Pj,seuil);
		    // if !noregenereration => point du bord du triangle deja genere 
		    bool addv = !pV;
		    if(!noregenereration && pV==0) addv = lmin > 1e-5;
		    if ( addv )			
		    { // new vertex
		      // cout << "    -- " << nv << "  New Vertices " << Pj << " n=" << n << " j=" << j << " " << PTj << endl;
			assert(nv<nvmax);
			vertices[nv]=Pj;
			(Label&) vertices[nv]=0; //  Internal vertices 
			pV = vertices + nv;
			quadtree->Add(*pV);
			nv++;                  
		    }  //  end of new vertex
		    if(noregenereration)
		    {
		      assert(pV);
		    //cout << j <<  "  " << *pV << " " << Pj << " " << number(pV) << " " << Norme2(Pj-*pV) << " " << nv << endl;
		      vt[j]=number(pV);
		    }
		}
		   
		    //  triangles[kt].SetVertex(j,pV);
		 // for(int j=0;j<3;j++)
		  // cout << kt << " " << n << endl;
		if(noregenereration)
		{
		    R2 A=vertices[vt[0]];
		    R2 B=vertices[vt[1]];
		    R2 C=vertices[vt[2]]; 
		    R a = (( B-A)^(C-A))*0.5;
		    
		    if (a>0) 
			triangles[kt].set(vertices,vt[0],vt[1],vt[2],T.lab);
		    else 
			triangles[kt].set(vertices,vt[0],vt[2],vt[1],T.lab);
		}
	    }   // end loop on all sub triangle
	    
	} //  end 
	if (verbosity>3 ) 
	{ 
	    cout << "  - regeneration = " << ! noregenereration <<endl; 
	    cout << "  - Nb of vertices       " << nv << endl; 
	    cout << "  - Nb of triangle       " << nt << endl; 
	    cout << "  - Nb of boundary edges " << neb << endl; 
	
	}
	//  
	if (!noregenereration )   // REGENRATION DU MAILLAGE 
	{ // 
	   // nebi=0;
	    long nba = neb+nebi;
	    // 
	    long nbsd = nbsdd; // bofbof 
			   //ok,  with correction of mshptg 
	    long  *sd;
	    sd=new long[2*nbsd+2];
	    
	    sd[0]=-1;
	    sd[1]=Th[0].lab;         

	    
	    long nbs=nv;
	    long nbsmax=nv;
	    long            err = 0;//, nbsold = nbs;       
		long           *c = 0;
		long           *tri = 0;
		long           *nu = 0;
		long           *reft = 0;
		typedef double Rmesh ;
		Rmesh          *cr = 0;
		Rmesh          *h = 0;
		long nbtmax = 2 * nbsmax;
		long * arete  = new long[2*nba]; 
		nu = new long[6*nbtmax];
		c = new long[2*nbsmax];
		tri = new long[Max((4 * nbsmax + 2 * nbsd),nba)];
		reft = new long[nbtmax];
		cr = new Rmesh[(2 * nbsmax + 2)];
		h = new Rmesh[nbsmax];
		for(int i=0,k=0; i<nv; i++)
		{ 
		    cr[k++]  =vertices[i].x;
		    cr[k++]=vertices[i].y;
		    h[i]=1; 
		}
		{
		    int k=0;
		for (int i=0 ;i<neb;i++)
		{
		    arete[k++] =number(bedges[i][0])+1;
		    arete[k++] =number(bedges[i][1])+1;
		}
		//  ajoute des aretes interne pour les forces ..  FH 
		for (int i=0 ;i<nebi;i++)
		{
		    arete[k++] =number(bedgesi[i][0])+1;
		    arete[k++] =number(bedgesi[i][1])+1;
		}
		}
		// construction du tableau sd
		
		for (int it=0,j=0;it<Th.nt;it++)
		    if (  split[it]) // un ssd par vieux triangle 
		    {
			assert(sdd[it]);
			sd[j++]=sdd[it];
			sd[j++] = Th[it].lab; 
		    }
			
		extern int 
		    mshptg8_ (Rmesh *cr, Rmesh *h, long *c, long *nu, long *nbs, long nbsmx, long *tri,
			      long *arete, long nba, long *sd,
			      long nbsd, long *reft, long *nbt, Rmesh coef, Rmesh puis, long *err);
		
		long nbt=0;
		if(verbosity>10)
		{
		    cout << " mshptg8_ " << endl;
		    cout << "    nbs =" << nbs << endl;
		    cout << "    nba =" << nba << endl;
		    cout << "    nbt =" << nbt << endl;
		    cout << "    nbsd =" << nbsd << endl;
		    cout << " sommets : " << endl;
		    for(int i=0;i<nbs; ++i)
			cout << " " << i+1 << " -- " << cr[2*i] << " " << cr[2*i+1] << endl;
		    cout << " aretes : " << endl;
		    for(int i=0;i<nba; ++i)
			cout << " " << i+1 << " -- " << arete[2*i] << " " << arete[2*i+1] << endl;
		    cout << " Sd : " << endl;
		    for(int i=0;i<nbsd; ++i)
			cout << " " << i+1 << " -- " << sd[2*i] << "  lab " << sd[2*i+1] << endl;
		    //nbsd=0;
		}
		mshptg8_ (cr, h, c, nu, &nbs, nbs, tri, arete, nba, (long *) sd, nbsd, reft, &nbt, .25, .75, &err);
		if( !(err==0 && nbt !=0))
		{
		    cerr << " Error mesh generation in mshptg : err = " << err << " nb triangles = " << nbt << endl;
		    ffassert(err==0 && nbt !=0);
		}
		// 
		
		delete [] triangles;
		// Correction FH bug  trunc  mesh with hole 25032005 +july 2005 
		int kt=0;
		R dmin=1e100;
		for(int i=0,k=0;i<nbt;i++,k+=3)
		{
		    //int ir = reft[i]; 
		    reft[i]=-1;
		    R2 A=vertices[nu[k]-1],B=vertices[nu[k+1]-1],C=vertices[nu[k+2]-1];
		    R2 G=(A+B+C)/3.,PHat;
		    double d=Area2(A,B,C);
		    dmin=min(d,dmin);
		    //cout << A << "\n" << B << "\n" << C << "\n" << A << "\n\n";
		    if(d<=1e-5 && 0)
		    {
			cout<< " T = "<<  i << " det= " << d << "  ::  " << A << " " << B << " " << C  << endl;
			
		    }
		    bool outside;
		    const Triangle * t=Th.Find(G,PHat,outside,0);
		    if(!outside ) { 
			int k=Th(t);
			if( split[k] )
			{
			    reft[i] = k;
			    kt++;
			}
		    }
		}
		nt=kt;
		// cout << " " << dmin << endl;
		if(verbosity>3)      
		    cout << "  - Nb Triangles = " << nt <<  " remove triangle in hole :" <<  nbt - nt 
			<<endl;
		triangles = new Triangle[nt];
		kt=0;
		for(int i=0,k=0;i<nbt;i++,k+=3)
		    if(reft[i]>=0) 
			triangles[kt++].set(vertices,nu[k]-1,nu[k+1]-1,nu[k+2]-1,Th[reft[i]].lab);
		assert(kt==nt);
		// END  Correction FH bug  trunc  mesh with hole 25032005 + july 2005 
		
		delete [] arete;
		delete [] nu;
		delete [] c;
		delete [] tri;
		delete [] reft;
		delete [] cr;
		delete [] h;
		delete [] sd;
		
	}
	delete [] bedgesi;   
	delete [] sdd;   
	ConsAdjacence();
}


const char Mesh::magicmesh[8]="Mesh 2D";
int  Mesh::kthrough =0;
int  Mesh::kfind=0;

Mesh::Mesh(const  Serialize &serialized) 
{
    TriangleConteningVertex =0;
    BoundaryAdjacencesHead=0;
    BoundaryAdjacencesLink=0;
    BoundaryEdgeHeadLink=0;
    quadtree =0;
    volume=0;
    NbMortars=0;
    dim=0;
    tet=0;
    edges=0;
    ntet=0;
    ne=0;
    
    mortars=0;
    TheAdjacencesLink =0;
    nv=0;
    nt =0;
    neb=0;
    triangles=0;
    vertices=0;
    bedges=0;
    area=0;
    bnormalv=0;
    //  ---  assert(serialized.samewhat(magicmesh));
    size_t  pp=0;;
    long long l;
    serialized.get(pp,l);
    serialized.get( pp,nt);
    serialized.get( pp,nv);
    serialized.get( pp,neb);
    if (verbosity>2) 
	cout << " mesh serialized : " << l << " " << nt << " " << nv << " " << neb << endl;
    assert ( nt > 0 && nv >0 && neb >=0);
    triangles = new Triangle [nt];
    vertices  = new Vertex[nv];
    bedges    = new BoundaryEdge[neb];
    area=0;
    ffassert(triangles && vertices && bedges);
    
    for (int i=0;i<nv;i++)    
    {
        serialized.get(pp,vertices[i].x);
        serialized.get(pp,vertices[i].y);
        serialized.get(pp,vertices[i].lab);
    }
    for (int i=0;i<nt;i++) { 
        int i0,i1,i2,ir;
        serialized.get(pp,i0);
        serialized.get(pp,i1);
        serialized.get(pp,i2);
        serialized.get(pp,ir);
	
        triangles[i].set(vertices,i0,i1,i2,ir); 
        area += triangles[i].area;}
    
    for (int i=0;i<neb;i++) { 
        int i0,i1,ir;
        serialized.get(pp,i0);
        serialized.get(pp,i1);
        serialized.get(pp,ir);
        bedges[i] = BoundaryEdge(vertices,i0,i1,ir); }
    assert( pp ==  serialized.size() );   
    if(verbosity>2) 
	cout << "   End of un serialize: area on mesh = " << area <<endl;  
    ConsAdjacence();
    
}

Serialize  Mesh::serialize() const
{
    
    long long  l=0;
    l += sizeof(long long);
    l += 3*sizeof(int);
    l += nt*(sizeof(int)*4);
    l += nv*( sizeof(int) + sizeof(double)*2);
    l += neb*(sizeof(int)*3);
    
    // cout << l << magicmesh << endl;
    Serialize  serialized(l,magicmesh);
    // cout << l << magicmesh << endl;
    size_t pp=0;
    serialized.put(pp, l); 
    serialized.put( pp,nt);
    serialized.put( pp,nv);
    serialized.put( pp,neb);
    if (verbosity>2) 
      cout << " mesh Serialized : " << l << " "  << nt << " " << nv << " " << neb << endl;
    for (int i=0;i<nv;i++)
    {
	serialized.put(pp, vertices[i].x);
	serialized.put(pp, vertices[i].y);
	serialized.put(pp, vertices[i].lab);
    }
    //cout << " ---  " << endl;
    for (int i=0;i<nt;i++)
    {
	const Triangle & K(triangles[i]);
	int i0= &K[0]-vertices;
	int i1= &K[1]-vertices;
	int i2= &K[2]-vertices;
	serialized.put(pp, i0);
	serialized.put(pp, i1);
	serialized.put(pp, i2);
	serialized.put(pp, K.lab);
    }
    // cout << " ---  " << endl;
    for (int i=0;i<neb;i++)
    {
	const BoundaryEdge & K(bedges[i]);
	int i0= &K[0]-vertices;
	int i1= &K[1]-vertices;
	serialized.put(pp, i0);
	serialized.put(pp, i1);
	serialized.put(pp, K.lab);
    }
    //  cout << " ---  " << endl;
    // cout << pp << " == " << serialized.size() << endl;
    assert(pp==serialized.size());
    return serialized; 
}

void Mesh::Buildbnormalv()
{
    if (bnormalv) 
    {assert(0);return;}
    int nb=0;
    for (int k=0;k<nt;k++)
	for (int i=0;i<3;i++)
	{  
	    int ii(i),kk;
	    kk=ElementAdj(k,ii);
	    if (kk<0 || kk==k) nb++;
	}
	    if(verbosity>2)
		cout << " number of real boundary edges " << nb << endl;
    bnormalv= new R2[nb];
    R2 *n=bnormalv;
    for (int k=0;k<nt;k++)
	for (int i=0;i<3;i++)
	{  
	    int ii(i),kk;
	    kk=ElementAdj(k,ii);
	    if (kk<0 || kk==k) {
		Triangle & K(triangles[k]);
		R2 N=K.n(i);
		Vertex & v0= K.Edge(0,i);
		Vertex & v1= K.Edge(1,i);
		v0.SetNormal(n,N);
		v1.SetNormal(n,N);         
	    }
	}
	    // cout << nb << " == " << n-bnormalv << endl;
	    assert(n - bnormalv <= nb );
}
}
//  static const R2 Triangle::Hat[3]= {R2(0,0),R2(1,0),R2(0,1)} ;