fespace.hpp

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

   TypeOfMortar const *  const tom; 
  FMortar(const FESpace * VVh,int k) ;
  public:
  int  operator[](int i) const ;
  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));}
  int NbDoF(int i) const ;//{return tom->NbDoF(Vh.Th,M,i);};  // number of DF 
  int NbOfNodes()const {return nbn;}
  int NbDoF() const ;
  int NodeOfDF(int i) const ;
  int DFOfNode(int i) const ;
  
  int nbsm; // nb of submortar
  aSubFMortar * sm;
  ~FMortar() { 
    delete [] datai;
    delete [] dataf;}
  
  int *datai;
   R (**dataf)(const FESpace *,const aSubFMortar *,R);
  private:
   FMortar(const FMortar &);
   void operator=(const FMortar &);

};

extern TypeOfFE & P1Lagrange;
extern TypeOfFE & P2Lagrange;
extern TypeOfFE & RTLagrange;
extern TypeOfFE & RTLagrangeOrtho;
extern TypeOfFE & P0Lagrange;
extern TypeOfFE & P1ncLagrange;

  class FESpace : public RefCounter , public UniqueffId
{
  public:
  typedef Fem2D::Mesh Mesh;
  typedef Fem2D::FElement FElement;
  typedef Mesh::Element  Element;
  typedef Mesh::Rd  Rd;
  typedef Fem2D::TypeOfFE TypeOfFE;
  typedef Fem2D::QuadratureFormular QFElement;
  typedef Fem2D::QuadratureFormular1d QFBorderElement;


  const Mesh &Th;
  TypeOfFE const * const ptrTFE;
  KN<const TypeOfFE *>  TFE; 
  private:
  ConstructDataFElement * cdef; //  juste pour les constantes 
  public:
  CountPointer<const Mesh> cmesh;
  const int N; // dim espace d'arrive
  const int Nproduit; // dim de l'espace produit generalement 1
  const int NbOfDF;
  const int NbOfElements;
  const int NbOfNodes;
  const int nb_sub_fem; // nb de sous elements finis tensorise (independe au niveau des composantes)
  int const* const dim_which_sub_fem;// donne les dependant des composantes liee a un meme sous element fini
   //   exemple si N=5,  
   // dim_which_sub_fem[0]=0;
   // dim_which_sub_fem[1] =1;
   // dim_which_sub_fem[2]= 2
   // dim_which_sub_fem[3] =2 
   // dim_which_sub_fem[4] =3
   //  =>
   //  le  sous  elements fini 0 est lie a la composante 0 
   //  le  sous  elements fini 1 est lie a la composante 1 
   //  le  sous  elements fini 2 est lie aux composantes 2,3 
   //  le  sous  elements fini 3 est lie a la composante 4 
   //  donc pour les CL. les composante 2 et 3 sont lie car elle sont utiliser pour definir un
   //  meme degre de libert\u00e9.
  
  int const*const NodesOfElement;
  int const*const FirstNodeOfElement;
  int const*const FirstDfOfNodeData;
  const TypeOfMortar * tom; 
  const int MaxNbNodePerElement;
  const int MaxNbDFPerElement;

//  par defaut P1                
  FESpace(const Mesh & TTh) 
    :
    Th(TTh),
    ptrTFE(0),
    TFE(1,0,&P1Lagrange),
    cdef(0),
    cmesh(TTh),
    N(1),
    Nproduit(1),
    NbOfDF(TTh.nv),
    NbOfElements(TTh.nt),
    NbOfNodes(TTh.nv),
    nb_sub_fem(TFE[0]->nb_sub_fem),
    dim_which_sub_fem(TFE[0]->dim_which_sub_fem),
    NodesOfElement(0),
    FirstNodeOfElement(0),
    FirstDfOfNodeData(0),
    tom(0),
    MaxNbNodePerElement(3),
    MaxNbDFPerElement(3*Nproduit)
 {}


  int FirstDFOfNode(int i) const {return FirstDfOfNodeData ? FirstDfOfNodeData[i] : i*Nproduit;}
  int LastDFOfNode(int i)  const {return FirstDfOfNodeData ? FirstDfOfNodeData[i+1] : (i+1)*Nproduit;}
  int NbDFOfNode(int i)  const {return FirstDfOfNodeData ? FirstDfOfNodeData[i+1]-FirstDfOfNodeData[i] : Nproduit;}
  int MaximalNbOfNodes() const {return MaxNbNodePerElement;};
  int MaximalNbOfDF() const {return MaxNbDFPerElement;};
  const int * PtrFirstNodeOfElement(int k) const {
      return NodesOfElement 
               ? NodesOfElement + (FirstNodeOfElement ? FirstNodeOfElement[k] : k*MaxNbNodePerElement)                 
               : 0;}   
               
  int SizeToStoreAllNodeofElement() const {  
       return  FirstNodeOfElement 
                ?  FirstNodeOfElement[NbOfElements] 
                : MaxNbNodePerElement*NbOfElements;}   
                      
  int NbOfNodesInElement(int k)   const {             
      return FirstNodeOfElement 
               ?  FirstNodeOfElement[k+1] - FirstNodeOfElement[k] 
               : MaxNbNodePerElement ;}
  int  esize() const { return  MaxNbDFPerElement*N*last_operatortype;}   // size to store all value of B. function        
      
    FESpace(const FESpace &,int k );
    FESpace(const FESpace **,int k ); 
    FESpace(const Mesh & TTh,const TypeOfFE **,int k,int nbdfv=0,const int *ndfv=0,int nbdfe=0,const int *ndfe=0 );//int NbDfOnSommet,int NbDfOnEdge,int NbDfOnElement,int NN=1); 
 
    FESpace(const Mesh & TTh,const TypeOfFE & ,
    int nbdfv=0,const int *ndfv=0,int nbdfe=0,const int *ndfe=0);//int NbDfOnSommet,int NbDfOnEdge,int NbDfOnElement,int NN=1); 
    
    FESpace(const Mesh & TTh,const TypeOfFE &,const TypeOfMortar & );//int NbDfOnSommet,int NbDfOnEdge,int NbDfOnElement,int NN=1); 
  ~FESpace();    
 // FESpace(Mesh & TTh,int NbDfOnSommet,int NbDfOnEdge,int NbDfOnElement,int NN=1); 
  int  renum();
      
  FElement operator[](int k) const { return FElement(this,k);} 
  FElement operator[](const Triangle & K) const { return FElement(this,Th.number(K));} 
  int  operator()(int k)const {return NbOfNodesInElement(k);}
  int  operator()(int k,int i) const { //  the node i of element k
     return NodesOfElement ?  *(PtrFirstNodeOfElement(k) + i)  : Th(k,i)  ;}
  
  void Draw(const KN_<R>& U,const KN_<R>& Viso,int j=0,float *colors=0,int nbcolors=0,bool hsv=true,bool drawborder=true) const ; // Draw iso line
  void Drawfill(const KN_<R>& U,const KN_<R>& Viso,int j=0,double rapz=1,float *colors=0,int nbcolors=0,bool hsv=true,bool drawborder=true) const ; // Draw iso line
 
  KN<R>  newSaveDraw(const KN_<R> & U,int composante,int & lg,int & nsb) const   ; 
  KN<R>  newSaveDraw(const KN_<R> & U,const KN_<R> & V,int iU,int IV,int & lg,int & nsb) const   ; 
  void Draw(const KN_<R>& U,const KN_<R> & Viso, R coef,int j0=0,int j1=1,float *colors=0,int nbcolors=0,bool hsv=true,bool drawborder=true) const  ; // Arrow
  void Draw(const KN_<R>& U,const KN_<R>& V,const KN_<R> & Viso, R coef,int iu=0,int iv=0,float *colors=0,int nbcolors=0,bool hsv=true,bool drawborder=true) const  ; // Arrow
  R2 MinMax(const KN_<R>& U,const KN_<R>& V,int j0,int j1,bool bb=true) const ;
  R2 MinMax(const KN_<R>& U,int j0, bool bb=true) const ;
 // void destroy() {RefCounter::destroy();}
  bool isFEMesh() const { return !cdef && ( N==1) ;} // to make optim
  void Show() const {  
 // cout << " Show: FESpace " << this << " " <<  N << " ";  if(cdef) cout << cdef->NodesOfElement << endl;else cout << endl;
  }

 private: // for gibbs  
  int gibbsv (long* ptvoi,long* vois,long* lvois,long* w,long* v);
    
};

inline baseFElement::baseFElement(  const FESpace &aVh, int k) 
    : Vh(aVh),T(Vh.Th[k]),tfe(aVh.TFE[k]),N(aVh.N),number(k){}
    
inline baseFElement::baseFElement(const   baseFElement & K,  const TypeOfFE & atfe) 
    : Vh(K.Vh),T(K.T),tfe(&atfe),N(Vh.N),number(K.number){}


inline FElement::FElement(const FESpace * VVh,int k) 
  : baseFElement(*VVh,k) ,
    p(Vh.PtrFirstNodeOfElement(k)),
    nb(Vh.NbOfNodesInElement(k))
    
     {} 
  
inline   int  FElement::operator[](int i) const {
   return  p ? p[i] :  ((&T[i])-Vh.Th.vertices);}  
   
inline   int  FElement::operator()(int i,int df) const {
   return  Vh.FirstDFOfNode(p ? p[i] :  ((&T[i])-Vh.Th.vertices)) + df;}  
   
inline   int  FMortar::operator()(int i,int df) const {throwassert(p);
   return  Vh.FirstDFOfNode(p[i]) + df;} 
    
inline   int  FMortar::operator[](int i) const {throwassert(p);
   return  p[i];}  
   
inline   int  FElement::NbDoF(int i) const {
   int node =p ? p[i] :  ((&T[i])-Vh.Th.vertices);
   return  Vh.LastDFOfNode(node)-Vh.FirstDFOfNode(node);}  

void  SetDefaultIsoValue(const KN_<R>& U,KN_<R> & Viso);
void  SetDefaultIsoValue(const KN_<R>& u,const KN_<R>& v,KN_<R> & Viso);
void MoveTo(R2 P); 
void LineTo(R2 P) ;

/*
void operator=(  KN_<R> & u,const FElementGlobalToLocal & x) 
{
  int n=u.N();
  throwassert(n==x.S.NbDoF());
  for (int i=0;i<n;i++) // get the local value
     v[i] = x.U[x.S(i)];
}
*/
inline  int FMortar::NbDoF(int i) const {
   int node = p[i];
   return  Vh.LastDFOfNode(node)-Vh.FirstDFOfNode(node);
};  // number of DF 
inline  int FMortar::NbDoF() const { return tom->NbDoF(Vh.Th,M);}
//inline  int FMortar::NbOfNodes()const {return }
inline  int FMortar::NodeOfDF(int i) const { return tom->NodeOfDF(Vh,M,i);}
inline  int FMortar::DFOfNode(int i) const { return tom->DFOfNode(Vh,M,i);}

inline ostream & operator << (ostream & f,const FElement & FE)
   {
     f << FE.number << "," <<FE.nb << ":" ;
     for (int i=0;i<FE.nb;i++) f << "\t"<<FE.p[i];
   return f;
     
   }
inline ostream & operator << (ostream & f,const FESpace & Vh)
   {
     cout << " list of nodes per element :" << endl;
     for (int k=0;k< Vh.NbOfElements;k++)
      { f <<setw(3) <<  k << ":";
      for (int j=0;j<Vh(k);j++)
        f << " " << setw(3) << Vh(k,j);
        cout << endl;
      }
      
      f << endl << " FirstDFOfNode :" ;       
      for (int i=0;i<=Vh.NbOfNodes;i++)
        {if (i%10==0)  cout << "\n" << setw(3) << i << " : ";
        cout << setw(3) << Vh.FirstDFOfNode(i) << " ";}
    
   return f;
     
   }
   
inline void FElement::BF(const R2 & P,RNMK_ & val) const {
 static bool whatdold[last_operatortype]={true,true,true,false,false,false,false,false,false,false};
 tfe->FB(whatdold,Vh.Th,T,P,val);}
inline void FElement::BF(const bool * whatd,const R2 & P,RNMK_ & val) const { tfe->FB(whatd,Vh.Th,T,P,val);}
//inline  void FElement::D2_BF(const R2 & P,RNMK_ & val) const { tfe->D2_FB(Vh.Th,T,P,val);}


//  -------
 extern const TypeOfMortar & TheMortarCas1P2; 
 
void PlotValue(const RN_ & Viso,int  k0,const char * cmm);
 
// to store all the type of TFE
// the problem is the TFE can be define on lot of file.cpp
struct ListOfTFE { 
  const char * name;
  TypeOfFE * tfe;
  ListOfTFE * next;

  static ListOfTFE * all ; // list of all object of this type 
  ListOfTFE (const char * n,TypeOfFE *t);
};
// to get a unique list of TypeOfFE 
// local variable of TypeOfFE
ListOfTFE & GetListOfTFE() ;


inline   R FElement::operator()(const R2 & PHat,
                                const KN_<R> & u,int i,int op)  const
{
 return (*tfe)(*this,PHat,u,i,op);
}


inline  complex<R> FElement::operator()(const R2 & PHat,const KN_<complex<R> > & u,int i,int op)  const 
{
 complex<double> * pu=u; // pointeur du tableau
  double *pr = static_cast<double*>(static_cast<void*>(pu));

  const KN_<R>  ur(pr,u.n,u.step*2);
  const KN_<R>  ui(pr+1,u.n,u.step*2);
  
   return complex<R>((*tfe)(*this,PHat,ur,i,op),(*tfe)(*this,PHat,ui,i,op));
}

}


#endif