fem.cpp

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

	    kthrough++;
	    if (k++>=10000) 
	    {
		/* cout << P << endl;
		reffecran();
		Draw(0);
		triangles[its].Fill(2);
		DrawMark(P,0.01);
		rattente(1);  */
		ffassert(k++<10000);
	    }
	    int kk,n=0,nl[3];
	    
	    
	    
	    R2 & A(K[0]), & B(K[1]), & C(K[2]);
	    R l[3]={0,0,0};
	    R area2= K.area*2;
	    R eps =  -area2*1e-6;
	    l[0] = Area2(P,B,C);
	    l[1] = Area2(A,P,C);
	    l[2] = area2-l[0]-l[1];
	    if (l[0] < eps) nl[n++]=0;
	    if (l[1] < eps) nl[n++]=1;
	    if (l[2] < eps) nl[n++]=2;
	    
	    if (n==0)
	    {  // interior => return
		outside=false; 
		Phat=R2(l[1]/area2,l[2]/area2);
		return &K;
	    }
	    else if (n==1) 
		kk=0;
	    else  
		kk=BinaryRand() ? 1 : 0;
	    
	    j= nl[ kk ];
	    
	    int itt =  ElementAdj(it,j);
	    if(itt!=it && itt >=0)  
	    {
		dP=DBL_MAX;
		it=itt;
		continue;
	    }  
	    
	    //  edge j on border
	    l[j]=0;
	    
	    if ( n==2 ) 
	    {
		kk = 1-kk;
		j= nl[ kk ];
		itt =  ElementAdj(it,j);                  
		if (itt && itt != it)
		{
		    dP=DBL_MAX;
		    it=itt;
		    continue;
		}
		//  on a corner of the mesh 
		l[j]=0;
		l[3-nl[0]+nl[1]]=1;
		Phat=R2(l[1],l[2]);
		return triangles +it;
	    }
	    //   on the border 
	    //   projection Ortho
	    
	    kout++;
	    
	    int j0=(j+1)%3,i0= &K[j0]-vertices;
	    int j1=(j+2)%3,i1= &K[j1]-vertices;
	    int ii,jj,iii;
	    bool ret=false;
	    
	    R2 AB=R2(K[j0],K[j1]),  AP(K[j0],P), BP(K[j1],P);
	    R la=  (AB,AP);
	    R lb= -(AB,BP);
	    if(la<0)     
		ii= i0, jj = j0,iii=i1;
	    else if ( lb <0)
		ii= i1, jj = j1,iii=i0;
	    else // PROJECTION between A,B
		ret = true;
	    if( ! ret)
	    { //  VERIF THE DISTANCE**2 Dicrease
		R2 Pjj(P,K[jj]);
		R dd = (Pjj,Pjj);
		if (dd >= dP ) {
		    Phat=PPhat;
		    // if(kout>1) cout << "        @ " << tt-triangles  << " " << Phat << " " << outside << endl; 
		    
		    return tt;
		}
		else
		{ 
		    l[0]=l[1]=l[2]=0;
		    l[jj]=1;
		    PPhat.x=l[1];
		    PPhat.y=l[2];                
		    dP=dd;
		    tt = triangles + it ;
		}
	    }  
	    /*
	     if(kout>1)
	     cout << "Find --- "<< P << " :  k=" << k << "  la lb " << la << " " << lb 
	     << " ii =" << ii << " iii= " << iii << " " << it << " dP= " 
	     << dP << " ret = " << ret <<endl;
	     */
	    if (ret || ii == iib) 
	    {  
		l[j]=0;
		
		l[j0]= Max(0.,Min(+lb/(la+lb),1.));
		l[j1]= 1-l[j0];
		Phat=R2(l[1],l[2]);
		//if(kout>1) cout << "        # " << it << " " << Phat << " " << outside << endl; 
		return triangles +it;
	    }
	    bool ok=false;
	    // next edge on true boundary 
	    for (int p=BoundaryAdjacencesHead[ii];p>=0;p=BoundaryAdjacencesLink[p])
	    { int e=p/2, ie=p%2;// je=2-ie;
				// cout << number(bedges[e][0]) << " " << number(bedges[e][1]) << endl;
		if (!bedges[e].in( vertices+iii)) //  edge not equal  to i0 i1 
		{  
		    ok=true;
		    iib = ii;
		    it= BoundaryElement(e,ie);  //  next triangle                      
						 // cout << "  ------ " << it << " " << Phatt <<  endl;
		    break;
		}
	    }
	    ffassert(ok); 
	}
	
	
}


int  WalkInTriangle(const Mesh & Th,int it, double *lambda,
		    R u, 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]);
    
    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;
}        
Mesh::~Mesh()
{
    //(cout << "   -- delete mesh " << this << endl);
    delete  quadtree;
    delete [] triangles;
    delete [] vertices;
    delete [] bedges;
    delete [] mortars;
    delete [] datamortars;
    // delete [] adj;
    delete [] TheAdjacencesLink;
    delete [] BoundaryEdgeHeadLink;
    delete [] BoundaryAdjacencesHead;
    delete [] BoundaryAdjacencesLink;
    delete []  TriangleConteningVertex;
    delete [] bnormalv;
    delete [] tet;
    delete [] edges;
}
//  for the  mortar elements
 int NbOfSubTriangle(int k)
{  
    if(k>0) return  k*k;
    else if(k<0) return 3*(k*k);
    ffassert(0);
    return 0;
}
    

 int NbOfSubInternalVertices(int kk)
{ 
    assert(kk);
    int k=Abs(kk);
    int  r= (k+2)*(k+1)/2;
    assert(r>=0);
    
    return kk<0 ? 3*r : r;
}



Mesh::Mesh(int nbv,R2 * P)
{
    
    TheAdjacencesLink=0;
    BoundaryEdgeHeadLink=0;
    BoundaryAdjacencesHead=0;
    BoundaryAdjacencesLink=0; 
    TriangleConteningVertex=0;              
    TriangleConteningVertex=0;
    dim=2;
    tet=0;
    edges=0;
    ntet=0;
    ne=0;
    volume=0;
    
    quadtree =0;
    NbMortars=0;
    NbMortarsPaper=0;
    nt=0;
    mortars=0;
    TheAdjacencesLink =0;
    datamortars=0;
    quadtree =0;
    bedges=0;
    neb =0;
    bnormalv=0;
    nv=nbv;
    vertices  = new Vertex[nv];
    triangles=0;
    area=0;
    for (int i=0;i<nv;i++)
	vertices[i]=P[i];
    bedges    =  0; 
    area=0;
    
    
    BuilTriangles(false) ;  
    
    ConsAdjacence();    
    
}
void Mesh::BuilTriangles(bool empty,bool removeouside)
{
    long nba = neb;
    // 
    long nbsd = 0; // bofbof 
    if(!removeouside) nbsd=1;
    // faux just pour un test 
    long  *sd;
    sd=new long[2];
    sd[0]=-1;
    sd[1]=0;        
    nbsd=removeouside?0:-1;
    
    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  = nba ? new long[2*nba] : 0; 
	nu = new long[6*nbtmax];
	c = new long[2*nbsmax];
	tri = new long[(4 * nbsmax + 2 * nbsd)];
	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; 
	}
	for (int i=0,k=0 ;i<neb;i++)
	{
	    arete[k++] =number(bedges[i][0])+1;
	    arete[k++] =number(bedges[i][1])+1;
	}
	
	
	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;
	mshptg8_ (cr, h, c, nu, &nbs, nbs, tri, arete, nba, (long *) sd, nbsd, reft, &nbt, .25, .75, &err);
	if(err) {
	    cerr << " Sorry Error build delaunay triangle   error = " << err << endl;
	    delete [] arete;
	    delete [] nu;
	    delete [] c;
	    delete [] tri;
	    delete [] reft;
	    delete [] cr;
	    delete [] h;
	    delete [] sd;	   
	    throw(ErrorExec("Error mshptg8_",1));
	   	}
	assert(err==0 && nbt !=0);
	delete [] triangles;
	nt = nbt;
	if(verbosity>1)
	    
	    cout << " Nb Triangles = " << nbt << endl;
	triangles = new Triangle[nt];
	for(int i=0,k=0;i<nt;i++,k+=3)
        {
	    triangles[i].set(vertices,nu[k]-1,nu[k+1]-1,nu[k+2]-1,reft[i]);
	    area += triangles[i].area;
        }  
	
	delete [] arete;
	delete [] nu;
	delete [] c;
	delete [] tri;
	delete [] reft;
	delete [] cr;
	delete [] h;
	delete [] sd;
}
Mesh::Mesh(const Mesh & Th,int * split,bool WithMortar,int label)
{ //  routine complique 
  //  count the number of elements
    area=Th.area;
    volume=0;
    BoundaryAdjacencesHead=0;
    BoundaryAdjacencesLink=0;
    BoundaryEdgeHeadLink=0;
    quadtree =0;
    NbMortars=0;
    ntet=0;
    ne=0;
    dim=2;
    tet=0;
    edges=0;
    mortars=0;
    TheAdjacencesLink =0;
    quadtree =0;
    bedges=0;
    neb =0;
    bnormalv=0;
    R2 Pmin,Pmax;
    Th.BoundingBox(Pmin,Pmax);
    nt=0;
    int nebi=0; // nb arete interne 
    int nbsdd =0;
    int splitmin=100,splitmax=0;
    for (int i=0;i<Th.nt;i++)
	if(split[i]) 
	{
	    splitmin=Min(splitmin, split[i]);
	    splitmax=Max(splitmax, split[i]);	    
	    nt += NbOfSubTriangle(split[i]);
	}
    
    bool constsplit=splitmin==splitmax;
    bool noregenereration = constsplit || WithMortar;

    if(verbosity>2) 
	cout << "  -  Mesh construct : from " << &Th << " split min " << splitmin 
	    << "  max " << splitmax << ", recreate " << !noregenereration 
	    << " label =" << label << endl;
    
    triangles = new Triangle[nt];
    assert(triangles);
    //  computation of thee numbers of vertices
    //  on decoupe tous betement
    // et on recolle le points ensuite 
    // avec le quadtree qui sont sur les aretes ou les sommets 
    //  calcul du magorant du nombre de sommets 
    int nvmax= 0;// Th.nv;
    { //int v;
	KN<bool> setofv(Th.nv,false);
	for (int i=0;i<Th.nt;i++)
	    if ( split[i] ) 
	    {   nbsdd++;
		for (int j=0;j<3;j++)
		{
		    int jt=j,it=Th.ElementAdj(i,jt);
		    if(it==i || it <0) neb += split[i];  //on est sur la frontiere
		    else if  (!split[it]) neb += split[i];//le voisin ne doit pas etre decoupe
		    else  //on est dans le domaine et le voisin doit etre decoupe
		    {
			int ie0,ie1;
			Th.VerticesNumberOfEdge(Th[i],j,ie0,ie1);
			BoundaryEdge * pbe = Th.TheBoundaryEdge(ie0,ie1);
			if(pbe && &(*pbe)[0] == &Th(ie0))   
			    neb += max(split[i],split[it]); // aretes frontiere (FH juillet 2005)
			if (!pbe && (ie0 < ie1))
			{
			    nebi += max(split[i],split[it]); // arete interne a force ...  (FH jan 2007) 
			   // cout << nebi << " zzzz  t" << i << " a=" << j << " taj=" << it << " sa =" <<    ie0 << " " << ie1 << " + = " << max(split[i],split[it]) << endl;
			}
		    }
		}
		for (int j=0;j<3;j++)
		    //   if ( setofv.insert(Th(i,j) ).second ) 
		    if ( !setofv[Th(i,j)]) 
		    { 
			setofv[Th(i,j)]=true;
			//  cout << Th(i,j)  << " " ;
			nvmax++; 
		    }
	    }
		if(verbosity>4)
		    cout << "  - nv old " << nvmax << endl;
    }
		
	int nebmax=neb;
	int nebimax=nebi;   
	int nbsddmax=nbsdd;
	for (int i=0;i<Th.nt;i++)
	    if(split[i])
		nvmax += NbOfSubInternalVertices(split[i]) -3;
	
	//  compute the minimal Hsize of the new mesh  
	R hm = Th[0].h_min();
	//  cout << " hm " << hm << endl;
	// change h() in h_min() bug correct  july 2005 FH  ----
	for (int it=0;it<Th.nt;it++)
	{ 
	    assert(split[it]>=0 && split[it]<=64);
	    //      cout << " it = " <<it << " h " <<  Th[it].h() << " " << split[it] << " hm = " << hm <<  endl;
	    if (split[it]) 
		hm=Min(hm,Th[it].h_min()/(R) split[it]);
	}
	R seuil=hm/400.0;
	if(verbosity>5)   
	    cout << " seuil = " <<  seuil << " hmin = " << hm <<  endl; 
	assert(seuil>1e-15);
	vertices = new Vertex[nvmax];
	assert( vertices );