fem.cpp

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

// -*- Mode : c++ -*-
//
// SUMMARY  :      
// USAGE    :        
// ORG      : 
// AUTHOR   : Frederic Hecht
// E-MAIL   : [e-mail]
//

/*
 
 This file is part of Freefem++
 
 Freefem++ is free software; you can redistribute it and/or modify
 it under the terms of the GNU Lesser General Public License as published by
 the Free Software Foundation; either version 2.1 of the License, or
 (at your option) any later version.
 
 Freefem++  is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU Lesser General Public License for more details.
 
 You should have received a copy of the GNU Lesser General Public License
 along with Freefem++; if not, write to the Free Software
 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 */

extern long verbosity ;
#include <cmath>
#include  <cfloat>
#include <cstdlib>
#include "error.hpp"
#include <iostream>
#include <fstream>
//#include <strstream.h>
//using namespace std;  //introduces namespace std
#include "RNM.hpp"
#include "rgraph.hpp"
#include "Serialize.hpp"
#include "fem.hpp"
#include <set>
namespace Fem2D {
    
    
class SubMortar { 
	friend class Mesh;
	friend ostream & operator<<(ostream & f,const SubMortar & m);
	R  alpha; // angle in radian 
	R2 from,to;
	int k; //  triangle number
	int i; //  edge on triangle
	int sens; //       
		  //  int head;
public:
	    SubMortar() :alpha(0),k(0),i(0),sens(0){}
	SubMortar(const Vertex & s,const Vertex & ss,int kk,int ii,int se) 
	    : alpha(Theta((R2) ss - s)),from(s),to(ss),k(kk),i(ii),sens(se) {}
	R len2() const { return Norme2_2(to-from);}
	R len2(const R2 & A) const{ return (to-from,(A-from));}
	
	
	bool operator<(const SubMortar &  b){ return alpha < b.alpha ;}  // to sort
									 //  bool side(const Triangle &K) {}
    };
    
    ostream & operator<<(ostream & f,const SubMortar & m)
    { f << " a=" << m.alpha << " " << m.k << " " << m.i << " " << m.sens << " l^2=" <<m.len2() <<  ";" ;
	return f;}
    
    void Mesh::BuildBoundaryAdjacences()
    {
	if(!BoundaryAdjacencesHead)
	{
	    BoundaryAdjacencesHead = new int[nv];
	    BoundaryAdjacencesLink = new int[neb+neb];
	    for (int i=0;i<nv;i++) 
		BoundaryAdjacencesHead[i]=-1;
	    int j2=0;
	    for (int j=0;j<neb;j++)
		for (int k=0;k<2;k++,j2++)
		{  
		    int v = number(bedges[j][k]);
		    assert(v >=0 && v < nv);
		    BoundaryAdjacencesLink[j2]=BoundaryAdjacencesHead[v];
		    BoundaryAdjacencesHead[v]=j2;
		}
		    
	}
    }  
	void Mesh::ConsAdjacence()
    {
	    //  warning in the paper a mortar is the whole edge of the coarse triangle
	    //  here a mortar is a connected componand of he whole edge of the coarse triangle 
	    //   minus  the extremite of mortar
	    //  -----------
	    int NbCollision=0,NbOfEdges=0,NbOfBEdges=0,NbOfMEdges=0;
	    const char MaskEdge[]={1,2,4};
	    const char AddMortar[]={8,16,32};
	    //    reffecran();
	    //    cadreortho(0.22,0.22,.1);
	    if (neb) BuildBoundaryAdjacences();
	    
	    if(TheAdjacencesLink) return; //
	    TheAdjacencesLink = new int[3*nt];
	    const int NbCode = 2*nv;
	    char * TonBoundary = new char [nt]; // the edge is 1 2 4   AddMortar = 8
	    
	    { int * Head = new int[NbCode];
		
		//  make the list
		int i,j,k,n,j0,j1;
		for ( i=0; i<NbCode; i++ )
		    Head[i]=-1; // empty list
		n=0; // make all the link
		for (i=0;i<nt;i++)
		{ 
		    Triangle & T=triangles[i];
		    for( j=0; j<3; j++,n++ )
		    { 
			VerticesNumberOfEdge(T,j,j0,j1);
			k = j0+j1;
			TheAdjacencesLink[n]=Head[k]; 
			Head[k]=n; //            
		    }
		    
		}
		//
		if (neb==0) { // build boundary
		    for(int step=0;step<2;step++)
		    {
			neb=0;
			for (i=0;i<nt;i++)
			{ 
			    Triangle & T=triangles[i];
			    for( j=0; j<3; j++,n++ )
			    { 
				VerticesNumberOfEdge(T,j,j0,j1);
				int kk = 0,im=Min(j0,j1);             
				for (int n=Head[j0+j1]; n>=0; n=TheAdjacencesLink[n])
				{ int jj=n%3,ii=n/3, jj0,jj1;
				    VerticesNumberOfEdge(triangles[ii],jj,jj0,jj1);
				    if(im==Min(jj0,jj1)) // same edge 
					kk++;
				}
				if (kk==1) { 
				    if(step) bedges[neb].set(vertices,j0,j1,1);
				    neb++;
				}
				
			    }
			}
			if (step==0) {
			    if (verbosity) cout << " we build " << neb << " boundary edges" << endl;
			    bedges = new BoundaryEdge[neb];
			}
		    }
		    BuildBoundaryAdjacences(); 
		}
		for (int k=0;k<nt;k++) TonBoundary[k]=0;
		
		BoundaryEdgeHeadLink = new int[neb];       
		for (i=0;i<neb;i++)
		{  
		    BoundaryEdge & be(bedges[i]);
		    int n;
		    int i0=number(be.vertices[0]);
		    int i1=number(be.vertices[1]);
		    throwassert(i0 >=0 && i0 < nv);
		    throwassert(i1 >=0 && i1 < nv);
		    throwassert(i1 != i0) ;
		    int im=Min(i0,i1);
		    BoundaryEdgeHeadLink[i]=-1; 
		    for ( n=Head[i0+i1]; n>=0; n=TheAdjacencesLink[n])
		    {
			int jj=n%3,ii=n/3, jj0,jj1;
			VerticesNumberOfEdge(triangles[ii],jj,jj0,jj1);
			if(im==Min(jj0,jj1)) // same edge 
			{
			    TonBoundary[n/3] += MaskEdge[n%3];
			    BoundaryEdgeHeadLink[i]=n;                  
			    if(i0==jj0) break; // FH 01072005 bon cote de l'arete 
					       // sinon on regard si cela existe?
			}
		    } 
		    if ( BoundaryEdgeHeadLink[i] <0 && verbosity) 
			cout << "   Attention l'arete frontiere " << i 
			    << " n'est pas dans le maillage " <<i0 << " " << i1 <<  endl;
		}
		
		//  find adj
		// reffecran();    
		
		for (i=0;i<nt;i++)
		{ 
		    Triangle & T=triangles[i];
		    for( j=0; j<3; j++,n++ )
		    { 
			VerticesNumberOfEdge(T,j,j0,j1);
			k = j0+j1; // code of current edge 
			int jm = Min(j0,j1), NbAdj=0, He=-1;
			int *pm=Head+k;
			while (*pm>=0) // be carefull  
			{                                 
			    int m=*pm,jj=m%3,ii=m/3, jj0,jj1;
			    VerticesNumberOfEdge(triangles[ii],jj,jj0,jj1);
			    if(jm==Min(jj0,jj1)) // same edge 
			    {                  
				NbAdj++;   
				// remove from  the liste 
				*pm=TheAdjacencesLink[m];
				TheAdjacencesLink[m]=He;  // link to He
				He = m;                  
			    }
			    else 
			    {
				NbCollision++;
				pm=TheAdjacencesLink+*pm; // next 
			    }
			}
			//  make a circular link 
			if (NbAdj>0) 
			{ 
			    int m=He;
			    while(TheAdjacencesLink[m]>=0)
				m=TheAdjacencesLink[m]; // find end of list     
							// close the List of common edges               
			    TheAdjacencesLink[m] = He;               
			}
			if (NbAdj >2) 
			{
			    int m=He;
			    do { 
				m=TheAdjacencesLink[m];
			    } while(TheAdjacencesLink[m]!=He);
			}
			
			if (NbAdj) NbOfEdges++;
			if(NbAdj==1)
			    if (! (TonBoundary[i]& MaskEdge[j])) 
			    { NbOfMEdges++;
				if(verbosity>99) 
				    cout << " Edge (" << j0 << " "<< j1 << ") : "  << j  << " of Triangle " << &T-triangles << " on mortar \n"
				    <<" --- > " << number(T[0]) << " " << number(T[1]) << " " << number(T[2]) << " /" << int(TonBoundary[i])<< "\n" ;
				TonBoundary[i]+= AddMortar[j];
			    }
				else { NbOfBEdges++; }
		    }
		}
		    
		    if (verbosity>1 ) {
			cout << "    Nb of Vertices " << nv <<  " ,  Nb of Triangles " 
			<< nt << endl ;
			cout << "    Nb of edge on user boundary  " << neb 
			<< " ,  Nb of edges on true boundary  " << NbOfBEdges << endl;
			if(NbOfMEdges) cout << "    Nb of edges on Mortars  = " << NbOfMEdges << endl;

		    }
		    delete [] Head; // cleanning memory
		    NbMortars =0;
		    NbMortarsPaper=0;
	    }
		{
		    //  construct the mortar 
		    int * linkg = new int [nv]; //  to link
		    int * linkd = new int [nv]; //  to link
		    int * NextT3= new int[3*nt];
		    int * headT3= new int[nv];
		    ffassert( linkg && linkd);
		    for (int i=0;i<nv;i++) 
			headT3[i]=linkg[i]=linkd[i]=-1; // empty
							//   create the 2 link 
							//  reffecran();
							//  Draw(0);
		    for (int k=0;k<nt;k++)
			for (int j=0;j<3;j++)
			    if (TonBoundary[k] & AddMortar[j]) 
			    { 
				//   triangles[k].Draw(j,0.8);
				int s0,s1;
				VerticesNumberOfEdge(triangles[k],j,s0,s1);
				linkg[s0] = (linkg[s0] == -1) ? s1 : -2 ;
				linkd[s1] = (linkd[s1] == -1) ? s0 : -2 ; 
				//              throwassert(linkg[s0] != -1 && linkg[s1] != -1 );
				//   cout << "On Mortar " << s0 << " " << s1 << " link " <<  linkg[s0]  << " " << linkd[s1] <<endl  ;                              
			    } 
				//  we remove the  boundary link       
				for (int k=0;k<nt;k++)
				    for (int j=0;j<3;j++)
					if (TonBoundary[k] & MaskEdge[j]) 
					{ 
					    int s0,s1;
					    VerticesNumberOfEdge(triangles[k],j,s0,s1);
					    // cout << s0 << " " << s1 << " ld " << linkd[s0]  << " " << linkd[s1] << " lg " << linkg[s0]  << " " << linkg[s1] << " apres " ;              
					    linkg[s0] = linkg[s0] != -1 ?  -2 : -1;
					    linkg[s1] = linkg[s1] != -1 ?  -2 : -1;
					    
					    linkd[s1] = linkd[s1] != -1 ?  -2 : -1;                                                 
					    linkd[s0] = linkd[s0] != -1 ?  -2 : -1; 
					    // cout  << " ld " << linkd[s0]  << " " << linkd[s1] << " lg" << linkg[s0]  << " " << linkg[s1] << endl;
					    
					} 
					    
					    //    remark if   linkd[i]  == -2  extremities of mortars (more than 2 mortars)
					    //    if  ((linkd[i] != -1) || (linkd[i] != -1))   =>   i is on mortars 
					    //    make a link for all triangles contening a  mortars
					    const int k100=100;
		    SubMortar  bmortars[k100];
		    int k3j=0;
		    for (int k=0;k<nt;k++) {
			const  Triangle & T=triangles[k];
			for (int j=0;j<3;j++,k3j++) 
			{ 
			    NextT3[k3j]=-2; 
			    
			    int is= number(T[j]);
			    // if (linkg[is]<-1 || linkd[is]<-1) // extremite of bmortars
			    int j0 =EdgesVertexTriangle[j][0];  //  the 2 edges contening j 
			    int j1 =EdgesVertexTriangle[j][1];
			    if  (TonBoundary[k] & (AddMortar[j0]|AddMortar[j1]))  // (linkg[is]!=-1) 
			    { 
				throwassert(linkd[is]!=-1 || linkg[is]!=-1);
				NextT3[k3j]=headT3[is];
				headT3[is] =  k3j; 
			    }}
		    }
		    //    loop on extremite of bmortars 
		    // calcule du nombre de joint 
		    /*      rattente(1);      
		    for (int is=0;is<nv;is++)
			if ( (linkg[is]<-1 || linkd[is]<-1) )
		    {MoveTo( vertices[is]);
			    ostringstream ss;
			    ss << is ;
			    plotstring(ss.str().c_str());
		    }  */  
		    datamortars=0;
		    int kdmg = 0,kdmd=0; 
		    int kdmgaa = 0,kdmdaa=0; 
		    int step=0;
		    mortars=0;
		    NbMortars = 0;
		    NbMortarsPaper=0;
		    do { // two step  
			 // one to compute the NbMortars
			 // one to store in mortars and kdm 
			int * datag=0,*datad=0;
			ffassert(step++<2);
			if (NbMortars)
			{ // do allocation 
			    int kdm=kdmgaa+kdmdaa;
			    if (verbosity>2)
				cout << "   sizeof datamortars" << kdm << " NbMortars=" << NbMortars << endl;
			    mortars     = new Mortar [NbMortars];
			    datamortars = new int [kdm];
			    for (int i=0;i<kdm;i++)
				datamortars[i]=0;
			    datag=datamortars;
			    datad=datamortars+kdmgaa;
			    ffassert(kdmg && kdmd);
			}
			int onbm=NbMortars;
			// cout << "begin " << NbMortars << " g " << kdmg << " d " <<kdmd << endl;
			int kdmgo=kdmgaa;
			int kdmdo=kdmdaa;
			int kdm=kdmdaa+kdmgaa;
			//  reset 
			NbMortars =0;
			kdmg =0;
			kdmd =0;
			
			
			for (int is=0;is<nv;is++)
			    if (linkg[is] == -2) 
			    { // for all extremity of mortars 
				if(linkd[is] != -2)
				{
				  cout <<" Bug in mortar constrution : close to vertex "<< is << endl; 
				  ffassert(linkd[is] == -2);
				}
				const Vertex & S = vertices[is];
				R2  A(S);
				int km=0;
				int p;
				for ( p=headT3[is] ;p>=0; p=NextT3[p])
				{  //  for all nonconformitie around sm            
				    int k=p/3;
				    int i=p%3;
				    const Triangle & T(triangles[k]);
				    //throwassert( vertices + sm == &T[i]);
				    // for the 2 egdes contening the vertex j of the triangle k
				    
				    for (int jj=0;jj<2;jj++)   // 2 edge j contening i              
				    { 
					int j = EdgesVertexTriangle[i][jj];                     
					int s0,s1;
					VerticesNumberOfEdge(T,j,s0,s1);
					throwassert (s0  == is || s1 == is);
					if ( TonBoundary[k] & AddMortar[j])
					{
					    int ss,sens;
					    if ( s0 == is)
					    { ss=s1;sens=1;}
					    else
					    { ss=s0;sens=-1;}
					    const Vertex & vss( vertices[ss]);
					    bmortars[km++] = SubMortar(S,vss,k,i,sens);
					    throwassert(km<k100);
					} 
				    }
				}
				ffassert(p!=-2);
				
				// throwassert(km % 2 == 0);
				HeapSort(bmortars,km); 
				// cout <<" Nb of mortars around vertex "<< is << " = " << km<< endl;
				//  searh the same mortar (left and right)
				//  with same angle
				int im=km-1; // the last
				
				double eps = 1e-6;
				double thetai=bmortars[im].alpha-2*Pi;
				
				for (int jm=0;jm<km;im=jm++) 
				{ //  for all break of vertex 
				    
				    double theta= (bmortars[jm].alpha-thetai);
				    thetai=bmortars[jm].alpha;
				    if (theta < eps  || theta+eps > 2*Pi)
				    {//  same angle mod 2 * Pi
				     //  beginning of a mortars 
					if(mortars)  { //  set the pointeur
					    ffassert(NbMortars< onbm);
					    ffassert(datag && datad);
					    ffassert(datag< datamortars+ kdm);
					    mortars[NbMortars].left  = datag;
					    mortars[NbMortars].right = datad;
					    mortars[NbMortars].Th = this;
					}
					// cout << "   Begin mortar " << is << " " << im << " " << jm << " theta =" << thetai << endl; 
					R2 AM(A,bmortars[im].to);
					AM = AM/Norme2(AM);
					int ig(im),id(jm);
					if ( bmortars[im].sens <0) ig=jm,id=im;
					//SubMortar & mg(bmortars[ig]);
					//SubMortar & md(bmortars[id]);
					//  loop sur droite gauche
					//  meme sommet
					//   0 = gauche 1=droite 
					int sgd[2]={0,0};
					int *link[2];                   
					int nm[2]={0,0}; 
					R  ll[2]={0.0,0.0}; 
					// 
					sgd[0]=sgd[1]=is; 
					link[0] = linkg;
					link[1] = linkd;
					int gd=0; //  gd = 0 => left side  an gd=1 => right side
					
					
					int kkkk=0;
					do { //   for all 
					    
					    
					    int sm = sgd[gd]; 
					    int  dg = 1-gd;
					    // cout << " sm = " << sm << " h=" << headT3[sm] << " gb=" << gd << " autre " << sgd[dg ] << " " << kkkk << endl;
					    
					    R lAV,avam;
					    Vertex *pV=0;
					    int p,k,i,j;
					    //  search the the first start ( sens = gd )
					    throwassert(headT3[sm]>=0);// on a mortar ?? 
						for ( p=headT3[sm] ;p>=0; p=NextT3[p])
						{                
						    k=p/3;
						    i=p%3;
						    const Triangle & T(triangles[k]);
						    // couleur(1);T.Draw(0.3);
						    throwassert( vertices + sm == &T[i]);
						    // for the 2 egdes contening the vertex j of the triangle k
						    j = EdgesVertexTriangle[i][dg];
						    Vertex &V = T[VerticesOfTriangularEdge[j][dg]];
						    throwassert( &T[VerticesOfTriangularEdge[j][gd]] == vertices + sm);                
						    if ( TonBoundary[k] & AddMortar[j])
						    {  // check the sens and the direction
							
							//  cout << number(T[VerticesOfTriangularEdge[j][gd]])  << " pv=" 
							//      << number(V)  << " sm=" << sm << " " << headT3[number(V)] << endl; 
							
							R2 AV(A,V);
							lAV = Norme2(AV);
							avam = (AV,AM);
							// cout << " ---  " << avam << " > " << ll[gd] <<  " " << Abs((AM.perp(),AV) )
							//     << " " << ( avam > ll[gd] && Abs((AM.perp(),AV)) < lAV * 1e-6 ) << endl;
							// go ahead in direction AM 
							if ( (avam > ll[gd])  && Abs((AM.perp(),AV)) < lAV * 1e-6 )  
							{pV = &V;break;} //  ok good                         
						    }
						}
						if ( ! (p>=0 && pV))
						    throwassert(p>=0 && pV); //  PB reach the end without founding 
						if ( ! ( Abs((AM.perp(),A-*pV)) < 1e-5) )
						    // cout << Abs((AM.perp(),A-*pV)) <<*pV << endl, 
						    throwassert( Abs((AM.perp(),A-*pV)) < 1e-5);
						
						throwassert(sm != number(pV));