element_rt.cpp

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

	    R2 B(TriangleHat[VerticesOfTriangularEdge[i][1]]);
	     Pt[k++]=A*(QF_GaussLegendre2[j].x)+ B*(1-QF_GaussLegendre2[j].x)
	  }
	      
        int other[6]= { 1,0, 3,2,5,4 };
	k=0;
	  for (int i=0;i<NbDoF;i++) {
	    pij_alpha[i]= IPJ(i,i,0);
	    if(other[i]>=0)
		pij_alpha[kk++]= IPJ(i,other[i],0);
	       P_Pi_h[i]=Pt[i]; }
	}
	
	void FB(const bool * whatd, const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
	
    } ;
    //                     on what     nu df on node node of df    
    int TypeOfFE_ConsEdge::Data[]={3,4,5,       0,0,0,       0,1,2,       0,0,0,        0,1,2,       0, 0,3};
    double TypeOfFE_ConsEdge::Pi_h_coef[]={1.,1.,1.};
    void TypeOfFE_ConsEdge::FB(const bool * whatd,const Mesh & ,const Triangle & K,const R2 & P,RNMK_ & val) const
    {
	//  const Triangle & K(FE.T);
	R2 A(K[0]), B(K[1]),C(K[2]);
	R l0=1-P.x-P.y,l1=P.x,l2=P.y; 
	
	if (val.N() <3) 
	    throwassert(val.N() >=3);
	throwassert(val.M()==1 );
	
	val=0; 
	if (whatd[op_id])
	{
	    
	    RN_ f0(val('.',0,0)); 
	    //      
	    f0[0] =  double(l0 <= min(l1,l2) ); // arete  
	    f0[1] =  double(l1 <= min(l0,l2) );
	    f0[2] =  double(l2 <= min(l0,l1) );
	}
	
    }
    
*/
    
void TypeOfFE_P1ncLagrange::FB(const bool * whatd,const Mesh & ,const Triangle & K,const R2 & P,RNMK_ & val) const
{
  //  const Triangle & K(FE.T);
  R2 A(K[0]), B(K[1]),C(K[2]);
  //  l1(  cshrink1*(cshrink*((1,0)-G)+G)-G)+G  = 1
  R l0=1-P.x-P.y,l1=P.x,l2=P.y; 
  
  if (val.N() <3) 
    throwassert(val.N() >=3);
  throwassert(val.M()==1 );
  // throwassert(val.K()==3 );
  
  val=0; 
  if (whatd[op_id])
   {
  RN_ f0(val('.',0,0)); 
  f0[0] = 1-l0*2;  
  f0[1] = 1-l1*2;
  f0[2] = 1-l2*2;
  }
  if (whatd[op_dx] || whatd[op_dy] )
   {
    R2 Dl0(K.H(0)), Dl1(K.H(1)), Dl2(K.H(2));
    if (whatd[op_dx]) 
      {
      RN_ f0x(val('.',0,op_dx)); 
      f0x[0] = -Dl0.x*2;
      f0x[1] = -Dl1.x*2;
      f0x[2] = -Dl2.x*2;
      }
    if (whatd[op_dy]) 
     {  
     RN_ f0y(val('.',0,op_dy)); 
     f0y[0] = -Dl0.y*2;
     f0y[1] = -Dl1.y*2;
     f0y[2] = -Dl2.y*2;
     }
   }
}

// The RT orthogonal FE

class TypeOfFE_RTortho : public  TypeOfFE { public:  
  static int Data[];
   TypeOfFE_RTortho(): TypeOfFE(0,1,0,2,Data,1,1,6,3) 
     {const R2 Pt[] = { R2(0.5,0.5), R2(0.0,0.5), R2(0.5,0.0) };
      for (int p=0,kk=0;p<3;p++)
       { P_Pi_h[p]=Pt[p];   
        for (int j=0;j<2;j++) 
        pij_alpha[kk++]= IPJ(p,p,j);
       }}

   void FB(const bool * watdd, const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
   void Pi_h_alpha(const baseFElement & K,KN_<double> & v) const ;
} ; 
//                     on what     nu df on node node of df    
  int TypeOfFE_RTortho::Data[]={3,4,5,       0,0,0,       0,1,2,       0,0,0,        0,1,2,   0,0, 0,0, 3,3};



 void TypeOfFE_RTortho::FB(const bool *whatd,const Mesh & Th,const Triangle & K,const R2 & PHat,RNMK_ & val) const
{ //  
//  const Triangle & K(FE.T);
  R2 P(K(PHat));
  R2 A(K[0]), B(K[1]),C(K[2]);
  //R l0=1-P.x-P.y,l1=P.x,l2=P.y; 
 // R2 Dl0(K.H(0)), Dl1(K.H(1)), Dl2(K.H(2));
  if (val.N() <3) 
   throwassert(val.N() >=3);
  throwassert(val.M()==2 );
//  throwassert(val.K()==3 );
  val=0;  
  R a=1./(2*K.area);
  R a0=   K.EdgeOrientation(0) * a ;
  R a1=   K.EdgeOrientation(1) * a  ;
  R a2=   K.EdgeOrientation(2) * a ;
 // if (Th(K)< 2) cout << Th(K) << " " <<  A << " "  << B << " " << C << "; " <<  a0 << " " << a1 << " "<< a2 << endl;;

  //  ------------
  if (whatd[op_id])
   {
   assert(val.K()>op_id);
  RN_ f0(val('.',0,0)); 
  RN_ f1(val('.',1,0)); 
  f1[0] =  (P.x-A.x)*a0;
  f0[0] = -(P.y-A.y)*a0;
  
  f1[1] =  (P.x-B.x)*a1;
  f0[1] = -(P.y-B.y)*a1;
  
  f1[2] =  (P.x-C.x)*a2;
  f0[2] = -(P.y-C.y)*a2;
  }
  // ----------------
    if (whatd[op_dx])
   {
   assert(val.K()>op_dx);
   val(0,1,op_dx) =  a0;  
   val(1,1,op_dx) =  a1;  
   val(2,1,op_dx) =  a2; 
  } 
    if (whatd[op_dy])
   {
   assert(val.K()>op_dy);
    val(0,0,op_dy) =  -a0;  
    val(1,0,op_dy) =  -a1;  
    val(2,0,op_dy) =  -a2;  
  }
}

void TypeOfFE_RTortho::Pi_h_alpha(const baseFElement & K,KN_<double> & v) const 
{
  const Triangle & T(K.T);

   for (int i=0,k=0;i<3;i++)
     {  
        R2 E(T.Edge(i));
        R signe = &T[ (i+1)%3] < &T[ (i+2)%3] ? 1.0 : -1.0 ;
        v[k++]= signe*E.x;
        v[k++]= signe*E.y;
     }   
}
// -------------------
// ttdc finite element fully discontinue. 
// -------------------
class TypeOfFE_P1ttdc : public  TypeOfFE { public:  
  static int Data[];
  static double Pi_h_coef[];
    static const R2 G;
    static const R cshrink;
    static const R cshrink1;
    //  (1 -1/3)*
    
    static R2 Shrink(const R2& P){ return (P-G)*cshrink+G;}
    static R2 Shrink1(const R2& P){ return (P-G)*cshrink1+G;}
    
   TypeOfFE_P1ttdc(): TypeOfFE(0,0,3,1,Data,1,1,3,3,Pi_h_coef)
    { const R2 Pt[] = { Shrink(R2(0,0)), Shrink(R2(1,0)), Shrink(R2(0,1)) }; 
      for (int i=0;i<NbDoF;i++) {
       pij_alpha[i]= IPJ(i,i,0);
       P_Pi_h[i]=Pt[i];
       // cout << Pt[i] << " " ;
      }
      //	cout <<" cshrink: " << cshrink << " cshrink1 : "<< cshrink1 <<endl;
     }

   void FB(const bool * whatd,const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
   
  
  virtual R operator()(const FElement & K,const  R2 & PHat,const KN_<R> & u,int componante,int op) const ;
   
} ;
    const R2 TypeOfFE_P1ttdc::G(1./3.,1./3.);   
    const R TypeOfFE_P1ttdc::cshrink=1-1e-2;   
    const R TypeOfFE_P1ttdc::cshrink1=1./TypeOfFE_P1ttdc::cshrink;   
    
class TypeOfFE_P2ttdc : public  TypeOfFE { public:  
  static int Data[];
  static double Pi_h_coef[];
    static const R2 G;
    static const R cshrink;
    static const R cshrink1;
    
    static R2 Shrink(const R2& P){ return (P-G)*cshrink+G;}
    static R2 Shrink1(const R2& P){ return (P-G)*cshrink1+G;}
    
   TypeOfFE_P2ttdc(): TypeOfFE(0,0,6,1,Data,3,1,6,6,Pi_h_coef)
    { const R2 Pt[] = { Shrink(R2(0,0)), Shrink(R2(1,0)), Shrink(R2(0,1)),
	                Shrink(R2(0.5,0.5)),Shrink(R2(0,0.5)),Shrink(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 bool * whatd,const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
 
 
} ;
//                          on what   nu df on node  node of df    
  int TypeOfFE_P1ttdc::Data[]={6,6,6,       0,1,2,       0,0,0,       0,0,0,         0,1,2,       0, 0,3};
  int TypeOfFE_P2ttdc::Data[]={6,6,6,6,6,6, 0,1,2,3,4,5, 0,0,0,0,0,0,  0,0,0,0,0,0,  0,1,2,3,4,5, 0, 0,6};
double TypeOfFE_P1ttdc::Pi_h_coef[]={1.,1.,1.};
double TypeOfFE_P2ttdc::Pi_h_coef[]={1.,1.,1.,1.,1.,1.};
    
const R2 TypeOfFE_P2ttdc::G(1./3.,1./3.);   
const R TypeOfFE_P2ttdc::cshrink=1-1e-2;   
const R TypeOfFE_P2ttdc::cshrink1=1./TypeOfFE_P2ttdc::cshrink;   
    
 
R TypeOfFE_P1ttdc::operator()(const FElement & K,const  R2 & P1Hat,const KN_<R> & u,int componante,int op) const 
{ 

  R2 PHat=Shrink1(P1Hat);  
  R u0(u(K(0))), u1(u(K(1))), u2(u(K(2)));
  R r=0;
  if (op==0)
    {
      R l0=1-PHat.x-PHat.y,l1=PHat.x,l2=PHat.y; 
      r = u0*l0+u1*l1+l2*u2;
    }
   else
    { 
       const Triangle & T=K.T;
       R2 D0 = T.H(0)*cshrink1 , D1 = T.H(1)*cshrink1  , D2 = T.H(2)*cshrink1 ;
       if (op==1)
         r =  D0.x*u0 + D1.x*u1 + D2.x*u2 ;
        else 
         r =  D0.y*u0 + D1.y*u1 + D2.y*u2 ;
    }
 //  cout << r << "\t";
   return r;
}

void TypeOfFE_P1ttdc::FB(const bool *whatd,const Mesh & ,const Triangle & K,const R2 & P1,RNMK_ & val) const
{
  R2 P=Shrink1(P1);  
   
  //  const Triangle & K(FE.T);
  R2 A(K[0]), B(K[1]),C(K[2]);
  R l0=1-P.x-P.y,l1=P.x,l2=P.y; 
  
  if (val.N() <3) 
    throwassert(val.N() >=3);
  throwassert(val.M()==1 );
  //  throwassert(val.K()==3 );
  
  val=0; 
  RN_ f0(val('.',0,op_id)); 
  
  if (whatd[op_id]) 
   {
     f0[0] = l0;
    f0[1] = l1;
    f0[2] = l2;}
 if (whatd[op_dx] || whatd[op_dy])
  {
  R2 Dl0(K.H(0)*cshrink1), Dl1(K.H(1)*cshrink1), Dl2(K.H(2)*cshrink1);
  
  if (whatd[op_dx]) 
   {
    RN_ f0x(val('.',0,op_dx)); 
   f0x[0] = Dl0.x;
   f0x[1] = Dl1.x;
   f0x[2] = Dl2.x;
  }
  
  if (whatd[op_dy]) {
    RN_ f0y(val('.',0,op_dy)); 
   f0y[0] = Dl0.y;
   f0y[1] = Dl1.y;
   f0y[2] = Dl2.y;
  }
  }
}



 void TypeOfFE_P2ttdc::FB(const bool *whatd,const Mesh & ,const Triangle & K,const R2 & P1,RNMK_ & val) const
{
    R2 P=Shrink1(P1);  
   
//  const Triangle & K(FE.T);
  R2 A(K[0]), B(K[1]),C(K[2]);
  R l0=1-P.x-P.y,l1=P.x,l2=P.y; 
  R l4_0=(4*l0-1),l4_1=(4*l1-1),l4_2=(4*l2-1); 
  
//  throwassert(FE.N == 1);  
  throwassert( val.N()>=6);
  throwassert(val.M()==1);
//  throwassert(val.K()==3 );
  
  val=0; 
// --     
 if (whatd[op_id])
  {
   RN_ f0(val('.',0,op_id)); 
  f0[0] = l0*(2*l0-1);
  f0[1] = l1*(2*l1-1);
  f0[2] = l2*(2*l2-1);
  f0[3] = 4*l1*l2; // oppose au sommet 0
  f0[4] = 4*l0*l2; // oppose au sommet 1
  f0[5] = 4*l1*l0; // oppose au sommet 3
  }
 if(  whatd[op_dx] || whatd[op_dy] || whatd[op_dxx] || whatd[op_dyy] ||  whatd[op_dxy])
 {
   R2 Dl0(K.H(0)*cshrink1), Dl1(K.H(1)*cshrink1), Dl2(K.H(2)*cshrink1);
  if (whatd[op_dx])
  {
    RN_ f0x(val('.',0,op_dx)); 
  f0x[0] = Dl0.x*l4_0;
  f0x[1] = Dl1.x*l4_1;
  f0x[2] = Dl2.x*l4_2;
  f0x[3] = 4*(Dl1.x*l2 + Dl2.x*l1) ;
  f0x[4] = 4*(Dl2.x*l0 + Dl0.x*l2) ;
  f0x[5] = 4*(Dl0.x*l1 + Dl1.x*l0) ;
  }

 if (whatd[op_dy])
  {  
    RN_ f0y(val('.',0,op_dy)); 
  f0y[0] = Dl0.y*l4_0;
  f0y[1] = Dl1.y*l4_1;
  f0y[2] = Dl2.y*l4_2;
  f0y[3] = 4*(Dl1.y*l2 + Dl2.y*l1) ;
  f0y[4] = 4*(Dl2.y*l0 + Dl0.y*l2) ;
  f0y[5] = 4*(Dl0.y*l1 + Dl1.y*l0) ;
  }
 
 if (whatd[op_dxx])
  {  
    RN_ fxx(val('.',0,op_dxx)); 

    fxx[0] = 4*Dl0.x*Dl0.x;
    fxx[1] = 4*Dl1.x*Dl1.x;
    fxx[2] = 4*Dl2.x*Dl2.x;
    fxx[3] =  8*Dl1.x*Dl2.x;
    fxx[4] =  8*Dl0.x*Dl2.x;
    fxx[5] =  8*Dl0.x*Dl1.x;
  }

 if (whatd[op_dyy])
  {  
    RN_ fyy(val('.',0,op_dyy)); 
    fyy[0] = 4*Dl0.y*Dl0.y;
    fyy[1] = 4*Dl1.y*Dl1.y;
    fyy[2] = 4*Dl2.y*Dl2.y;
    fyy[3] =  8*Dl1.y*Dl2.y;
    fyy[4] =  8*Dl0.y*Dl2.y;
    fyy[5] =  8*Dl0.y*Dl1.y;
  }
 if (whatd[op_dxy])
  {  
    assert(val.K()>op_dxy);
    RN_ fxy(val('.',0,op_dxy)); 
  
    fxy[0] = 4*Dl0.x*Dl0.y;
    fxy[1] = 4*Dl1.x*Dl1.y;
    fxy[2] = 4*Dl2.x*Dl2.y;
    fxy[3] =  4*(Dl1.x*Dl2.y + Dl1.y*Dl2.x);
    fxy[4] =  4*(Dl0.x*Dl2.y + Dl0.y*Dl2.x);
    fxy[5] =  4*(Dl0.x*Dl1.y + Dl0.y*Dl1.x);
  }
 
 }
 
}

 
//
// end ttdc
// ------------------

static TypeOfFE_RTortho The_TypeOfFE_RTortho;
static TypeOfFE_RT The_TypeOfFE_RT;
static TypeOfFE_P0 The_TypeOfFE_P0;
static TypeOfFE_P1ttdc The_TypeOfFE_P1ttdc;
static TypeOfFE_P2ttdc The_TypeOfFE_P2ttdc;
static TypeOfFE_RTmodif The_TypeOfFE_RTmodif;
static TypeOfFE_P1ncLagrange The_TypeOfFE_P1nc;
static TypeOfFE_ConsEdge The_TypeOfFE_ConsEdge; // add FH 
TypeOfFE  & RTLagrangeOrtho(The_TypeOfFE_RTortho);
TypeOfFE  & RTLagrange(The_TypeOfFE_RT);
TypeOfFE  & RTmodifLagrange(The_TypeOfFE_RTmodif);
TypeOfFE  & P0Lagrange(The_TypeOfFE_P0);
TypeOfFE  & P1ncLagrange(The_TypeOfFE_P1nc);
TypeOfFE  & P1ttdc(The_TypeOfFE_P1ttdc);
TypeOfFE  & P2ttdc(The_TypeOfFE_P2ttdc);
TypeOfFE  & P0edge(The_TypeOfFE_ConsEdge);

}