element_p4.cpp

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

#include "error.hpp"
#include "AFunction.hpp"
#include "rgraph.hpp"
using namespace std;  
#include "RNM.hpp"
#include "fem.hpp"
#include "FESpace.hpp"
#include "AddNewFE.h"
// Attention probleme de numerotation des inconnues
// -------------------------------------------------
// dans freefem, il y a un noeud par objets  sommet, arete, element.
// et donc la numerotation des dl dans l'element depend 
// de l'orientation des aretes
// 
/// ---------------------------------------------------------------
namespace  Fem2D {

 // ------ P4  Hierarchical (just remove P1 node of the P2 finite element)  --------
 class TypeOfFE_P4Lagrange : public  TypeOfFE { public:  
   static const int k=4; 
   static const int ndf = (k+2)*(k+1)/2;
   static int Data[];
   static double Pi_h_coef[];
   static const int nn[15][4] ;
   static const int aa[15][4] ;
   static const int ff[15];
   static const int il[15];
   static const int jl[15];
   static const int kl[15];
    
   
   TypeOfFE_P4Lagrange(): TypeOfFE(3+3*3+3,1,Data,4,1,15+6,15,0)
   {  
     static const  R2 Pt[15] =  {
       R2( 0/4. , 0/4. ) , 
       R2( 4/4. , 0/4. ) , 
       R2( 0/4. , 4/4. ) , 
       R2( 3/4. , 1/4. ) , 
       R2( 2/4. , 2/4. ) , 
       R2( 1/4. , 3/4. ) , 
       R2( 0/4. , 3/4. ) , 
       R2( 0/4. , 2/4. ) , 
       R2( 0/4. , 1/4. ) , 
       R2( 1/4. , 0/4. ) , 
       R2( 2/4. , 0/4. ) , 
       R2( 3/4. , 0/4. ) , 
       R2( 1/4. , 2/4. ) , 
       R2( 2/4. , 1/4. ) , 
       R2( 1/4. , 1/4. ) } 
     ;
     
                       //    3,4,5, 6,7,8, 9,10,11, 
     int other[15]= { 0,1,2, 5,4,3, 8,7,6, 11,10,9,12,13,14}; 
     
     int kk=0;
     for (int i=0;i<NbDoF;i++) 
       {
	 pij_alpha[kk++]= IPJ(i,i,0);
	 if(other[i]!=i)
	   pij_alpha[kk++]= IPJ(i,other[i],0);
	 P_Pi_h[i]=Pt[i]; 
       }

      assert(P_Pi_h.N()==NbDoF);
      assert(pij_alpha.N()==kk);
   }
  
   void FB(const bool * whatd, const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
   void Pi_h_alpha(const baseFElement & K,KN_<double> & v) const
   {
     for (int i=0;i<15+6;++i)
       v[i]=1;
     int e0=K.EdgeOrientation(0);
     int e1=K.EdgeOrientation(1);
     int e2=K.EdgeOrientation(2);
     int ooo[6]={e0,e0,e1,e1,e2,e2};
     /*   3,4
	  5,
	  6,7
	  8,9,
	  10,
	  11,12,
	  13,14,
	  15
	  16,17
     */     
     int iii[6]={3,6,8,11,13,16}; 
     int jjj[6];
     for(int i=0;i<6;++i)
       { 
         jjj[i]= iii[i]+1; // si orient = -1
       }
     for(int i=0;i<6;++i)
       if(ooo[i]==1) v[jjj[i]]=0;
       else v[iii[i]]=0;
          
   }
 } ;
  //                     on what     nu df on node node of df    
  int TypeOfFE_P4Lagrange::Data[]={
    0,1,2,3,3,3,4,4,4,5,5,5,6,6,6,     // the support number  of the node of the df 
    0,0,0,0,1,2,0,1,2,0,1,2,0,1,2,     // the number of the df on  the node  
    0,1,2,3,3,3,4,4,4,5,5,5,6,6,6,    // the node of the df 
    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,     //  the df come from which FE (generaly 0)
    0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,    //  which are de df on sub FE
    0,                      // for each compontant $j=0,N-1$ it give the sub FE associated 
    0,15 };
  
  void TypeOfFE_P4Lagrange::FB(const bool * whatd,const Mesh & ,const Triangle & K,const R2 & P,RNMK_ & val) const
  {
 
    R2 A(K[0]), B(K[1]),C(K[2]);
    R l0=1.-P.x-P.y,l1=P.x,l2=P.y; 
    R L[3]={l0*k,l1*k,l2*k};
    throwassert( val.N()>=10);
    throwassert(val.M()==1);
    // Attention il faut renumeroter les fonction de bases
    //   car dans freefem++, il y a un node par sommet, arete or element
    //   et la numerotation naturelle  mais 2 noud pas arete
    //   donc p est la perumation
    //   echange de numerotation si les arete sont dans le mauvais sens 
    int p[15];
    for(int i=0;i<15;++i)
      p[i]=i;
    
    if(K.EdgeOrientation(0) <0) Exchange(p[3],p[5]);// 3,4
    if(K.EdgeOrientation(1) <0) Exchange(p[6],p[8]);// 5,6 
    if(K.EdgeOrientation(2) <0) Exchange(p[9],p[11]);// 7,8
    //cout << KN_<int>(p,10) <<endl;
    val=0; 
    /*
    //  les fonction de base du Pk Lagrange sont 
    //  
    //
    */
// -- 
    
    
    if (whatd[op_id])
      {
	RN_ f0(val('.',0,op_id)); 
	for (int df=0;df<ndf;df++)
	  {
	    int pdf=p[df];
	    R f=1./ff[df];
	    
	    for( int i=0;i<k;++i)
	      {
		f *= L[nn[df][i]]-aa[df][i];
		//cout <<  L[nn[df][i]]-aa[df][i]<< " ";
	      }
	    f0[pdf] = f;
	    //cout << pdf<< " " << df << " f " <<f <<endl;
	  }
	//cout <<" L " << L[0] << " " << L[1] << " " << L[2] << endl;
	//cout << ndf << " nbf = "<< f0 <<endl;
      }
    
    
    if(  whatd[op_dx] || whatd[op_dy] || whatd[op_dxx] || whatd[op_dyy] ||  whatd[op_dxy])
      {
	R2 D[]={K.H(0)*k, K.H(1)*k,K.H(2)*k };
	if (whatd[op_dx] || whatd[op_dy] )
	  {
	    for (int df=0;df<ndf;df++)
	      {
		int pdf=p[df];
		R fx=0.,fy=0.,f=1./ff[df];
		for( int i=0;i<k;++i) 
		  {
		    int n= nn[df][i];
		    R Ln=L[n]-aa[df][i];
		    fx= fx*Ln+f*D[n].x;
		    fy= fy*Ln+f*D[n].y;
		    f = f*Ln;
		  } 
		if(whatd[op_dx]) val(pdf,0,op_dx)=fx;
		if(whatd[op_dy]) val(pdf,0,op_dy)=fy;		
	      } 
	  }
	
	if (whatd[op_dyy] ||whatd[op_dxy]|| whatd[op_dxx] )
	  {  
	    for (int df=0;df<ndf;df++)
	      {
		int pdf=p[df];		
		R fx=0.,fy=0.,f=1./ff[df];
		R fxx=0.,fyy=0.,fxy=0.;
		for( int i=0;i<k;++i) 
		  {
		    int n= nn[df][i];
		    R Ln=L[n]-aa[df][i];
		    fxx=fxx*Ln+2.*fx*D[n].x;
		    fyy=fyy*Ln+2.*fy*D[n].y;
		    fxy=fxy*Ln+fx*D[n].y+fy*D[n].x;
		    fx= fx*Ln+f*D[n].x;
		    fy= fy*Ln+f*D[n].y;
		    f = f*Ln;
		  } 
		if(whatd[op_dxx]) val(pdf,0,op_dxx)=fxx;
		if(whatd[op_dyy]) val(pdf,0,op_dyy)=fyy;
		if(whatd[op_dxy]) val(pdf,0,op_dxy)=fxy;
	      }
	  }
	
      }
  }
#include "Element_P4.hpp"
  


// link with FreeFem++ 
static TypeOfFE_P4Lagrange P4LagrangeP4;
// a static variable to add the finite element to freefem++
static AddNewFE  P4Lagrange("P4",&P4LagrangeP4); 
} // FEM2d namespace 


// --- fin --