fespace.hpp

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
 */

#ifndef FESpace_h_
#define FESpace_h_
#include "FESpacen.hpp"

#include <cstdarg>

namespace  Fem2D {
  class FESpace; 
class TypeOfFE;
  /*
// numbering of derivative 
enum operatortype { op_id=0,
   op_dx=1,op_dy=2,
   op_dxx=3,op_dyy=4,
   op_dyx=5,op_dxy=5,   
   op_dz=6,
   op_dzz=7,     
   op_dzx=8,op_dxz=8, 
   op_dzy=9,op_dyz=9   
   }; 
    
const int last_operatortype=10;
 

inline void initwhatd(bool *whatd,int k)
{
  for (int i=0;i<last_operatortype;i++)
    whatd[i]=false;    
  whatd[k]=true;
}
  */
typedef KN_<R> RN_;
typedef KN<R>  RN;
typedef KNM_<R> RNM_;
//typedef KNM<R>  RNM;
typedef KNMK_<R> RNMK_;
typedef KNMK<R>  RNMK;

inline int FindIn(int v,int *a,int n)
{
  for (int i=0;i<n;i++)
    if( v==a[i]) return i;
  return -1;
};


//  Methode boeuf 
//    on suppose que pour chaque elment fini
// on sait calculer 
//   un tableau de point x,y 
// les valeur de la fonction et des derives
// FB(RN_<R2> 
// i  =  df
// j  = [0 N[  (ou N est la dim de l'espace d'arrive N=3 
// k = 0,1,2   f,fx,fy 


class FElement;
class baseFElement;
class FMortar;
class TypeOfMortar;

//  for FE
//typedef void  (* basicFunction)(const FElement & FE,const R2 &P, RNMK_ & val);
typedef void  (* InterpolFunction)(R* v, const R2 & P,const baseFElement &,int where,const R2 & Phat, void * arg);
//typedef void (*InterpolantFunction)(const FElement & FE,RN_ & val, InterpolFunction f, R* v);
//  for FM   ( Finite Mortar  Mortar + interpolation)
typedef void  (* basicFunctionMortar)(const FMortar & FE,const R2 &P, RNMK_ & val);
typedef void  (* InterpolFunctionMortar)(R* v, const R2 & P,const Mortar &,int where,const R2 & Phat);
typedef void  (*InterpolantFunctionMortar)(const FMortar & FE,RN_ & val, InterpolFunctionMortar f, R* v);

//void P1Functions(const FElement &FE,const R2 & P,RNMK_ & val);
//void P2Functions(const FElement &FE,const R2 & P,RNMK_ & val);


class VofScalarFE { 
  R v[3];
  public:
  VofScalarFE() {v[0]=v[1]=v[2]=0;}
  VofScalarFE(R f,R fx,R fy) {v[0]=f;v[1]=fx;v[2]=fy;}

  R & operator[](int i){ return v[i];}
  const R & operator[](int i) const { return v[i];}  
  R f()  { return v[0];}
  R fx() { return v[1];}
  R fy() { return v[2];}
 VofScalarFE &operator+=(const VofScalarFE & a) { v[0]+=a.v[0]; v[1]+=a.v[1]; v[2]+=a.v[2];return *this;}  
};

class ConstructDataFElement {
  friend class FESpace;
  //int thecounter; //  chang 09/2008 the counter is a pointer because 
    //  if you remove before the master the counter become invalide. 
  int * counter;
  int MaxNbNodePerElement;
  int MaxNbDFPerElement;
  int *NodesOfElement;
  int *FirstNodeOfElement;
  int *FirstDfOfNode;
  int NbOfElements;
  int NbOfDF;
  int NbOfNode;
  int Nproduit;
  ConstructDataFElement(const Mesh &Th,/*int NbDfOnSommet,int NbDfOnEdge,int NbDfOnElement,*/
  const  KN<const TypeOfFE *> & TFEs,const TypeOfMortar *tm=0,
  int nbdfv=0,const int *ndfv=0,int nbdfe=0,const int *ndfe=0);
  ConstructDataFElement(const ConstructDataFElement *,int k);
  ConstructDataFElement(const FESpace ** l,int k,const  KN<const TypeOfFE *> & TFEs) ;
  void renum(const long *r,int l)  ; 
  ~ConstructDataFElement();
 void Make(const Mesh &Th,/*int NbDfOnSommet,int NbDfOnEdge,int NbDfOnElement*/const  KN<const TypeOfFE *> & TFEs,const TypeOfMortar *tm=0,
   int nbdfv=0,const int *ndfv=0,int nbdfe=0,const int *ndfe=0);
private:
    static int *NewCounter() { int * p=new int; *p=0;return p;}// add the build thecounter. 
};

template<class T>
inline int sum(const T ** l,int const T::*p,int n)
  {
    int r=0;
    for (int i=0;i<n;i++)
     r += l[i]->*p;
     return r;
  }
  
template<class T>
inline int max(const T ** l,int const T::*p,int n)
  {
    int r=0;
    for (int i=0;i<n;i++)
     r =Max( l[i]->*p,r);
     return r;
  }
 
 
class TypeOfFE { public:
//  The FEM is in R^N
//  The FEM is compose from nb_sub_fem
//  dim_which_sub_fem[N] give
  friend class FESpace;
  friend class FElement;
  friend class FEProduitConstruct;
  const int NbDoF;
  const int NbNodeOnVertex, NbNodeOnEdge, NbNodeOnElement;
  const int NbDfOnVertex, NbDfOnEdge, NbDfOnElement, N,nb_sub_fem;
  const  int NbNode;
  //  remark 
 // virtual void FB(const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const =0;
  virtual void FB(const bool *,const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const =0;
//  virtual void D2_FB(const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const=0;
 // virtual void Pi_h(const baseFElement & K,RN_ & val, InterpolFunction f, R* v,int , void *) const=0;
  
  // soit 
  //  $(U_pj)_{{j=0,N-1}; {p=0,nbpoint_Pi_h-1}}$ les valeurs de U au point Phat[i];
  //   p est le numero du point d'integration
  //   j est la composante 
  //  l'interpole est  defini par
  //  Pi_h u = \sum_k u_k w^k , et u_i = \sum_pj alpha_ipj U_pj
  //  la matrice de alpha_ipj est tres creuse
  struct IPJ { 
    int i,p,j; // i is DoF, p is Point, j is componante 
    IPJ(int ii=0,int pp=0,int jj=0):i(ii),p(pp),j(jj) {}
  };
  
  friend KN<IPJ> Makepij_alpha(const TypeOfFE **,int );
  friend KN<R2> MakeP_Pi_h(const TypeOfFE **,int );
 
   const KN<IPJ > & Ph_ijalpha() const {return pij_alpha;} // la struct de la matrice
   const KN<R2> & Pi_h_R2() const { return P_Pi_h;}   // les points
   virtual void Pi_h_alpha(const baseFElement & K,KN_<double> & v) const
     {  assert(coef_Pi_h_alpha);
        v=KN_<double>(coef_Pi_h_alpha,pij_alpha.N());}

   //  ----  
  
  const int nbsubdivision; // nb of subdivision for plot
   int const * const DFOnWhat; //  0,1,2 vertex 3,4,5 edge 6  triangle
   int const * const DFOfNode; // n\u00c9 du df on Node
   int const * const NodeOfDF; // 
   int const * const fromFE;   //  the df  come from df of FE
   int const * const fromDF;   //  the df  come from with FE
   int const * const dim_which_sub_fem; // from atomic sub FE for CL (size N)
   KN<IPJ > pij_alpha ;
   KN<R2 > P_Pi_h ;
  double *coef_Pi_h_alpha;
   KN<TypeOfFE *> Sub_ToFE; //  List of atomic sub TFE avril 2006 
    //  form Atomic Sub FE 
   int const * const fromASubFE; //  avril 2006 for CL
   int const * const fromASubDF; //  avril 2006 for CL
  int const * const  begin_dfcomp; //  mai 2008 for optimiation 
  int const  * const end_dfcomp; // mai 2008 
   
  // if 0  no plot 
  public: 
  
    TypeOfFE(const TypeOfFE & t,int k,const int * data,const int * data1) 
    : 
      NbDoF(t.NbDoF*k),
      NbNodeOnVertex(t.NbNodeOnVertex),NbNodeOnEdge(t.NbNodeOnEdge),NbNodeOnElement(t.NbNodeOnElement),
      NbDfOnVertex(t.NbDfOnVertex*k),NbDfOnEdge(t.NbDfOnEdge*k),NbDfOnElement(t.NbDfOnElement*k),
      N(t.N*k),nb_sub_fem(t.nb_sub_fem*k),
      NbNode(t.NbNode), 
     nbsubdivision(t.nbsubdivision),
    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(t.pij_alpha.N()*k),P_Pi_h(t.P_Pi_h),
    coef_Pi_h_alpha(0),
    Sub_ToFE(nb_sub_fem),
    fromASubFE(data1+0*NbDoF),
      fromASubDF(data1+1*NbDoF),
      begin_dfcomp(data1+2*NbDoF),
      end_dfcomp(data1+2*NbDoF+N)

    { 
       for(int i=0,kk=0;i<k;i++)
        for (int j=0;j<t.nb_sub_fem;j++)
          Sub_ToFE[kk++]=t.Sub_ToFE[j];
       assert(begin_dfcomp[0]==0 && end_dfcomp[N-1]==NbDoF);
       
      throwassert(dim_which_sub_fem[N-1]>=0 && dim_which_sub_fem[N-1]< nb_sub_fem);
     // Warning the componant is moving first 
        for (int j=0,l=0;j<t.pij_alpha.N();j++) // for all sub DF
          for(int i=0; i<k; i++,l++) // for componate
          {
          pij_alpha[l].i=t.pij_alpha[j].i*k+i;   //  DoF number
          pij_alpha[l].p=t.pij_alpha[j].p;       //  point of interpolation
          pij_alpha[l].j=t.pij_alpha[j].j+i*t.N; //  componante of interpolation
          }                         
    } 
       
    TypeOfFE(const TypeOfFE ** t,int k,const int * data,const int * data1) 
    : 
      NbDoF(sum(t,&TypeOfFE::NbDoF,k)),
      NbNodeOnVertex(NbNodebyWhat(data,NbDoF,0)),
      NbNodeOnEdge(NbNodebyWhat(data,NbDoF,3)),
      NbNodeOnElement(NbNodebyWhat(data,NbDoF,6)),
      NbDfOnVertex(sum(t,&TypeOfFE::NbDfOnVertex,k)),
      NbDfOnEdge(sum(t,&TypeOfFE::NbDfOnEdge,k)),
      NbDfOnElement(sum(t,&TypeOfFE::NbDfOnElement,k)),
      N(sum(t,&TypeOfFE::N,k)),nb_sub_fem(sum(t,&TypeOfFE::nb_sub_fem,k)),
      NbNode( (NbDfOnVertex ? 3 :0) + (NbDfOnEdge ? 3 :0 ) +(NbDfOnElement? 1 :0) ),      
      nbsubdivision(max(t,&TypeOfFE::nbsubdivision,k)),
      
    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(Makepij_alpha(t,k)),
    P_Pi_h(MakeP_Pi_h(t,k)),
    coef_Pi_h_alpha(0),
    Sub_ToFE(nb_sub_fem),
    fromASubFE(data1+0*NbDoF),
      fromASubDF(data1+1*NbDoF),
      begin_dfcomp(data1+2*NbDoF),
      end_dfcomp(data1+2*NbDoF+N)    
  {
    for(int i=0,kk=0;i<k;i++)
      for (int j=0;j<t[i]->nb_sub_fem;j++)