fespace-v0.cpp

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

    }
 }     
 
 
/*
 void TypeOfFEProduit::D2_FB(const Mesh & Th,const Triangle & K,const R2 & P,RNMK_ & val) const
 {
   int n=teb.NbDoF;
   int m=teb.N;   
   val=0.0;
   SubArray t(3);
   RNMK_ v(val(SubArray(n,0,k),SubArray(m),t));
   teb.D2_FB(Th,K,P,v);
   for (int i=1;i<k;i++)
     val(SubArray(n,i,k),SubArray(m,m*i),t)=v; 
 } 
*/ 
/*
 void TypeOfFEProduit::FB(const Mesh & Th,const Triangle & K,const R2 & P,RNMK_ & val) const
 {
   int n=teb.NbDoF;
   int m=teb.N;   
   val=0.0;
   SubArray t(3);
   RNMK_ v(val(SubArray(n,0,k),SubArray(m),t));
   teb.FB(Th,K,P,v);
   for (int i=1;i<k;i++)
     val(SubArray(n,i,k),SubArray(m,m*i),t)=v; 
 }
 */
 void TypeOfFEProduit::FB(const bool * whatd,const Mesh & Th,const Triangle & K,const R2 & P,RNMK_ & val) const
 {
   int n=teb.NbDoF;
   int m=teb.N;   
   val=0.0;
   SubArray t(val.K());
   RNMK_ v(val(SubArray(n,0,k),SubArray(m),t));
   teb.FB(whatd,Th,K,P,v);
   for (int i=1;i<k;i++)
     val(SubArray(n,i,k),SubArray(m,m*i),t)=v; 
 }
 
/* 
 void TypeOfFEProduit::Pi_h(const baseFElement & K,RN_ & val, InterpolFunction f, R* v,int j, void * arg) const
 {
   const baseFElement  KK(K,teb);
   int m=teb.N;   
    for(int i=0;i<k;i++)
     { RN_  vv(val(SubArray(m,m*i)));
     teb.Pi_h(KK,vv,f,v,j+i*m,arg);}
 }
 
/* 
 void TypeOfFESum::D2_FB(const Mesh & Th,const Triangle & K,const R2 & P,RNMK_ & val) const
 {
   val=0.0;
   SubArray t(3);
   for (int i=0;i<k;i++)
    {
     int j=comp[i];
     int ni=NN[i];
     int di=DF[i];  
     int i1=i+1; 
     int nii=NN[i1];
     int dii=DF[i1];
     throwassert(ni<nii && di < dii);
     RNMK_ v(val(SubArray(dii-di,di),SubArray(nii-ni,ni),t));     
     if (j<=i)
       teb[i]->D2_FB(Th,K,P,v);       
     else
       v=val(SubArray(DF[j+1]-DF[j],DF[j]),SubArray(NN[j+1]-NN[j],NN[j]),t);     
    } }
*/    
/*
 void TypeOfFESum::FB(const Mesh & Th,const Triangle & K,const R2 & P,RNMK_ & val) const
 {
   val=0.0;
   SubArray t(3);
   for (int i=0;i<k;i++)
    {
     int j=comp[i];
     int ni=NN[i];
     int di=DF[i];  
     int i1=i+1; 
     int nii=NN[i1];
     int dii=DF[i1];
     throwassert(ni<nii && di < dii);
     RNMK_ v(val(SubArray(dii-di,di),SubArray(nii-ni,ni),t));     
     if (j<=i)
       teb[i]->FB(Th,K,P,v);       
     else
       v=val(SubArray(DF[j+1]-DF[j],DF[j]),SubArray(NN[j+1]-NN[j],NN[j]),t);     
    }
 }
*/ 
  void TypeOfFESum::FB(const bool * whatd,const Mesh & Th,const Triangle & K,const R2 & P,RNMK_ & val) const
 {
   val=0.0;
   SubArray t(val.K());
   for (int i=0;i<k;i++)
    {
     int j=comp[i];
     int ni=NN[i];
     int di=DF[i];  
     int i1=i+1; 
     int nii=NN[i1];
     int dii=DF[i1];
     throwassert(ni<nii && di < dii);
     RNMK_ v(val(SubArray(dii-di,di),SubArray(nii-ni,ni),t));     
     if (j<=i)
       teb[i]->FB(whatd,Th,K,P,v);       
     else
       v=val(SubArray(DF[j+1]-DF[j],DF[j]),SubArray(NN[j+1]-NN[j],NN[j]),t);     
    }
 }
/* 
 void TypeOfFESum::Pi_h(const baseFElement & K,RN_ & val, InterpolFunction f, R* v,int jjj, void * arg) const
 {
    for(int i=0;i<k;i++)
     { 
       const baseFElement  KK(K,*teb[i]);
       int dfii=DF[i+1],dfi=DF[i];
       RN_  vv(val(SubArray(dfii-dfi,dfi)));
       teb[i]->Pi_h(KK,vv,f,v,jjj+NN[i],arg);
     }
    // cout << val(SubArray(NbDoF)) << endl;
     
 }
 /*
 void TypeOfFE_P1Lagrange::D2_FB(const Mesh & ,const Triangle & ,const R2 & ,RNMK_ & val) const
{ //  
  val=0;
}
*/
/*
 void TypeOfFE_P2Lagrange::D2_FB(const Mesh & ,const Triangle & K,const R2 & P,RNMK_ & val) const
{ // 2 times derivatives  for error indicator
//  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; 
  R2 Dl0(K.H(0)), Dl1(K.H(1)), Dl2(K.H(2));
  R l4_0=(4*l0-1),l4_1=(4*l1-1),l4_2=(4*l2-1); 
  
  throwassert(val.N() >=6);
  throwassert(val.M()==1 );
  throwassert(val.K()==3 );
  
  val=0; 
  RN_ fxx(val('.',0,0)); 
  RN_ fxy(val('.',0,1)); 
  RN_ fyy(val('.',0,2)); 
  
  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;

  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;

  fxy[0] = 4*Dl0.y*Dl0.y;
  fxy[1] = 4*Dl1.y*Dl1.y;
  fxy[2] = 4*Dl2.y*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);

}
*/
 R TypeOfFE_P1Lagrange::operator()(const FElement & K,const  R2 & PHat,const KN_<R> & u,int componante,int op) const 
{ 
   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) , D1 = T.H(1)  , D2 = T.H(2) ;
       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_P1Lagrange::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 );
//  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)), Dl1(K.H(1)), Dl2(K.H(2));
  
  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;
  }
  }
}

///////////////////////////////////////////////////////////////////////////////
///////////////////////////////// NEW /////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
/*
void TypeOfFE_P1Bubble::FB(const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const
{
  assert(0);
}


void TypeOfFE_P1Bubble::Pi_h(const baseFElement & K,RN_ & val, InterpolFunction f, R* v,int, void *) const
{
  assert(0);
}

/*
R TypeOfFE_P1Bubble::operator()(const FElement & K,const  R2 & PHat,const KN_<R> & u,int componante,int op) const
{
  assert(0);
}
*/

void TypeOfFE_P1Bubble::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, lb=l0*l1*l2*9.; 
  
  if (val.N() <4) 
   throwassert(val.N() >=4);
  throwassert(val.M()==1 );
//  throwassert(val.K()==3 );
  
  val=0; 
  RN_ f0(val('.',0,op_id)); 
  
  if (whatd[op_id]) 
   {
    f0[0] = l0-lb;
    f0[1] = l1-lb;
    f0[2] = l2-lb;
    f0[3] = 3.*lb;
   }
  if(  whatd[op_dx] || whatd[op_dy] || whatd[op_dxx] || whatd[op_dyy] ||  whatd[op_dxy])
 {
  R2 Dl0(K.H(0)), Dl1(K.H(1)), Dl2(K.H(2)), 
     Dlb((Dl0*l1*l2+Dl1*l0*l2+Dl2*l0*l1)*9.);
  
  if (whatd[op_dx]) 
   {
    RN_ f0x(val('.',0,op_dx)); 
   f0x[0] = Dl0.x-Dlb.x;
   f0x[1] = Dl1.x-Dlb.x;
   f0x[2] = Dl2.x-Dlb.x;
   f0x[3] = 3.*Dlb.x;
  }
  
  if (whatd[op_dy]) {
    RN_ f0y(val('.',0,op_dy)); 
   f0y[0] = Dl0.y-Dlb.y;
   f0y[1] = Dl1.y-Dlb.y;
   f0y[2] = Dl2.y-Dlb.y;
   f0y[3] = 3.*Dlb.y;
  }
 if (whatd[op_dxx])
  {  
    RN_ fxx(val('.',0,op_dxx)); 
    R lbdxx= 2*((Dl0.x*Dl1.x)*l2+(Dl1.x*Dl2.x)*l0+(Dl2.x*Dl0.x)*l1);
    fxx[0] = -lbdxx;
    fxx[1] = -lbdxx;
    fxx[2] = -lbdxx;
    fxx[3] = 3*lbdxx;
  }

 if (whatd[op_dyy])
  {  
    RN_ fyy(val('.',0,op_dyy));
    R lbdyy= 2*((Dl0.y*Dl1.y)*l2+(Dl1.y*Dl2.y)*l0+(Dl2.y*Dl0.y)*l1);
     
    fyy[0] =  -lbdyy;
    fyy[1] =  -lbdyy;
    fyy[2] =  -lbdyy;
    fyy[3] =  3*lbdyy;
  }
 if (whatd[op_dxy])
  {  
    assert(val.K()>op_dxy);
    RN_ fxy(val('.',0,op_dxy)); 
    R lbdxy= (Dl0.x*Dl1.y+ Dl0.y*Dl1.x)*l2+(Dl1.x*Dl2.y+Dl1.y*Dl2.x)*l0+(Dl2.x*Dl0.y+Dl2.y*Dl0.x)*l1;  
    fxy[0] = 4*Dl0.x*Dl0.y-9.*(l0-l1-l2);
    fxy[1] = 4*Dl1.x*Dl1.y-9.*(l0-l1-l2);
    fxy[2] = 4*Dl2.x*Dl2.y-9.*(l0-l1-l2);
    fxy[3] = 27.*(l0-l1-l2);
  }
  }

}
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

/* void TypeOfFE_P1Lagrange::FB(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; 
  R2 Dl0(K.H(0)), Dl1(K.H(1)), Dl2(K.H(2));
  
  if (val.N() <3) 
   throwassert(val.N() >=3);
  throwassert(val.M()==1 );
  throwassert(val.K()==3 );
  
  val=0; 
  RN_ f0(val('.',0,0)); 
  RN_ f0x(val('.',0,1)); 
  RN_ f0y(val('.',0,2)); 
  
  f0[0] = l0;
  f0[1] = l1;
  f0[2] = l2;
  
  f0x[0] = Dl0.x;
  f0x[1] = Dl1.x;
  f0x[2] = Dl2.x;
  
  f0y[0] = Dl0.y;
  f0y[1] = Dl1.y;
  f0y[2] = Dl2.y;
}

 void TypeOfFE_P2Lagrange::FB(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; 
  R2 Dl0(K.H(0)), Dl1(K.H(1)), Dl2(K.H(2));
  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; 
  RN_ f0(val('.',0,0)); 
  RN_ f0x(val('.',0,1)); 
  RN_ f0y(val('.',0,2)); 
// --     
  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
  
  
  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) ;
  
  
  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) ;
  
}
*/
void TypeOfFE_P2Lagrange::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; 
  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);