fespace-v0.cpp

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

#include <cmath>
#include <cstdlib>
#include "error.hpp"
#include <iostream>
#include <fstream>
#include <map>
#include "rgraph.hpp"
using namespace std;  

#include "RNM.hpp"
#include "fem.hpp"
#include "FESpace.hpp"

namespace  Fem2D {

 int  Make(const TypeOfFE ** t,int k,KN<R2> & P,KN<int> & I)
 {
   typedef  TypeOfFE::IPJ IPJ;

   int n=0,nn=0;

   for (int i=0;i<k;i++)
    {
     const KN<R2> p(t[i]->P_Pi_h);
     for (int j=0;j<p.N();j++,nn++)
     {
       P[n]=p[j]; // par defaut un nouveau => ajout
       I[nn]=n++;       
       for (int l=0;l<n-1;l++)
        {
          R2 QP(p[j],P[l]);
          if ( (QP,QP) < 1.0e-12 ) { 
            I[nn]=l;
            n--; // on a retrouver =>  detruit
            break;}
        }
       
     }
    }
   return n; // nombre de point commun
 }
 
 KN<TypeOfFE::IPJ > Makepij_alpha(const TypeOfFE ** t,int k)
 {
 // Attention les df est numerote de facon croissant 
 // en faisant une boucle sur les TypeOfFE
 // comme dans la class TypeOfFESum
  typedef TypeOfFE::IPJ IPJ;
  int n=0,m=0;
  for (int i=0;i< k;i++) {
    n += t[i]->pij_alpha.N();
    m += t[i]->P_Pi_h.N();
    }
  KN<TypeOfFE::IPJ > ij(n);
  KN<int> I(m);
  KN<R2> P(m);
  Make(t,k,P,I);
  int p0=0,i0=0,N0=0,nn=0;
  for (int i=0;i< k;i++) {
    const KN<IPJ > p(t[i]->pij_alpha);
     for (int j=0;j<p.N();j++,nn++)  
      {
        ij[nn].i= p[j].i+i0; //  comme dans TypeOfFESum
        ij[nn].j= N0+ p[j].j;
        ij[nn].p= I[p[j].p+p0];
       
      //   cout << nn << "Makepij_alpha: " << ij[nn].i << " " << ij[nn].p << " " << ij[nn].j << endl;
        
      }
    i0+=t[i]->NbDoF;
    p0+=t[i]->P_Pi_h.N();
    N0+=t[i]->N;}
  return ij; 
 }
 KN<R2 > MakeP_Pi_h(const TypeOfFE **t,int k)
 {
  int np=0;
  for (int i=0;i< k;i++) 
    np += t[i]->P_Pi_h.N();
    
  KN< R2 >  yy(np);
  KN<int> zz(np);
  int kk=Make(t,k,yy,zz);
 // cout << " MakeP_Pi_h: " << kk << " from " << np << endl; 
  return yy(SubArray(kk));
 
 }

ListOfTFE * ListOfTFE::all ; // list of all object of this type 

ListOfTFE::ListOfTFE (const char * n,TypeOfFE *t) : name(n),tfe(t) 
{
  if(!t)
  assert(t);
  static int count=0;
  if (count++==0) 
    all=0; // init of all in dependant of the ordre of the objet file   
  next=all;
  all=this;
}

const TypeOfFE ** Make(const FESpace **l,int k) {
  const TypeOfFE** p=new const TypeOfFE*[k];
  for (int i=0;i<k;i++)
    p[i]=l[i]->TFE[0];
  return p;
}
const TypeOfFE ** Make(const TypeOfFE **l,int k) {
  const TypeOfFE** p=new const TypeOfFE*[k];
  for (int i=0;i<k;i++)
    p[i]=l[i];
  return p;
}

  
bool Same(const FESpace **l,int k)
{
   for (int i=1;i<k;i++)
    if (l[0] != l[i] ) return false;
   return true;
}   

class FESumConstruct { protected:
   const TypeOfFE ** teb;
   const int k;
   int nbn; // nb of node 
   int  *  data;
   int  * const NN; //  NN[ i:i+1[ dimension de l'element i
   int  * const DF; // DF[i:i+1[  df associe a l'element i
   int  * const comp; //  
   FESumConstruct(int kk,const TypeOfFE **t);   
   virtual ~FESumConstruct(){
     delete [] DF;
     delete [] NN;
     delete [] comp;
     delete [] data;}   
};

  void FElement::Pi_h(RN_ val,InterpolFunction f,R *v, void * arg=0) const {
   // routine: a  tester FH.
    FElement::aIPJ ipj(Pi_h_ipj()); 
    FElement::aR2  PtHat(Pi_h_R2()); 
    KN<R>   Aipj(ipj.N());
    KNM<R>  Vp(N,PtHat.N());
        
     Pi_h(Aipj);
     for (int p=0;p<PtHat.N();p++)
          { 
            f(v,T(PtHat[p]),*this,T.lab,PtHat[p],arg);
            KN_<double> Vpp(Vp('.',p));
            for (int j=0;j<N;j++)          
               Vpp[j]=v[j];              
           }
           
         for (int i=0;i<Aipj.N();i++)
          { 
           const FElement::IPJ &ipj_i(ipj[i]);
           val[ipj_i.i] += Aipj[i]*Vp(ipj_i.j,ipj_i.p);           
          }
 }

class TypeOfFESum: public FESumConstruct, public  TypeOfFE { public:
   TypeOfFESum(const FESpace **t,int kk): 
     FESumConstruct(kk,Make(t,kk)),TypeOfFE(teb,kk,data) {}
       TypeOfFESum(const TypeOfFE **t,int kk): 
     FESumConstruct(kk,Make(t,kk)),TypeOfFE(teb,kk,data) {}

  // void FB(const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
   void FB(const bool * whatd,const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
//   void D2_FB(const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
//  void Pi_h(const baseFElement & K,RN_ & val, InterpolFunction f, R* v,int, void * arg ) const; 
   virtual void Pi_h_alpha(const baseFElement & K,KN_<double> & v) const
    { 
      for (int i=0,k0=0;i<k;i++) {
          int n=teb[i]->NbDoF;
          KN_<R> sv(v(SubArray(n,k0)));
          teb[i]->Pi_h_alpha(K,sv);
          k0+= n;}
    }
   ~TypeOfFESum(){  delete []  teb;}
} ;

class FEProduitConstruct { protected:
   const TypeOfFE & teb;
   int k;
   int * data;
   FEProduitConstruct(int kk,const TypeOfFE &t)  ;   
   ~FEProduitConstruct(){delete [] data;}   
};

class TypeOfFEProduit: protected FEProduitConstruct, public  TypeOfFE { public:  
   TypeOfFEProduit(int kk,const TypeOfFE &t): 
     FEProduitConstruct(kk,t),TypeOfFE(t,kk,data)  {}
     
  // void FB(const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
   void FB(const bool * whatd,const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
//   void D2_FB(const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
//   void Pi_h(const baseFElement & K,RN_ & val, InterpolFunction f, R* v,int, void * arg ) const;  
   virtual void Pi_h_alpha(const baseFElement & K,KN_<double> & v) const
    { int nbof=teb.NbDoF;
      for (int i=0,k0=0;i<k;i++,k0+=nbof)
        {          
          KN_<R> sv(v(SubArray(nbof,k0)));
          teb.Pi_h_alpha(K,sv);
        }
    }
 
   ~TypeOfFEProduit(){}
} ;

FEProduitConstruct::FEProduitConstruct(int kk,const TypeOfFE &t)
 :k(kk),teb(t) 
{  
  int m= teb.NbDoF;
  KN<int> nn(teb.NbNode);
  nn=0; // nb de dl par noeud 
  for (int i=0;i<m;i++)
    nn[teb.NodeOfDF[i]]++;   
    
  int n= m*kk;
  int N= teb.N*kk;
  data = new int [n*5+N];
  int c=0;
  
  for (int i=0;i<m;i++)
   for (int j=0;j<kk;j++)
    data[c++] = teb.DFOnWhat[i];
  
  for (int i=0;i<m;i++)
   for (int j=0;j<kk;j++) // num of df on node  for the df = j 
      data[c++] = teb.DFOfNode[i]+j*nn[teb.NodeOfDF[i]];
    
  for (int i=0;i<m;i++)
   for (int j=0;j<kk;j++) 
    data[c++] = teb.NodeOfDF[i]; //  node of df

  for (int i=0;i<m;i++)
   for (int j=0;j<kk;j++) 
    data[c++] = j; //  node from of FE
    
  for (int i=0;i<m;i++)
   for (int j=0;j<kk;j++) 
    data[c++] = i; //  node from of df in FE
    
   for (int j=0;j<kk;j++) 
    for (int i=0;i<teb.N;i++) 
      data[c++]= teb.dim_which_sub_fem[i] + teb.nb_sub_fem*j ;
}

FESumConstruct::FESumConstruct(int kk,const TypeOfFE **t)
 :k(kk),teb(t),NN(new int[kk+1]),DF(new int[kk+1]) , comp(new int[kk])
{  
   map<const TypeOfFE *,int> m;
   int i=k,j;    
   while(i--) // on va a l'envert pour avoir comp[i] <=i 
      m[teb[i]]=i;
    // l'ordre comp est important comp est croissant  mais pas de pb. 
   i=k;    
   while(i--) 
     comp[i]=m[teb[i]]; //  comp[i] <=i
        
  // reservatition des intervalles en espaces
  int n=0,N=0;
   for ( j=0;j<kk;j++)
     {NN[j]=N;N+=teb[j]->N;}
   NN[kk] = N;
 //  reservation des interval en df   
   n=0;
   for ( j=0;j<kk;j++)
    { DF[j]=n;n+=teb[j]->NbDoF;}
   DF[kk] = n;
//  n = nb de DF total   
//  N the fem is in R^N 
   
  data = new int [n*5 + N];
  int c=0;
  int ki= 0; 
// recherche des noeuds
   KN<int> w(7),nn(7);
   w=0;
   nn=0; 
   
         
   for ( j=0;j<kk;j++)
     for ( i=0;i<teb[j]->NbDoF;i++)
         nn[teb[j]->DFOnWhat[i]]++;
   nbn=0;      
   for( j=0;j<7;j++)
     if (nn[j]) nn[j]=nbn++;
     else nn[j]=-1;
   KN<int> dln(7);
   dln=0;
  // nn donne numero de noeud sur what            
   for ( j=0;j<kk;j++)
     for ( i=0;i<teb[j]->NbDoF;i++)
       data[c++] = teb[j]->DFOnWhat[i];
    
   for ( j=0;j<kk;j++)
    {
     int  cc=c;
     for ( i=0;i<teb[j]->NbDoF;i++)
       data[c++] = teb[j]->DFOfNode[i]+dln[teb[j]->DFOnWhat[i]];
     for ( i=0;i<teb[j]->NbDoF;i++)
       dln[teb[j]->DFOnWhat[i]]=Max(dln[teb[j]->DFOnWhat[i]],data[cc++]+1);      
    }
        
       
   for ( j=0;j<kk;j++)
    { 
     //  w renumerotation des noeuds 
     //  Ok si un noeud par what 
     for ( i=0;i<teb[j]->NbDoF;i++)
       data[c++] = nn[teb[j]->DFOnWhat[i]];
    }
    
   for ( j=0;j<kk;j++)
     for ( i=0;i<teb[j]->NbDoF;i++)
       data[c++] = j; //  node from of FE

    
   for ( j=0;j<kk;j++)
     for ( i=0;i<teb[j]->NbDoF;i++)
       data[c++] = i; //  node from of df in FE
       
   int xx=0;
   for (j=0;j<kk;j++)
     { 
      int xxx=xx;
      for (i=0;i<teb[j]->N;i++)
       { 
         data[c] = teb[j]->dim_which_sub_fem[i]+xx;
         xxx=Max(xxx,data[c]+1);
         c++;
       }
       xx=xxx;
     }
  
  throwassert(c== 5*n+N);      
/*  int cc=0;
   cout << " Data : " << endl;
  for ( i=0;i<5;i++)    {
    for (j=0;j<n;j++)
      cout << " " << data[cc++];
     cout << endl;}
 cout << " which " ;
 for (i=0;i<N;i++)
   cout << " " << data[cc++];
  cout << endl;*/
}

class TypeOfFE_P1Lagrange : public  TypeOfFE { public:  
  static int Data[];
  static double Pi_h_coef[];
   TypeOfFE_P1Lagrange(): TypeOfFE(1,0,0,1,Data,1,1,3,3,Pi_h_coef)
    { const R2 Pt[] = { R2(0,0), R2(1,0), R2(0,1) }; 
      for (int i=0;i<NbDoF;i++) {
       pij_alpha[i]= IPJ(i,i,0);
       P_Pi_h[i]=Pt[i]; }
     }
  // void FB(const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
   void FB(const bool * whatd,const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
   
//   void D2_FB(const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
  // void Pi_h(const baseFElement & K,RN_ & val, InterpolFunction f, R* v,int, void *) const;
virtual R operator()(const FElement & K,const  R2 & PHat,const KN_<R> & u,int componante,int op) const ;
   
} ;

///////////////////////////////////////////////////////////////////////////////
////////////////////////////////// NEW ////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

class TypeOfFE_P1Bubble : public  TypeOfFE { public:  
  static int Data[];
  static double Pi_h_coef[];
   TypeOfFE_P1Bubble(): TypeOfFE(1,0,1,1,Data,1,1,4,4,Pi_h_coef)
    { const R2 Pt[] = { R2(0,0), R2(1,0), R2(0,1), R2(1./3.,1./3.) }; 
      for (int i=0;i<NbDoF;i++) {
       pij_alpha[i]= IPJ(i,i,0);
       P_Pi_h[i]=Pt[i]; }
     }
  // void FB(const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
   void FB(const bool * whatd,const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
   
//   void D2_FB(const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
 //  void Pi_h(const baseFElement & K,RN_ & val, InterpolFunction f, R* v,int, void *) const;
  //virtual R operator()(const FElement & K,const  R2 & PHat,const KN_<R> & u,int componante,int op) const ;
   
} ;
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

class TypeOfFE_P2Lagrange : public  TypeOfFE { public:  
  static int Data[];
  static double Pi_h_coef[];
  
   TypeOfFE_P2Lagrange(): TypeOfFE(1,1,0,1,Data,3,1,6,6,Pi_h_coef)
    { const R2 Pt[] = { R2(0,0), R2(1,0), R2(0,1),R2(0.5,0.5),R2(0,0.5),R2(0.5,0) };
      for (int i=0;i<NbDoF;i++) {
       pij_alpha[i]= IPJ(i,i,0);
       P_Pi_h[i]=Pt[i]; }
     }
   
  // void FB(const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
   void FB(const bool * whatd,const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
 //  void D2_FB(const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
  // void Pi_h(const baseFElement & K,RN_ & val, InterpolFunction f, R* v,int, void *) const;
} ;

int TypeOfFE_P1Lagrange::Data[]={0,1,2,       0,0,0,       0,1,2,       0,0,0,        0,1,2,       0};
int TypeOfFE_P1Bubble::Data[]={0,1,2,6,     0,0,0,0,     0,1,2,3,     0,0,0,0,        0,1,2,3,     0};
int TypeOfFE_P2Lagrange::Data[]={0,1,2,3,4,5, 0,0,0,0,0,0, 0,1,2,3,4,5, 0,0,0,0,0,0,  0,1,2,3,4,5, 0};
double TypeOfFE_P1Lagrange::Pi_h_coef[]={1.,1.,1.};
double TypeOfFE_P1Bubble::Pi_h_coef[]={1.,1.,1.,1.};
double TypeOfFE_P2Lagrange::Pi_h_coef[]={1.,1.,1.,1.,1.,1.};

inline void dump(char *m,int n,int * p)
{
  cout << m ;
  for (int i=0;i<n;i++) cout << " " << p[i] ;
  cout << endl;
}



ConstructDataFElement::~ConstructDataFElement()
{
  if(*counter==0) 
   {
    delete [] NodesOfElement;
    delete []  FirstNodeOfElement;