fespace.hpp

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

	Sub_ToFE[kk++]=t[i]->Sub_ToFE[j];
    assert(begin_dfcomp[0]==0 && end_dfcomp[N-1]==NbDoF);
       
    Sub_ToFE= this;
    throwassert(dim_which_sub_fem[N-1]>=0 && dim_which_sub_fem[N-1]< nb_sub_fem);} 

  TypeOfFE(const int i,const int j,const int k,const int NN,const  int  *   data,int nsub,int nbsubfem,
    int kPi,int npPi,double * coef_Pi_h_a=0) 
   : NbDoF(3*(i+j)+k),
     NbNodeOnVertex(NbNodebyWhat(data,NbDoF,0)),
     NbNodeOnEdge(NbNodebyWhat(data,NbDoF,3)),
     NbNodeOnElement(NbNodebyWhat(data,NbDoF,6)),     
/*     NbDfOnVertex(Count(data,NbDoF,0)),
     NbDfOnEdge(Count(data,NbDoF,3)),
     NbDfOnElement(Count(data,NbDoF,6)),
 */  
     NbDfOnVertex(i),NbDfOnEdge(j),NbDfOnElement(k),N(NN),nb_sub_fem(nbsubfem),
     NbNode( (NbDfOnVertex ? 3 :0) + (NbDfOnEdge ? 3 :0 ) +(NbDfOnElement? 1 :0) ),
    nbsubdivision(nsub),
    DFOnWhat(data),
    DFOfNode(data+NbDoF),
    NodeOfDF(data+2*NbDoF),
    fromFE(data+3*NbDoF),
    fromDF(data+4*NbDoF),
    dim_which_sub_fem(data+5*NbDoF),
    pij_alpha(kPi),P_Pi_h(npPi),
    coef_Pi_h_alpha(coef_Pi_h_a),
    Sub_ToFE(nb_sub_fem),
     fromASubFE(data+3*NbDoF),
     fromASubDF(data+4*NbDoF),
     begin_dfcomp(data+5*NbDoF+N),
     end_dfcomp(data+5*NbDoF+2*N)    
     { 
      Sub_ToFE= this;
      assert(begin_dfcomp[0]==0 && end_dfcomp[N-1]==NbDoF);
     // cout << "TypeOfFE " <<NbDoF << " : " << NbDfOnVertex << " " << NbDfOnEdge << " " << NbDfOnElement << 
     // " : " << NbNodeOnVertex << " " << NbNodeOnEdge << " " << NbNodeOnElement << endl;
      assert(NbDfOnVertex==Count(data,NbDoF,0));
      assert(NbDfOnVertex==Count(data,NbDoF,1));
      assert(NbDfOnVertex==Count(data,NbDoF,2));
      assert(NbDfOnEdge==Count(data,NbDoF,3));
      assert(NbDfOnEdge==Count(data,NbDoF,4));
      assert(NbDfOnEdge==Count(data,NbDoF,5));
      assert(NbDfOnElement==Count(data,NbDoF,6));
      
      throwassert(dim_which_sub_fem[N-1]>=0 && dim_which_sub_fem[N-1]< nb_sub_fem);} 
      
  TypeOfFE(const int nbdf,const int NN,const  int  *   data,int nsub,int nbsubfem,
    int kPi,int npPi,double * coef_Pi_h_a=0) 
   : 
     NbDoF(nbdf),
     NbNodeOnVertex(NbNodebyWhat(data,NbDoF,0)),
     NbNodeOnEdge(NbNodebyWhat(data,NbDoF,3)),
     NbNodeOnElement(NbNodebyWhat(data,NbDoF,6)),     
     NbDfOnVertex(Count(data,NbDoF,0)),
     NbDfOnEdge(Count(data,NbDoF,3)),
     NbDfOnElement(Count(data,NbDoF,6)),
     N(NN),
     nb_sub_fem(nbsubfem),
     NbNode( (NbDfOnVertex ? 3 :0) + (NbDfOnEdge ? 3 :0 ) +(NbDfOnElement? 1 :0) ),
    nbsubdivision(nsub),
    DFOnWhat(data),
    DFOfNode(data+NbDoF),
    NodeOfDF(data+2*NbDoF),
    fromFE(data+3*NbDoF),
    fromDF(data+4*NbDoF),
    dim_which_sub_fem(data+5*NbDoF),
    pij_alpha(kPi),P_Pi_h(npPi),
    coef_Pi_h_alpha(coef_Pi_h_a),
    Sub_ToFE(nb_sub_fem) ,   
    fromASubFE(data+3*NbDoF),
     fromASubDF(data+4*NbDoF),
     begin_dfcomp(data+5*NbDoF+N),
     end_dfcomp(data+5*NbDoF+2*N)    

     { 
      Sub_ToFE= this;
      assert(begin_dfcomp[0]==0 && end_dfcomp[N-1]==NbDoF);
      assert(NbDfOnVertex==Count(data,NbDoF,0));
      assert(NbDfOnVertex==Count(data,NbDoF,1));
      assert(NbDfOnVertex==Count(data,NbDoF,2));
      assert(NbDfOnEdge==Count(data,NbDoF,3));
      assert(NbDfOnEdge==Count(data,NbDoF,4));
      assert(NbDfOnEdge==Count(data,NbDoF,5));
      assert(NbDfOnElement==Count(data,NbDoF,6));
      
      throwassert(dim_which_sub_fem[N-1]>=0 && dim_which_sub_fem[N-1]< nb_sub_fem);} 

  virtual ~TypeOfFE() { }
  virtual R operator()(const FElement & K,const  R2 & PHat,const KN_<R> & u,int componante,int op)
 const ;
 
 private:
 static int Count(const int *data,int n,int which) 
   {
      int kk=0;
      for (int i=0;i<n;++i)
        if (which == data[i]) ++kk;
      return kk;}
      
 static int LastNode(const int *data,int n) 
   {
      int kk=0,i0=2*n;
      for(int i=0;i<n;i++)
        kk=Max(kk,data[i+i0]);
      return kk;}
      
static int NbNodebyWhat(const int *data,int n,int on) 
  { 
      const int nmax=100;
      int t[nmax];
      for (int i=0;i<n;i++)
      	t[i]=0;
       int k=0,i0=2*n;
       for(int i=0;i<n;i++)
          if ( data[i]==on)
           { int node= data[i+i0];
           //  cout << " node " << node << endl;
             assert(node < nmax);
            if( ! t[node] )
               t[node] = 1,++k;
            } 
           
//     cout << " on " << on << " k = " << k << endl;
     return k;
  }      
      
} ; 



class aSubFMortar;

class TypeOfMortar { 
  friend class FESpace;
  friend class FMortar;
  friend class ConstructDataFElement;
  virtual int NbLagrangeMult(const Mesh &Th,const Mortar &M) const =0;
  virtual int NbDoF(const Mesh &Th,const Mortar &M,int i) const =0;
  virtual int NbOfNodes(const Mesh &Th,const Mortar &M) const =0;
  virtual int NbDoF(const Mesh &Th,const Mortar &M) const =0;
  virtual int NodeOfDF(const FESpace &,const Mortar &M,int i) const =0;
  virtual int DFOfNode(const FESpace &,const Mortar &M,int i) const =0;
  virtual void ConstructionOfNode(const Mesh &Th,int im,int * NodesOfElement,int *FirstNodeOfElement,int &lastnodenumber) const=0; 
  virtual void ConsTheSubMortar(FMortar &) const =0; 
  protected:
  const int NbDfOnVertex, NbDfOnEdge , N;
  public: 
    TypeOfMortar (int i,int j): NbDfOnVertex(i),NbDfOnEdge(j),N(1){};
   
  virtual  ~TypeOfMortar(){}
     
} ;

class FElement; 

  class baseFElement { public:
      typedef Fem2D::FESpace FESpace;
    typedef Fem2D::Mesh Mesh;
    typedef Fem2D::Triangle Element;
    typedef Fem2D::TypeOfFE  TypeOfFE;

  const FESpace &Vh;  
  const Triangle &T;
  const TypeOfFE * tfe; 
  const int N;
  const int number;  
  baseFElement(  const FESpace &aVh, int k) ;
  baseFElement(const   baseFElement & K,  const TypeOfFE & atfe) ;
  R EdgeOrientation(int i) const {return T.EdgeOrientation(i);}
  //  : Vh(aVh),T(Vh.Th[k]),tfe(aVh.TFE[k]),N(aVh.N),number(k){}

};

class FElement : public baseFElement { public:
  typedef const KN<TypeOfFE::IPJ> &  aIPJ;
  typedef  TypeOfFE::IPJ   IPJ;
  typedef const KN<R2> &  aR2;
  
  
  friend class Fem2D::FESpace;
  const int *p;
  const int nb;
  FElement(const FESpace * VVh,int k) ;
  public:
  int NbOfNodes()const {return nb;}
  int  operator[](int i) const ;//{ return  p ? p+i :  ((&T[i])-Vh.Th.vertices);}  N\u00c9 du noeud 
  int  NbDoF(int i) const ;  // number of DF 
  int  operator()(int i,int df) const ;// { N\u00c9 du DoF du noeud i de df local df 
  int  operator()(int df) const { return operator()(NodeOfDF(df),DFOfNode(df));}
  void Draw(const KN_<R> & U, const  KN_<R> & VIso,int j=0) const ;  
    
  void SaveDraw(const KN_<R> & U,int j,R* Usave) const ;  
  void SaveDraw(const KN_<R> & U,const KN_<R> & V,int iU,int iV,R* UVsave) const ;  

  void Drawfill(const KN_<R> & U, const  KN_<R> & VIso,int j=0,double rapz=1) const ;  
  void Draw(const RN_& U,const RN_& V, const  KN_<R> & Viso,R coef,int i0,int i1) const;
  R2   MinMax(const RN_& U,const RN_& V,int i0,int i1) const  ;
  R2   MinMax(const RN_& U,int i0) const  ;
  void BF(const R2 & P,RNMK_ & val) const;// { tfe->FB(Vh.Th,T,P,val);}
  void BF(const bool * whatd, const R2 & P,RNMK_ & val) const;// { tfe->FB(Vh.Th,T,P,val);}
 // void D2_BF(const R2 & P,RNMK_ & val) const ; //{ tfe->D2_FB(Vh.Th,T,P,val);}
  void Pi_h(RN_ val,InterpolFunction f,R *v,  void * arg) const;//  {tfe->Pi_h(Vh.Th,T,val,f,v);}
  aIPJ Pi_h_ipj() const { return tfe->Ph_ijalpha();}
  aR2 Pi_h_R2() const { return tfe->Pi_h_R2();}
  void Pi_h(KN_<R> & v) const { return tfe->Pi_h_alpha(*this,v);}
  
  int NbDoF() const { return tfe->NbDoF;}
  int DFOnWhat(int i) const { return tfe->DFOnWhat[i];}
  int FromDF(int i) const { return tfe->fromDF[i];}
  int FromFE(int i) const { return tfe->fromFE[i];}
 // df is the df in element 
  int NodeOfDF(int df) const { return tfe->NodeOfDF[df];} // a node
  int FromASubFE(int i) const { return tfe->fromASubFE[i];}
  int FromASubDF(int i) const { return tfe->fromASubDF[i];}
   int DFOfNode(int df) const { return tfe->DFOfNode[df];} // the df number on the node 

  //  add Mai 2008 for optimization 
  int dfcbegin(int ic) const { return this->tfe->begin_dfcomp[ic];}
  int dfcend(int ic) const { return this->tfe->end_dfcomp[ic];}
 
  R operator()(const R2 & PHat,const KN_<R> & u,int i,int op)  const ;
  complex<R> operator()(const R2 & PHat,const KN_<complex<R> > & u,int i,int op)  const ;
  
 // FElementGlobalToLocal operator()(const KN_<R> & u ) const { return FElementGlobalToLocal(*this,u);}
  private:
  int nbsubdivision() const { return tfe->nbsubdivision;} // for draw 
  };
  
  
class aSubFMortar { public:
  R a,b; 
  R2 A,B; 
  int left,right; //
  int TLeft() const { return left/3;}
  int ELeft() const{ return left%3;}
  
  int TRight() const{ return right/3;}
  int ERight()const{ return right%3;}
  
  int * DfNumberOFmul;
  int  Nbmul; 
 // int *whatnode[2];
  R (**f)(const FESpace *,const aSubFMortar * ,R); // array of function of mul on aSubMortar 
  R lg1(){throwassert(a<b);return b-a;}  
  aSubFMortar():DfNumberOFmul(0),Nbmul(0),f(0){}
 
  R2 operator()(const Mesh & Th,R x,int side) const
   {
      int kk = side==0 ? left : right;
      int k=kk/3;
      int e=kk%3;
      const Triangle & K(Th[k]);
      int c = VerticesOfTriangularEdge[e][0];
      int d = VerticesOfTriangularEdge[e][1];
      
      const R2 &C=K[c];
      const R2 &D=K[d];
      // $ [AB] include in  [CD] $
      R2 CD(C,D);
      R  CD2=(CD,CD);
      R la= (CD,A-C)/CD2;
      R lb= (CD,B-C)/CD2;
      //  A = C*(1-la)+ la*D
      //  B = C*(1-lb)+ lb*D
      //  X = (1-x)*A+x*B;
      //  X = (1-x)*(C*(1-la)+ la*D) + x*(C*(1-lb)+ lb*D)
      //  X = C*((1-x)*(1-la)+x*(1-lb))
      R xx=((1-x)*(1-la)+x*(1-lb));
      throwassert(xx>=0 && xx <= 1);
      return TriangleHat[c]*(xx)+TriangleHat[d]*(1-xx);
   }
  
};  



class FMortar { public:
  friend class FESpace;
  const FESpace &Vh;  
  const Mortar &M;
  const int *p;
  const int nbn; //  nb of node
//  const int nbnl,nbnr; // nb of node left, right
  
  //   we supposse 
  //  the node numbering in a mortar is
  //  1   the lagrange mul.
  //  2   the node of left side
  //  3   the node of right side
  
  const int N;