fem.cpp

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

						int kkgd= 3*k + j;   
						ll[gd] = avam;
						//   find the SubMortar m of vertex  sm
						//   with same sens and direction 
						//   
						
						for (int s= number(pV);s>=0;s=link[gd][s],nm[gd]++) 
						{
						    //  on est sur l'arete kkgd
						    
						    throwassert( s == number(triangles[kkgd/3][VerticesOfTriangularEdge[kkgd%3][dg]]));
						    sgd[gd]=s;// save the last
							if ( ! ( Abs((AM.perp(),A-vertices[s])) < 1e-5) )
							    // cout << Abs((AM.perp(),A-vertices[s])) << vertices[s] << endl, 
							    throwassert( Abs((AM.perp(),A-vertices[s])) < 1e-5);
							//cout << " s=" << s << " h=" << headT3[s] << " " << link[gd][s] << " " <<  link[dg][s] << endl; ;
							throwassert(kkgd>=0 && kkgd < 3*nt);
							if (datamortars) 
							{
							    throwassert(datag - datamortars == nm[0] + kdmg);
							    throwassert(datad - datamortars == nm[1] + kdmd + kdmgo );
							    
							    if (gd == 0)  *datag++ = kkgd; // store 
							    else          *datad++ = kkgd; //
							    
							}  
							
							//                            cout  << " ++++ "<<  ll[gd] << " > " << ll[dg]  << " " << headT3[sgd[dg]] << " " <<sgd[dg] << endl;    
							
							int kk=kkgd,kkk=0,kkkk=0;
							if ( link[gd][s] >=0) { 
							    for (int pp=headT3[s] ;pp>=0; pp=NextT3[pp],kkk++)
							    {  throwassert(number(triangles[pp/3][pp%3]) == s);
								
								if( (pp/3)!= (kk/3))
								{
								    kkkk++,kkgd=pp - (pp%3) + EdgesVertexTriangle[pp%3][dg];                                  
								}
							    }
							    throwassert( kkgd/3 != kk/3);
							    throwassert(kkk==2 && kkkk==1);}  
						} 
						
						if (ll[gd]>ll[dg] &&  headT3[sgd[dg]]>=0) //changement de cote 
						    gd = dg; 
						throwassert(kkkk++<100);
						//cout <<" kkkk=" << kkkk << " " << sgd[0] << " " << sgd[1] << endl;
					} while (sgd[0] != sgd[1]); 
					
					kdmgaa = Max(kdmgaa,kdmg  + nm[0]);
					kdmdaa = Max(kdmdaa,kdmd  + nm[1]);                                                                    
					
					if (is < sgd[0]  &&  headT3[sgd[0]] >=0) {
					    //cout << "    Mortars from (on saute) " << is << " to " << sgd[0] << " " << nm[0] << " " << nm[1]<< " " <<  kdmg <<  " " << kdmd << endl;
					    if( mortars ) { //  restore 
						datag -= nm[0];
						datad -= nm[1];   } 
					    
					} 
					else {
					    // cout << "    Mortars from " << is << " to " << sgd[0] << " " << nm[0] << " " << nm[1]<< " " <<  kdmg <<  " " << kdmd << endl;
					    
					    if(mortars ) {
						ffassert(NbMortars< onbm);
						mortars[NbMortars].nleft  = nm[0];
						mortars[NbMortars].nright = nm[1];
						
						//  check
						for (int i=0;i<mortars[NbMortars].nleft;++i)
						    if ( mortars[NbMortars].left[i] <0 ||  mortars[NbMortars].left[i] < 3*nt)
							ffassert(0);
						for (int i=0;i<mortars[NbMortars].nright;++i)
						    if ( mortars[NbMortars].right[i] <0 ||  mortars[NbMortars].right[i] < 3*nt)
							ffassert(0);                             
						ffassert(datag <= datamortars + kdmgo + kdmdo); 
						ffassert(datad <= datamortars + kdmgo + kdmdo); 
					    }
					    kdmg += nm[0];
					    kdmd += nm[1]; 
					    NbMortars++;
					}
				    } // same angle 
				}//  for all break of vertex                  
			    } // for all extremity of mortars 
				if (verbosity>1 && NbMortars) 
				    cout << "    Nb Mortars " << NbMortars << /*" " << kdmg << " "<<  kdmd <<*/ endl;
			if (mortars) 
			{
			    // cout << "end " << NbMortars << " g " << kdmg << " d " <<kdmd << endl;
			    // cout << kdmgo << " " << kdmg << " " << kdmdo << " " << kdmd << endl;
			    ffassert(kdmgo == kdmgaa && kdmdo == kdmdaa);}
		    }  while (NbMortars && !mortars) ;
		    
		    //   rattente(1);
		    //   reffecran();
		    //    compute the NbMortarsPaper
		    int t3,nt3 = nt*3;
		    NbMortarsPaper=0;
		    for (int i=0;i<nt3;i++) 
			NextT3[i]=0;
		    for (int i=0;i<NbMortars;i++)
		    {
			
			if (mortars[i].nleft==1)
			{ t3=mortars[i].left[0];}
			else  if (mortars[i].nright==1)
			{ t3=mortars[i].right[0];}
			else { cout << "   -- Bizarre " << mortars[i].nleft << " " << mortars[i].nright << endl;
			}
			if (NextT3[t3]==0) NbMortarsPaper++;
			NextT3[t3]++;
		    }
		    delete [] linkg;
		    delete [] linkd;
		    delete [] NextT3;
		    delete [] headT3;
		    
		}//  end of mortar construction   
		
		
		
		delete [] TonBoundary;
		
		//  construction of TriangleConteningVertex
		TriangleConteningVertex = new int[nv];
		for (int it=0;it<nt;it++)
		    for (int j=0;j<3;j++)
			TriangleConteningVertex[(*this)(it,j)]=it;
		
		Buildbnormalv();
		if (verbosity>4) 
		{ 
		    cout << "    Number of Edges                 = " << NbOfEdges << endl;
		    cout << "    Number of Boundary Edges        = " << NbOfBEdges << " neb = " << neb << endl;
		    cout << "    Number of Mortars  Edges        = " << NbOfMEdges << endl;
		    cout << "    Nb Of Mortars with Paper Def    = " <<  NbMortarsPaper << " Nb Of Mortars = " << NbMortars;             
		    cout << "    Euler Number nt- NbOfEdges + nv = " 
			<< nt + NbMortars - NbOfEdges + nv << "= Nb of Connected Componant - Nb Of Hole " 
			<<endl;}
		
    }
	    void  Mesh::BoundingBox(R2 &Pmin,R2 &Pmax) const 
	    {
		ffassert(nv);
		Pmin=Pmax=vertices[0];
		for (int i=0;i<nv;i++)
		{ 
		    const R2 & P=vertices[i];
		    Pmin.x = Min(Pmin.x,P.x);
		    Pmin.y = Min(Pmin.y,P.y);
		    Pmax.x = Max(Pmax.x,P.x);
		    Pmax.y = Max(Pmax.y,P.y);
		} 
		//  cout << " Bounding Box = " << Pmin <<" " <<  Pmax << endl;       
	    }
    void Mesh::read(const char * filename)
    {
	ifstream f(filename);
	if (!f) {
	    cerr << "Erreur ouverture du fichier " << filename << endl;
	throw(ErrorExec("exit",1));}
	// ffassert(f);
	if(verbosity)
	    cout << " -- Mesh::read On file \"" <<filename<<"\""<<  endl;
	read(f);
    }
    void Mesh::read(ifstream & f , bool bin )
    { // read the mesh
	dim=2;
	ne=0;
	ntet=0;
	volume=0;
	TriangleConteningVertex =0;
	BoundaryAdjacencesHead=0;
	BoundaryAdjacencesLink=0;
	BoundaryEdgeHeadLink=0;
	quadtree =0;
	NbMortars=0;
	tet=0;
	edges=0;    
	mortars=0;
	TheAdjacencesLink =0;
	area=0;
	bnormalv=0;
	
	int i,i0,i1,i2,ir;
	
	
	f >> nv >> nt >> neb ;
	if(verbosity>10)
	    cout << "    Nb of Vertex " << nv << " " << " Nb of Triangles " 
	    << nt << " Nb of boundary edge " << neb <<  endl;
	ffassert(f.good() && nt && nv) ;
	triangles = new Triangle [nt];
	vertices  = new Vertex[nv];
	bedges    = new BoundaryEdge[neb];
	area=0;
	ffassert(triangles && vertices && bedges);
	
	for (i=0;i<nv;i++)    
	    f >> vertices[i],ffassert(f.good());
	
	for (i=0;i<nt;i++) { 
	    f >> i0 >> i1 >> i2 >> ir;
	    ffassert(f.good() && i0>0 && i0<=nv && i1>0 && i1<=nv && i2>0 && i2<=nv);
	    triangles[i].set(vertices,i0-1,i1-1,i2-1,ir); 
	area += triangles[i].area;}
	
	for (i=0;i<neb;i++) { 
	    f >> i0 >> i1 >> ir;
	bedges[i] = BoundaryEdge(vertices,i0-1,i1-1,ir); }
	
	if(verbosity)
	    cout << "   End of read: area on mesh = " << area <<endl;  
	ConsAdjacence();
	//   BoundingBox(cMin,cMax);//  Set of cMin,Cmax
    }
    
	    Mesh::Mesh(int nbv,int nbt,int nbeb,Vertex *v,Triangle *t,BoundaryEdge  *b)
	    {
		TriangleConteningVertex =0;
		BoundaryAdjacencesHead=0;
		BoundaryAdjacencesLink=0;
		BoundaryEdgeHeadLink=0;
		quadtree =0;
		NbMortars=0;
		mortars=0;
		dim=2;
		tet=0;
		volume=0;
		edges=0;
		ntet=0;
		ne=0;    
		TheAdjacencesLink =0;
		nv=nbv;
		nt =nbt;
		neb=nbeb;
		triangles=t;
		vertices=v;
		bedges=b;
		area=0;
		bnormalv=0;
		if (t && nt >0) 
		{
		    for (int i=0;i<nt;i++)  
			area += triangles[i].area;
		    ConsAdjacence();  
		}  
		else
		{
		    bool removeouside=nbt>=0;
		    nt=0;
		    BuilTriangles(true,removeouside);
		    /*    if( ! removeouside)
		    {
			neb=0; // remove the edge
			delete [] bedges;
			bedges=0;        
			
		    } */
		    ConsAdjacence();
		}
		// cout << " build: Mesh : " << this << endl;
	    }  
	    
	    inline int BinaryRand(){
#ifdef RAND_MAX
		const long HalfRandMax = RAND_MAX/2;
		return rand() <HalfRandMax;
#else
		return rand() & 16384; // 2^14 (
#endif
		
} 
int  WalkInTriangle(const Mesh & Th,int it, double *lambda,
                    const  KN_<R> & U,const  KN_<R> & V, R & dt)
{
    const Triangle & T(Th[it]);
    const R2 Q[3]={(const R2) T[0],(const R2) T[1],(const R2) T[2]};
    int i0=Th.number(T[0]);
    int i1=Th.number(T[1]);
    int i2=Th.number(T[2]);
    
    R u   = lambda[0]*U[i0] + lambda[1]*U[i1] + lambda[2]*U[i2];
    R v   = lambda[0]*V[i0] + lambda[1]*V[i1] + lambda[2]*V[i2];
    R2 P  = lambda[0]*Q[0]  + lambda[1]*Q[1]  + lambda[2]*Q[2];
    
    //  cout << " " << u << " " << v ;
    R2 PF = P + R2(u,v)*dt;
    
    //  couleur(15);MoveTo( P); LineTo( PF);
    R l[3];
    l[0] = Area2(PF  ,Q[1],Q[2]);
    l[1] = Area2(Q[0],PF  ,Q[2]);
    l[2] = Area2(Q[0],Q[1],PF  );
    R Det = l[0]+l[1]+l[2];
    l[0] /= Det;
    l[1] /= Det;
    l[2] /= Det;
    const R eps = 1e-5;
    int neg[3],k=0;
    int kk=-1;
    if (l[0]>-eps && l[1]>-eps && l[2]>-eps) 
    {
	dt =0;
	lambda[0] = l[0];
	lambda[1] = l[1];
	lambda[2] = l[2];
    }
    else 
    {
	
	if (l[0]<eps && lambda[0] != l[0]) neg[k++]=0;
	if (l[1]<eps && lambda[1] != l[1]) neg[k++]=1;
	if (l[2]<eps && lambda[2] != l[2]) neg[k++]=2;
	R eps1 = T.area     * 1.e-5;
	
	if (k==2) // 2  
	{
	    // let j be the vertex beetween the 2 edges 
	    int j = 3-neg[0]-neg[1];
	    // 
	    R S = Area2(P,PF,Q[j]);
	    
	    if (S>eps1)
		kk = (j+1)%3;
	    else if (S<-eps1)
		kk = (j+2)%3;
	    else if (BinaryRand())
		kk = (j+1)%3;
	    else
		kk = (j+2)%3;
	    
	} 
	else if (k==1)
	    kk = neg[0];
	if(kk>=0)
	{
	    R d=lambda[kk]-l[kk];
	    
	    throwassert(d);
	    R coef =  lambda[kk]/d;
	    R coef1 = 1-coef;
	    dt        = dt*coef1;
	    lambda[0] = lambda[0]*coef1 + coef *l[0];
	    lambda[1] = lambda[1]*coef1 + coef *l[1];
	    lambda[2] = lambda[2]*coef1 + coef *l[2];
	    lambda[kk] =0;
	}
    }
    int jj=0;
    R lmx=lambda[0];
    if (lmx<lambda[1])  jj=1, lmx=lambda[1];
    if (lmx<lambda[2])  jj=2, lmx=lambda[2];
    if(lambda[0]<0) lambda[jj] += lambda[0],lambda[0]=0;
    if(lambda[1]<0) lambda[jj] += lambda[1],lambda[1]=0;
    if(lambda[2]<0) lambda[jj] += lambda[2],lambda[2]=0;
    return kk;
}
R2  SubInternalVertex(int N,int k)
    {
	if(N<0)
	  {
	      R eps = 1e-08;
	      R2 p[3][3]= { 
		  { R2(1-eps,+eps),R2(eps,1-eps),R2(1./3.+eps,1./3.+eps) },
		  { R2(0,1-eps)  ,R2(0,+eps),R2(1./3.-eps,1./3.) },
		  { R2(eps,0),R2(1-eps,0),R2(1./3.,1./3.-eps) } };
	      
	      int j=k%3;
	      R2 P=SubInternalVertex(-N,k/3);
	      R l0=1.-P.x-P.y,l1=P.x,l2=P.y;
	      return p[j][0]*l0+ p[j][1]*l1+ p[j][2]*l2;
	  }
	int i,j;
	num1SubTVertex(N,k,i,j);
	return R2( (double) i/ (double)N,(double) j/(double)N);
    }
    
R2 SubTriangle(const int N,const int n,const int l)
{
    // compute the subdivision of a triangle in N*N
    // N number of sub division
    // n number of the sub triangle
    // l vertex of the sub triangle
    if(N<0)
    {
	R eps = 1e-08;
	R2 p[3][3]= { 
	    { R2(1-eps,+eps),R2(eps,1-eps),R2(1./3.+eps,1./3.+eps) },
	    { R2(0,1-eps)  ,R2(0,+eps),R2(1./3.-eps,1./3.) },
	    { R2(eps,0),R2(1-eps,0),R2(1./3.,1./3.-eps) } };
	    
	int j=n%3;
	R2 P=SubTriangle(-N,n/3,l);
	R l0=1.-P.x-P.y,l1=P.x,l2=P.y;
	return p[j][0]*l0+ p[j][1]*l1+ p[j][2]*l2;
    }
    throwassert(n < N*N);
    int i = n % N;
    int j = n / N;
    int k = N - i - j;
    if(k<=0)
      {
	if(l==1) i++;
	else if(l==2) j++;
      }
    else
      if(l==1) j++;
      else if(l==2) i++;
    // if(l==1) i++;
    // if(l==2) j++;
    // if ( k <= 0 )cout << " - " << endl;
    return k >0 
	? R2( (double) i/ (double)N,(double) j/(double)N)
	: R2( (double) (N-j)/ (double)N,(double) (N-i)/(double)N);
    
} 

int numSubTriangle(const int N,const int n,const int l)
    {
	// compute the subdivision of a triangle in N*N
	// N number of sub division
	// n number of the sub triangle
	// l vertex of the sub triangle
	if(N<0)
	  {
	    int j=n%3;
	    return numSubTriangle(-N,n/3,l)*3+j;
	  }
	throwassert(n < N*N);
	int i = n % N;
	int j = n / N;
	int k = N - i - j;
	if(k<=0)
	  {
	    if(l==1) i++;
	    else if(l==2) j++;
	  }
	else
	  if(l==1) j++;
	  else if(l==2) i++;

	// if ( k <= 0 )cout << " - " << endl;
	return k >0 
	? numSubTVertex(N,i, j)
	: numSubTVertex(N,N-j,N-i);
	
} 
    

int Walk(const Mesh & Th,int& it, R *l,
         const KN_<R>  & U,const KN_<R>  & V, R dt) 
{
    
    int k=0;
    int j; 
    while ( (j=WalkInTriangle(Th,it,l,U,V,dt))>=0) 
    { 
	//int jj  = j;
	throwassert( l[j] == 0);
	R a= l[(j+1)%3], b= l[(j+2)%3];
	int itt =  Th.ElementAdj(it,j);
	if(itt==it || itt <0)  return -1;
	it = itt;
	l[j]=0;
	l[(j+1)%3] = b;
	l[(j+2)%3] = a;
	ffassert(k++<1000);
    }
    return it;
}

const Triangle *  Mesh::Find( R2 P, R2 & Phat,bool & outside,const Triangle * tstart) const
{
    int it,j;
    if ( tstart )
	it =  (*this)(tstart);
    else  
    {  
	const Vertex * v=quadtree->NearestVertexWithNormal(P);
	if (!v) 
	{ 
	    v=quadtree->NearestVertex(P);
	    assert(v);
	}
	/*
	 if (verbosity>100) 
	 cout << endl << (*this)(v) << *v << " " << Norme2(P-*v) << endl; 
	 */
	it=Contening(v);
    }
    
    //     int itdeb=it;     
    //     int count=0;
    //     L1: 
    outside=true; 
    //int its=it;
    int iib=-1;//,iit=-1;
	R dP=DBL_MAX;
	R2 PPhat;
	const Triangle * tt;
	int k=0,kout=0;    
	kfind++;
	while (1)
	{ 
	    const Triangle & K(triangles[it]);