fem.hpp

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

template<class Rd>
class TBoundaryEdge: public Label {
public:
 typedef TVertex<Rd> Vertex;
    static const int NbWhat = 3; // 3+3+1 
    static const int NbV = 2; // 3+3+1 
    static const int NbE = 1; //
    
  Vertex *vertices[2];
  TBoundaryEdge(Vertex * v0,int i0,int i1,int r): Label(r) 
  { vertices[0]=v0+i0; vertices[1]=v0+i1; }
  void set(Vertex * v0,int i0,int i1,int r)
  { lab=r,vertices[0]=v0+i0; vertices[1]=v0+i1; }
  bool in(const Vertex * pv) const {return pv == vertices[0] || pv == vertices[1];}
  TBoundaryEdge(){}; // constructor empty for array 
  void Draw() const;
  Vertex & operator[](int i) const {return *vertices[i];}
  R length() const { return Norme2(R2(*vertices[0],*vertices[1]));}
  void Renum(Vertex   *v0, long * r) { 
    for (int i=0;i<2;i++) 
      vertices[i]=v0+r[vertices[i]-v0];}
  
  SortedTriplet what(int i,Vertex *v0,TBoundaryEdge * t0) { 
      if (i<0) ffassert(i>=0);
      else if  (i<2) return SortedTriplet(vertices[i]-v0);
      else if( (i==0) ) return SortedTriplet( vertices[0]-v0, vertices[1]-v0) ;
      else ffassert(0);
  }
  
};

typedef TBoundaryEdge<R2> BoundaryEdge;
typedef BoundaryEdge Edge;
  // typedef Tetraedre Tet;  // just to play

template<class Rd>
class TMortar { 
  public:
 typedef TVertex<Rd> Vertex;  
  friend class Mesh;
  friend class ConstructDataFElement;
   Mesh * Th;
   int nleft,nright;
   int *left,*right;
   TMortar(): Th(0),nleft(0),nright(0),left(0),right(0){}
   void Draw() const;
   public:
     int NbLeft() const{return nleft;} 
     int NbRight() const{return nright;} 
     int TLeft(int i) const { return left[i]/3;}
     //   
     int NbT() const {return nleft+nright;} 
     int T_e(int i,int & e)  const { // return the triangle number + the edge number 
        throwassert(i>=0 && i < nleft+nright);
        int k= (i<nleft ? left[i]: right[i-nleft]);
        e=k%3;
        return k/3;}
     int ELeft(int i) const { return left[i]%3;}
     int TRight(int i) const { return right[i]/3;}
     int ERight(int i) const { return right[i]%3;}
     
     //  warning  i is in [0, nleft]
     Vertex & VLeft(int i) const ;
     Vertex & VRight(int i) const;
     
     
};
typedef TMortar<R2> Mortar;

 
class FQuadTree;

class Mesh: public RefCounter { 
public:
  
  typedef TTriangle<R2> Triangle;
  typedef TTriangle<R2> Element;
  typedef BoundaryEdge  BorderElement; 
  typedef TVertex<R2>  Vertex;
  typedef  R2 Rd;
  typedef R2  RdHat;// for parametrization 
  typedef Rd::R R;
  typedef FQuadTree GTree;
  

  static const char magicmesh[8]  ;
  int dim; 
  int nt,nv,neb,ne,ntet;
  R area;
  R volume;
  R lenbord;
  static int kthrough,kfind;
  FQuadTree *quadtree; 
  Vertex *vertices;
  Triangle *triangles;
  BoundaryEdge  *bedges;
  Edge  *edges;  // edge element 
  Tetraedre * tet; // 
  
    
  int NbMortars,NbMortarsPaper;
  Mortar *mortars; //  list of mortar
  int    *datamortars;
  R2  * bnormalv; //  boundary vertex normal 
  //Triangle * adj;
  Triangle & operator[](int i) const {throwassert(i>=0 && i<nt);return triangles[i];}
 // const Triangle & operator[](int i) const {return triangles[i];}
  Vertex & operator()(int i) const {throwassert(i>=0 && i<nv);return vertices[i];}
  Mesh(const char * filename) {read(filename);} // read on a file
  Mesh(const string s) {read(s.c_str());}
  Mesh( const Serialize & ) ;
  Mesh(int nbv,R2 * P);

  R mesure(){ return area;}
  R bordermesure(){ return lenbord;}
  
  Serialize serialize() const;
  Mesh(int nbv,int nbt,int nbeb,Vertex *v,Triangle *t,BoundaryEdge  *b);  
  Mesh(const Mesh & Thold,int *split,bool WithMortar=true,int label=1);
  ~Mesh();
  int number(const Triangle & t) const {return &t - triangles;}
  int number(const Triangle * t) const {return t  - triangles;}
  int number(const Vertex & v)   const {return &v - vertices;}
  int number(const Vertex * v)   const {return v  - vertices;}
  int operator()(const Triangle & t) const {return &t - triangles;}
  int operator()(const Triangle * t) const {return t  - triangles;}
  int operator()(const Vertex & v)   const {return &v - vertices;}
  int operator()(const Vertex * v)   const {return v  - vertices;}
  int operator()(int it,int j) const {return number(triangles[it][j]);}
  BoundaryEdge &  be(int i) const { return bedges[i];}
  Element &  t(int i) const { return triangles[i];}
  Vertex &  v(int i) const { return vertices[i];}

        // Nu vertex j of triangle it
  void BoundingBox(R2 & Pmin,R2 &Pmax) const;
  void InitDraw() const ;
  void Draw(int init=2,bool fill=false) const;
  void DrawBoundary() const; 
  int Contening(const Vertex * v) const{ return TriangleConteningVertex[ v  - vertices];}
  int renum();
  int gibbsv (long* ptvoi,long* vois,long* lvois,long* w,long* v);
  int ElementAdj(int it,int &j) const 
      {int i=TheAdjacencesLink[3*it+j];j=i%3;return i/3;}
  int nTonEdge(int it,int e) const { int k=3*it+e;return k==TheAdjacencesLink[k] ? 1 : 2;}
      
  void VerticesNumberOfEdge(const Triangle & T,int j,int & j0,int & j1) const 
      {j0 =  number(T[(j+1)%3]),j1=  number(T[(j+ 2)%3]);}
  bool SensOfEdge(const Triangle & T,int j) const 
      { return  number(T[(j+1)%3]) <number(T[(j+ 2)%3]);}
      
  int GetAllElementAdj(int it,int *tabk) const
   { //  get the tab of all adj triangle to a traingle (max 3)
     //  and return the size of the tab 
      int i=0;
      tabk[i]=TheAdjacencesLink[3*it+0]/3;
      if(tabk[i] >=0 && tabk[i]!=it) i++; 
      tabk[i]=TheAdjacencesLink[3*it+1]/3;
      if(tabk[i] >=0 && tabk[i]!=it) i++; 
      tabk[i]=TheAdjacencesLink[3*it+2]/3;
      if(tabk[i] >=0 && tabk[i]!=it) i++; 
      return i;
   }
      
  int BoundaryElement(int be,int & edgeInT) const {
     int i= BoundaryEdgeHeadLink[be]; edgeInT = i%3; 
     return i/3;}

     
  Triangle * Find(const R2 & P) const ;
  const Triangle * Find(R2 P, R2 & Phat,bool & outside,const Triangle * tstart=0) const  ;
  
  BoundaryEdge * TheBoundaryEdge(int i,int j)  const
   {  
    int p2;
    for (int p=BoundaryAdjacencesHead[i];p>=0;p=BoundaryAdjacencesLink[p])
     { 
     if ( bedges[p2=p/2].in(vertices+j) )   return bedges+p2;
     }
     
    return 0;}
 void destroy() {RefCounter::destroy();}
 void MakeQuadTree() ;
      
private:
  void read(const char * filename); 
  void read(ifstream & f,bool bin=false); 
  void BuildBoundaryAdjacences();
  void ConsAdjacence();
  void Buildbnormalv(); 
  void BuilTriangles(bool empty,bool removeoutside=true);
  // to construct the adj triangle               
  int *TheAdjacencesLink;
  int *BoundaryEdgeHeadLink;
  int *BoundaryAdjacencesHead;
  int *BoundaryAdjacencesLink; 
  int *TriangleConteningVertex;       
  // no copy
  Mesh(const Mesh &);
  void operator=(const Mesh &);       
};


//  2 routines to compute the caracteristic
int WalkInTriangle(const Mesh & Th,int it, double *lambda,
                   const  KN_<R> & U,const  KN_<R> & V, R & dt);
int  WalkInTriangle(const Mesh & Th,int it, double *lambda,
                     R u, R v, R & dt); 
                   
int Walk(const Mesh & Th,int & it, R *l,
         const KN_<R>  & U,const KN_<R>  & V, R dt) ;
 
void DrawMark(R2 P,R k=0.02);

         
template<class T>
void  HeapSort(T *c,long n)
{
  long l,j,r,i;
  T crit;
  c--; // on decale de 1 pour que le tableau commence a 1
  if( n <= 1) return;
  l = n/2 + 1;
  r = n;
  while (1) { // label 2
    if(l <= 1 ) { // label 20
      crit = c[r];
      c[r--] = c[1];
    if ( r == 1 ) { c[1]=crit; return;}
    } else  crit = c[--l]; 
    j=l;
    while (1) {// label 4
      i=j;
      j=2*j;
      if  (j>r) {c[i]=crit;break;} // L8 -> G2
      if ((j<r) && (c[j] < c[j+1])) j++; // L5
      if (crit < c[j]) c[i]=c[j]; // L6+1 G4
      else {c[i]=crit;break;} //L8 -> G2
    }
  }
}

template<class T,class T1,class T2>
void  HeapSort(T *c,T1 *c1,T2 *c2,long n)
{
  long l,j,r,i;
  T crit;  T1 crit1;  T2 crit2;
  c--;c1--;c2--; // on decale de 1 pour que le tableau commence a 1
  if( n <= 1) return;
  l = n/2 + 1;
  r = n;
  while (1) { // label 2
    if(l <= 1 ) { // label 20
      crit = c[r];     crit1 = c1[r];  crit2 = c2[r];
      c2[r] = c2[1];   c1[r] = c1[1];  c[r--] = c[1];
    if ( r == 1 ) { c2[1] = crit2,   c1[1] = crit1; c[1]=crit; return;}
    } else  {crit = c[--l];crit1=c1[l];crit2=c2[l]; }
    j=l;
    while (1) {// label 4
      i=j;
      j=2*j;
      if  (j>r) {c[i]=crit;c1[i]=crit1;c2[i]=crit2;break;} // L8 -> G2
      if ((j<r) && (c[j] < c[j+1])) j++; // L5
      if (crit < c[j]) {c[i]=c[j];c1[i]=c1[j];c2[i]=c2[j];} // L6+1 G4
      else {c[i]=crit;c1[i]=crit1;c2[i]=crit2;break;} //L8 -> G2
    }
  }
}

 inline   int numSubTVertex(int N,int i,int j)
    {  //  i,j  coordonne barycentre * N dans l'ele鵨nt de reference.
	i=i+j; // numerotation / diag  
	// i,j 
	assert(j<=i && 0<= j); 
	return j+i*(i+1)/2;
    }
    
 inline    void  num1SubTVertex(int N,int l,int & i,int & j)
    {
	
      i= (int) ((-1 + sqrt(1.+8*l))/2); //  modif gcc 3.3.3  FH 100109
	j = l - i*(i+1)/2;
        // io=in+j;
	// in = io-j
	i=i-j;
	assert( l == numSubTVertex(N,i,j)); 
    }
    
  R2 SubTriangle(const int N,const int n,const int l);
    int numSubTriangle(const int N,const int n,const int l);
  int  NbOfSubTriangle(const int N);
  int NbOfSubInternalVertices(int kk);
  R2  SubInternalVertex(int N,int k);
    
      //  warning  i is in [0, nleft]
template<class Rd>
 inline     TVertex<Rd> & TMortar<Rd>::VLeft(int i) const 
  { throwassert(i>=0 && i <= nleft);
    return i< nleft ? (*Th)[TLeft( i  )][VerticesOfTriangularEdge[ELeft(i)][0]]
                    : (*Th)[TLeft( i-1)][VerticesOfTriangularEdge[ELeft(i-1)][1]];}
template<class Rd>
 inline TVertex<Rd> & TMortar<Rd>::VRight(int i) const 
  { throwassert(i>=0 && i <= nright);
    return i< nright ? (*Th)[TRight( i  )][VerticesOfTriangularEdge[ERight(i)][1]]
                     : (*Th)[TRight( i-1)][VerticesOfTriangularEdge[ERight(i-1)][0]];}

 void DrawCommentaire(const char *cm,float x=0.1,float y=0.97) ;
 
inline R2 minmax(const R2 & a,const R2 & b) 
  {return R2(Min(a.x,b.x),Max(a.y,b.y));}


}

using Fem2D::R; 

#include "FQuadTree.hpp"

namespace Fem2D {

inline void Mesh::MakeQuadTree() 
  {
    if (!quadtree) {
        R2 Pn,Px;
        BoundingBox(Pn,Px);
        quadtree=new FQuadTree(this,Pn,Px,nv);
    } 
 }

}

#endif