fespacen.hpp

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

// ********** DO NOT REMOVE THIS BANNER **********
// ORIG-DATE:     Jan 2008
// -*- Mode : c++ -*-
//
// SUMMARY  : Generic Fiinite Element header  1d, 2d, 3d  
// USAGE    : LGPL      
// ORG      : LJLL Universite Pierre et Marie Curi, Paris,  FRANCE 
// 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

 Thank to the ARN ()  FF2A3 grant
 ref:ANR-07-CIS7-002-01 
 */

#ifndef FESpacen_HPP_
#define FESpacen_HPP_
/*
 *  FEspacen.hpp
 *  EF23n
 *
 *  Created by Fr閐閞ic Hecht on 04/12/07.
 *  Copyright 2007 Universite Pierre et marie Curie  All rights reserved.
 *
 */
#include <cassert>
#include <cmath>
#include <cstdlib>
#include <fstream>
#include <iostream>
#include <complex>
#include <map>
using namespace std;
#include "error.hpp"
#include "ufunction.hpp"
#include "Mesh3dn.hpp"
#include "Mesh2dn.hpp"
#include "Mesh1dn.hpp"

#include "RNM.hpp"


#include "QuadratureFormular.hpp"

namespace Fem2D {

template<class Mesh> class GFESpace;
template<class Mesh> class GFElement;
template<class Mesh> class GbaseFElement;
template<class Mesh> class GTypeOfFE;
 
// 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   
}; 

typedef unsigned int What_d;
 
const unsigned int Fop_id= 1<< op_id;

const unsigned int Fop_dx= 1<< op_dx;
const unsigned int Fop_dy= 1<< op_dy;
const unsigned int Fop_dz= 1<< op_dz;

const unsigned int Fop_dxx= 1<< op_dxx;
const unsigned int Fop_dxy= 1<< op_dxy;
const unsigned int Fop_dxz= 1<< op_dxz;

const unsigned int Fop_dyx= 1<< op_dyx;
const unsigned int Fop_dyy= 1<< op_dyy;
const unsigned int Fop_dyz= 1<< op_dyz;

const unsigned int Fop_dzx= 1<< op_dzx;
const unsigned int Fop_dzy= 1<< op_dzy;
const unsigned int Fop_dzz= 1<< op_dzz;

const unsigned int Fop_D0 = Fop_id;
const unsigned int Fop_D1 = Fop_dx | Fop_dy | Fop_dz;
const unsigned int Fop_D2 = Fop_dxx | Fop_dyy | Fop_dzz | Fop_dxy | Fop_dxz | Fop_dyz;
const unsigned int Fop_Dall = Fop_D0| Fop_D1| Fop_D2;

inline What_d  Fwhatd(const operatortype op) { return 1<< op;}


const int last_operatortype=10;
const bool operatortypeValue[last_operatortype]= {true,false,false,false,false,false,false,false,false,false} ; 


inline void initwhatd(bool *whatd,int k)
{
    for (int i=0;i<last_operatortype;i++)
	whatd[i]=false;    
    whatd[k]=true;
}

typedef double R;
typedef KN_<R> RN_;
typedef KN<R>  RN;
typedef KNM_<R> RNM_;
typedef KNMK_<R> RNMK_;
typedef KNMK<R>  RNMK;

template<class Mesh> class GFElement; 
template<class Mesh> class GFESpace; 
template<class Mesh>class GTypeOfFE ;

class dataTypeOfFE {
    

private:
  const int * data; 
  const int * dataalloc; 
    
public:
    const int ndfonVertex; 
    const int ndfonEdge;
    const int ndfonFace;
    const int ndfonVolume;
     const int NbDoF;
    const int NbNode;
    int  N,nb_sub_fem;  
    const int nbsubdivision; // nb of subdivision for plot    
    const bool discontinue;

    int const * const DFOnWhat; 
    int const * const DFOfNode; // nu  du df on Node
    int const * const NodeOfDF; // nu du node du df
    int const * const fromFE;   //  the df  come from df of FE
    int const * const fromDF;   //  the df  come from with FE
    int const * const fromASubFE; //  avril 2006 for CL
    int const * const fromASubDF; //  avril 2006 for CL
    int const * const dim_which_sub_fem; // from atomic sub FE for CL
    const  int * ndfOn() const { return & ndfonVertex;}
    
  dataTypeOfFE(const int *nnitemdim,const int dfon[4],int NN,int nbsubdivisionn,int nb_sub_femm=1,bool discon=true);
    // pour evite un template 
    // nitemdim : nbitem : si d==2   3,3,1,0 , si d=3: 4,6,4,1 , si d==1 = 2,1,0,0 
    // dfon : nombre de df par item 
    // NN 
    dataTypeOfFE(const int nitemdim[4],const KN< dataTypeOfFE const *>  & tef);
    
    virtual ~dataTypeOfFE(){ if(dataalloc) delete [] dataalloc;}
 };

template<class RdHat>
class InterpolationMatrix {
public:
  const int N,np,ncoef; 
  bool invariant;
  int k; 
  KN<RdHat> P; 
  KN<R> coef;
  KN<int> comp;
  KN<int> p;
  KN<int> dofe; 

  template<class Mesh>
  InterpolationMatrix(const GFESpace<Mesh> &Vh);

  template<class Mesh>
  InterpolationMatrix(const GTypeOfFE<Mesh> & tef);

  void set(int k);


private: // copie interdit ...
  InterpolationMatrix(const InterpolationMatrix &);
  void operator=(const InterpolationMatrix &);
};

template <class RdHat>
ostream & operator<<(ostream& f,const InterpolationMatrix<RdHat> &M)
{ f<< M.N << " "<< M.k << " "<< M.np << " "<< M.ncoef << endl;
  f<< " = " << M.P ;
  f << "coef=" <<M.coef ;
  f << "comp " << M.comp;
  f << "p = " << M.p;
  f << "dofe = " << M.dofe;
  return f;
}

/*
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) {}
};
inline ostream & operator<<( ostream & f,const IPJ & ipj)
{ f << ipj.i << ' '<< ipj.p << ' '<< ipj.j << " ,  ";
  return f;}
*/
/*
template<class RdHat>
class GFEInterpole 
{
public:
  KN<IPJ> pij_alpha ;
  KN<RdHat> P_Pi_h ;
  double *coef_Pi_h_alpha;
  bool clean;
  
  GFEInterpole(int kPi,int npPi,double * coef_Pi_h_a,bool cc)
    : pij_alpha(kPi),
      P_Pi_h(npPi),
      coef_Pi_h_alpha(coef_Pi_h_a),
      clean(cc)
  {}
  
  GFEInterpole(const KN< GFEInterpole< RdHat> const  *>  & tef, const KN<int> &tefN)
    : pij_alpha(),P_Pi_h(),coef_Pi_h_alpha(0),clean(true)
  {
    
    map<RdHat,int,lessRd> mpt;

    int np=0,nc=0;
    for(int i=0;i<tef.N();++i)
      {
	np += tef[i]->P_Pi_h.N();
	nc += tef[i]->pij_alpha.N();
      }
    
    KN<int>  kpt(np);
    int npp=0,kkk=0;
    for(int i=0,k=0;i<tef.N();++i)
      for(int j=0;j<tef[i]->P_Pi_h.N();++j,++k)
	{	
	  RdHat P(tef[i]->P_Pi_h[j]);
	  if( mpt.find(P) == mpt.end())
	    mpt[P]=npp++; 
	  kpt[kkk++]=mpt[P]; 
	}
    
    pij_alpha.init(nc);
    P_Pi_h.init(npp);
    for(int i=0,k=0;i<tef.N();++i)
      for(int j=0;j<tef[i]->P_Pi_h.N();++j,++k)
	P_Pi_h[kpt[k]]=tef[i]->P_Pi_h[j];
    int op=0,oj=0,oi=0;
    for(int i=0,k=0;i<tef.N();++i)
      {
	for(int j=0;j<tef[i]->pij_alpha.N();++j,++k)
	  {
	    IPJ a =tef[i]->pij_alpha[j];
	    pij_alpha[k]=IPJ(a.i+oi,kpt[a.p+op],a.j+oj);
	  }
	op += tef[i]->P_Pi_h.N();
	oi += tef[i]->pij_alpha.N();
	oj += tefN[i];
      }
    
    cout << "GFEInterpole pij_alpha : n Pt " << npp << " " << pij_alpha.N() <<  endl;
    for(int i=0;i<nc;++i)
      cout << pij_alpha[i] << endl;
    coef_Pi_h_alpha = new double[nc]; 
    assert(op==np);
    assert(oi==nc);
  }
      
  const KN<IPJ > & Ph_ijalpha() const {return pij_alpha;} // la struct de la matrice
  const KN<RdHat> & Pi_h_Rd() const { return P_Pi_h;}   // les points
  
  
  virtual ~GFEInterpole() 
  {
    if(clean) delete[] coef_Pi_h_alpha;
  }
  
  };
  
*/


template<class Mesh>
class GTypeOfFE : public  dataTypeOfFE
{ 
public:
  typedef typename Mesh::Element Element;
  typedef  Element E;
  typedef typename Element::RdHat RdHat;
  typedef typename Element::Rd Rd;
  typedef GFElement<Mesh> FElement;

  // generic data build interpolation
  bool  invariantinterpolationMatrix;
  int NbPtforInterpolation;
  int NbcoefforInterpolation;
  KN<RdHat> PtInterpolation; 
  KN<int> pInterpolation,cInterpolation,dofInterpolation; 
  KN<R>  coefInterpolation;
  
    // 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
 


  KN<GTypeOfFE<Mesh> *> Sub_ToFE; 
  KN<int> begin_dfcomp, end_dfcomp;
  KN<int> first_comp,last_comp; // for each SubFE 
    
  virtual void init(InterpolationMatrix<RdHat> & M,FElement * pK,int odf,int ocomp,int *pp) const
  {
    assert(pK==0);
    assert(M.np==NbPtforInterpolation);
    assert(M.ncoef==NbcoefforInterpolation);
    M.P=PtInterpolation;
    M.coef=coefInterpolation;
    M.comp=cInterpolation;
    M.p=pInterpolation;
    M.dofe=dofInterpolation;
  }

  
  // the full constructor ..
  GTypeOfFE(const int dfon[4],const int NN,int  nsub,int kPi,int npPi,bool invar,bool discon) 
    : 
    dataTypeOfFE(Element::nitemdim,dfon, NN, nsub,1,discon),
    invariantinterpolationMatrix(invar),
    NbPtforInterpolation(npPi),
    NbcoefforInterpolation(kPi),
    PtInterpolation(NbPtforInterpolation), 
    pInterpolation(NbcoefforInterpolation),
    cInterpolation(NbcoefforInterpolation),
    dofInterpolation(NbcoefforInterpolation),
    coefInterpolation(NbcoefforInterpolation),

    Sub_ToFE(this->nb_sub_fem),
    begin_dfcomp(N,0),
    end_dfcomp(N,this->NbDoF),
    first_comp(this->nb_sub_fem,0),
    last_comp(this->nb_sub_fem,N)

    { 
      Sub_ToFE=this;
      assert(this->dim_which_sub_fem[this->N-1]>=0 && this->dim_which_sub_fem[this->N-1]< this->nb_sub_fem);
    } 

    //  simple constructeur of lagrange type Finite element   1 point par Node for interpolation
  GTypeOfFE(const int dfon[4],const int NN,int  nsub,bool invar,bool discon)
    : 
    dataTypeOfFE(Element::nitemdim,dfon, NN, nsub,1,discon),
    
    invariantinterpolationMatrix(invar),
    NbPtforInterpolation(this->NbDoF),
    NbcoefforInterpolation(this->NbDoF),
    PtInterpolation(NbPtforInterpolation),
    pInterpolation(NbcoefforInterpolation),
    cInterpolation(NbcoefforInterpolation),
    dofInterpolation(NbcoefforInterpolation),
    coefInterpolation(NbcoefforInterpolation),
    
    Sub_ToFE(this->nb_sub_fem),
    begin_dfcomp(N,0),
    end_dfcomp(N,this->NbDoF),
    first_comp(this->nb_sub_fem,0),
    last_comp(this->nb_sub_fem,N)    
  { 
    Sub_ToFE=this;
    assert(this->dim_which_sub_fem[this->N-1]>=0 && this->dim_which_sub_fem[this->N-1]< this->nb_sub_fem);
  } 

  //  simple constructeur pour sum direct d'espace
  GTypeOfFE(const KN<GTypeOfFE const  *>  & tef)
    : 
    dataTypeOfFE(Element::nitemdim,tef),

    invariantinterpolationMatrix(0),
    NbPtforInterpolation(0),
    NbcoefforInterpolation(ncoeftef(tef)),
    PtInterpolation(NbPtforInterpolation),
    pInterpolation(NbcoefforInterpolation),
    cInterpolation(NbcoefforInterpolation),
    dofInterpolation(NbcoefforInterpolation),
    coefInterpolation(NbcoefforInterpolation),


    Sub_ToFE(this->nb_sub_fem),
    begin_dfcomp(this->N,0),
    end_dfcomp(this->N,this->NbDoF),
    first_comp(this->nb_sub_fem,0),
    last_comp(this->nb_sub_fem,N)        
  { 
    Sub_ToFE=this;
    assert(this->dim_which_sub_fem[this->N-1]>=0 && this->dim_which_sub_fem[this->N-1]< this->nb_sub_fem);
  } 
  
  
  virtual R operator()(const GFElement<Mesh>  & K,const  RdHat & PHat,const KN_<R> & u,int componante,int op) const  ;
  virtual void FB(const What_d whatd,const Mesh & Th,const Element & K,const Rd &P, KNMK_<R> & val) const =0;
  
  static KN<int> tefN(const KN<GTypeOfFE const  *>  & tef) {
    int n=tef.N();
    KN<int> ntef(n);
    for(int i=0;i<n;++i) ntef[i]=tef[i]->N;
    return ntef;}

  static int ncoeftef(const KN<GTypeOfFE const  *>  & tef) {
    int k=0,n=tef.N();
    for(int i=0;i<n;++i)
      k+= tef[i]->NbcoefforInterpolation;
    return k;}


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++)