fespace-v0.cpp

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

  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)), Dl1(K.H(1)), Dl2(K.H(2));
  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);
  }
 
 }
 
}
/*
 void TypeOfFE_P1Lagrange::Pi_h(const baseFElement & K,RN_ & val, InterpolFunction f, R* v,int j,  void * arg) const
{
  const R2 Pt[] = { R2(0,0), R2(1,0), R2(0,1) };
   for (int i=0;i<3;i++)
     {  
     f(v,K.T(Pt[i]),K,i,Pt[i],arg),val[i]=*(v+j);}
 
}
 void TypeOfFE_P2Lagrange::Pi_h(const baseFElement & K,RN_ & val, InterpolFunction f, R* v,int j, void * arg) const
{
  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<6;i++)
     {
     f(v,K.T(Pt[i]),K,i,Pt[i],arg),val[i]=*(v+j);}
 
}
*/
 
 
//TypeOfFE  P1Lagrange(1,0,0,P1Functions,D2_P1Functions,P1Interpolant,DataP1Lagrange);
//TypeOfFE  P2Lagrange(1,1,0,P2Functions,D2_P2Functions,P2Interpolant,DataP2Lagrange,3);

//  case of   fine mesh   
class TypeOfMortarCas1: public TypeOfMortar { 
  friend class FESpace;
  friend class FMortar;
  friend class ConstructDataFElement;
  protected:
  int NbLagrangeMult(const Mesh &,const Mortar &M) const ;
 
   int NbDoF(const Mesh &,const Mortar &M,int i) const 
     { int l(M.NbLeft()),r(M.NbRight());
       int n =Max(l,r);
       int mn=Min(l,r);
       return (l+r)*(NbDfOnVertex + NbDfOnEdge) + (n+1)*NbDfOnVertex + n*NbDfOnEdge -mn-1; 
      }
  int NbOfNodes(const Mesh &,const Mortar &M) const // call one time  
     {int l(M.NbLeft()),r(M.NbRight()); return (l+r)*(vertex_is_node+edge_is_node)+1;}
  int NbDoF(const Mesh &,const Mortar &M) const 
     { int l(M.NbLeft()),r(M.NbRight());
       int n =Max(l,r);
       int mn=Min(l,r);
       return (l+r)*(NbDfOnVertex + NbDfOnEdge) + (n+1)*NbDfOnVertex + n*NbDfOnEdge -mn-1; 
      }
  
   int NodeOfDF(const FESpace &Vh,const Mortar &M,int i) const 
     {throwassert(0);return 0;}
   int DFOfNode(const FESpace &Vh,const Mortar &M,int i) const 
     {throwassert(0);return 0;}
   void ConstructionOfNode(const Mesh &Th,int im,int * NodesOfElement,int *FirstNodeOfElement,int &lastnodenumber) const;
   void ConsTheSubMortar(FMortar & ) const; 
     
   const int vertex_is_node,edge_is_node;
  public: 
    TypeOfMortarCas1 (int i,int j): TypeOfMortar(i,j),
      vertex_is_node(i?1:0),edge_is_node(j?1:0) {};
     
}  MortarCas1P2(1,1) ;

 const TypeOfMortar & TheMortarCas1P2(MortarCas1P2); 
 

void TypeOfMortarCas1::ConstructionOfNode(const Mesh &Th,int im,int * NodesOfElement,int *FirstNodeOfElement,int &lastnodenumber) const
{  
  // im   mortar number 
 // trop complique on change 
             const Mortar &M(Th.mortars[im]);
             int k = Th.nt+im;
             int  kk=FirstNodeOfElement[k]; //  begin   
             // lagrange  multiplicator one new node 
              NodesOfElement[kk++] = lastnodenumber++;
/*                               
             int il = M.NbLeft();
             int ir = M.NbRight();
             int ir1 = ir-1;
             //  left
             
             for( int j=0;j<il;j++)  //  numbering vertex0 edge vertex1
              { 
                int K = M.TLeft(j);  // triangle
                int e = M.ELeft(j);  //  edge
                int nbneK = FirstNodeOfElement[K];
                int oe = vertex_is_node ? 3 : 0;
                int i0 = VerticesOfTriangularEdge[e][0];
                int i1 = VerticesOfTriangularEdge[e][1];
                if (vertex_is_node && !j)   //  just the first time 
                   NodesOfElement[kk++]=NodesOfElement[nbneK +i0];
                if (edge_is_node)
                   NodesOfElement[kk++]=NodesOfElement[nbneK+oe+e];
                if (vertex_is_node )  
                   NodesOfElement[kk++]=NodesOfElement[nbneK +i1];
              }

             //  right 
             for( int j=0;j<ir;j++)  //  numbering vertex0 edge vertex1
              { 
                int K = M.TRight(j);  // triangle
                int e = M.ERight(j);  //  edge
                int nbneK = FirstNodeOfElement[K];
                int oe = vertex_is_node ? 3 : 0;
                
                int i0 = VerticesOfTriangularEdge[e][1]; //  change the sens because right side
                int i1 = VerticesOfTriangularEdge[e][0];
              //  if (vertex_is_node &&  !j)   // never 
               //    NodesOfElement[kk++]=NodesOfElement[nbneK +i0];
                if (edge_is_node) 
                   NodesOfElement[kk++]=NodesOfElement[nbneK+oe+e];
                if (vertex_is_node  && (j != ir1) )  //  skip the last 
                   NodesOfElement[kk++]=NodesOfElement[nbneK +i1];
              } */
              
              throwassert(FirstNodeOfElement[k+1]==kk);
}

 R  d1P1F0(const FESpace *,const aSubFMortar *,R x) {return 1-x;}// 1 on 0
 R  d1P1F1 (const FESpace *,const aSubFMortar *,R x) {return x;}//  1 on 1
 
 R  d1P2F0 (const FESpace *,const aSubFMortar *,R x) {return (1-x)*(1-2*x);}// 1 on x=0
 R  d1P2F1(const FESpace *,const aSubFMortar *,R x) {return (1-x)*x*4;} // 1 on x=1/2
 R  d1P2F2(const FESpace *,const aSubFMortar *,R x) {return x*(2*x-1);} // 1 on x=1
 
 void  TypeOfMortarCas1::ConsTheSubMortar(FMortar & sm) const
   { //  constuction of 
   /* 
     int nbsm; // nb of submortar
  aSubFMortar * sm;
  ~FMortar() { delete [] dataDfNumberOFmul; delete [] dataf;}
  private:
  
  int *dataDfNumberOFmul;
   R (**dataf)(const FESpace *,const aSubFMortar *,R);

   */
   //  typedef
     const Mesh &Th(sm.Vh.Th);
     const Mortar & M(sm.M);
     int nl=M.NbLeft();
     int nr=M.NbRight();
     int nbsm= nl+nr-1;
     sm.nbsm = nbsm;
     int ldata = 6*nbsm;// 3 gauche+ 3 droite 
     sm.sm = new aSubFMortar[nbsm];
     sm.datai = new int [ldata];
     sm.dataf = new (R (*[ldata])(const FESpace *,const aSubFMortar *,R))  ;
     int *dataDfNumberOFmul=sm.datai;
     
     R (**dataf)(const FESpace *,const aSubFMortar *,R) ;
     dataf=sm.dataf;
     for (int i=0;i<ldata;i++) sm.dataf[i]=0;
   //  int * data0=sm.data+ldata;
   //  int * data1=data0+ldata;
     
     //  now the construction 
     int l=0,g=0;
     R2 A(M.VLeft(0));
     R2 B(M.VLeft(nl));
     throwassert(&M.VLeft(0) == &M.VRight(0));
     throwassert(&M.VLeft(nl) == &M.VRight(nr));
    
     R2 AB(A,B);
     R lg=Norme2(AB);
    // cout << " Mortar from " << A << " to " << B << " = " <<lg << endl;
     R2 AB1(AB/lg);
     int il=0,ir=0;
     int k=0;
     R la=0.0;
     R lb=0.0;
     R2 AA(A),BB(A);
   //  cout << "lg : " <<lg ;
     do {
       sm.sm[k].left  = M.left[il];
       sm.sm[k].right =  M.right[ir];
       R2 Bl(M.VLeft(il+1));
       R2 Br(M.VRight(ir+1));
       R ll=(AB1,Bl-A), lr=(AB1,Br-A);
     //  throwassert ( ll >=0 && lr >= 0);
     //  throwassert ( ll <=lg  && lr <= lg);
       
   //    cout << "AA , BB = " << AA << "," << BB << endl;
    //   cout << " " << ll << " " << lr << " ll=" << sm.sm[k].left << ", ";
       if (ll<lr) {BB=Bl,lb=ll,il++;} else {BB=Br, lb=lr,ir++;}
  //     cout << k << " " << k << " " << la/lg << " " << lb/lg << endl;
       sm.sm[k].a = la/lg;
       sm.sm[k].b = lb/lg;
       sm.sm[k].A=AA;
       sm.sm[k].B=BB;       
       la=lb;
       AA=BB;
       k++;
       throwassert(k<=nbsm);
     } while (il<nl && ir < nr);
     
  //   cout << "k=" << k <<endl;
     throwassert(nbsm==Max(nl,nr)); 
     //throwassert(nbsm<=k);
    nbsm=k;
    sm.nbsm=k;
//   construction of interpolation 
//  1) on a P1 on P2 
//   P2  si les longueurs des aSubMortar precedende et suivant  sont les meme 
//   sinon P1 
   //  calcul de leps 
     R leps=1.0/(1048576.0) ; // 1/2^20 
     R lgp=0;
     R lgc=0;
     R lgs=sm.sm[0].lg1();
    // cout << lgp << " " << lgc << " " << lgs << endl;
     
     int nmul=0;
     for (int k=0;k<nbsm;k++) 
       {
         
         lgp=lgc;
         lgc=lgs;
         lgs=  k+1 == nbsm  ? 0 : sm.sm[k+1].lg1();
         sm.sm[k].DfNumberOFmul= dataDfNumberOFmul;
         sm.sm[k].f=dataf;
        // cout << lgp << " " << lgc << " " << lgs << " ";
         if ( Abs(lgp-lgc) < leps && Abs(lgs-lgc) < leps )
          { // P2
           sm.sm[k].Nbmul=3;
           *dataDfNumberOFmul++=nmul++;
           *dataDfNumberOFmul++=nmul++;
           *dataDfNumberOFmul++=nmul;
           *dataf++ = d1P2F0;
           *dataf++ = d1P2F1;
           *dataf++ = d1P2F2;
          // cout << "P2 " << nmul << " " ;
           
          }
         else 
          { // P1
                   sm.sm[k].Nbmul=2;
           *dataDfNumberOFmul++=nmul++;
           *dataDfNumberOFmul++=nmul;
           *dataf++ = d1P1F0;
           *dataf++ = d1P1F1;
          // cout << "P1 " << nmul << " " ;
           

          }
       }
      nmul++;
     // cout << " " << nmul << " " <<  sm.NbDoF(0) << endl;
      throwassert(nmul==sm.NbDoF(0));
     
   }
   
  int TypeOfMortarCas1::NbLagrangeMult(const Mesh &,const Mortar &M) const 
     { 
       int nl = M.NbLeft();
       int nr = M.NbRight();
       R2 A(M.VLeft(0));
       R2 B(M.VLeft(nl));
       R2 AB(A,B);
       R lg=Norme2(AB);
       R leps = lg/1048576.0;
       throwassert(nl==1 || nr==1);
       R lgp=0,lgc=0,lgs=0;
       int nbmul=3; 
       if (nr==1) 
        {
        R2 AA(M.VLeft(0)),BB(M.VLeft(1));
        lgp= Norme2(BB-AA); // preced
        AA=BB;
        BB=M.VLeft(2); 
        lgc= Norme2(BB-AA); // courant 
        
        for (int i=1;i<nl-1;i++)
         { 
            AA=BB;
            BB=M.VLeft(i+2); 
            lgs=Norme2(AA-BB); // le suivant 
            if ( Abs(lgp-lgc) < leps && Abs(lgs-lgc) < leps )
              nbmul+=2; // P2
            else 
              nbmul+=1;// P1;
            lgp=lgc;
            lgc=lgs;
            
         }
        }
        else
        {
        R2 AA(M.VRight(0)),BB(M.VRight(1));
        lgp= Norme2(BB-AA); // preced
        AA=BB;
        BB=M.VRight(2); 
        lgc= Norme2(BB-AA); // courant 
        
        for (int i=1;i<nr-1;i++)
         { 
            AA=BB;
            BB=M.VRight(i+2); 
            lgs=Norme2(AA-BB); // le suivant 
            if ( Abs(lgp-lgc) < leps && Abs(lgs-lgc) < leps )
              nbmul+=2; // P2
            else 
              nbmul+=1;// P1;
            lgp=lgc;
            lgc=lgs;
            
         }
        }
       throwassert(nbmul>2);
       return nbmul;
      }  

 
// --- 
 FMortar::FMortar(const FESpace * VVh,int k)
  : 
    Vh(*VVh),
    M(Vh.Th.mortars[k-Vh.Th.nt]),
    N(VVh->N),
    p(Vh.PtrFirstNodeOfElement(k)),
    nbn(Vh.NbOfNodesInElement(k)),
    tom(Vh.tom)
    
 { throwassert(k>=Vh.Th.nt && k <Vh.Th.nt + Vh.Th.NbMortars);
   VVh->tom->ConsTheSubMortar(*this);}
   
 R TypeOfFE::operator()(const FElement & K,const  R2 & PHat,const KN_<R> & u,int componante,int op) const 
  {
   R v[1000],vf[100];
   assert(N*3*NbDoF<=1000 && NbDoF <100 );
   KNMK_<R> fb(v,NbDoF,N,op+1); //  the value for basic fonction
   KN_<R> fk(vf,NbDoF);
   for (int i=0;i<NbDoF;i++) // get the local value
    fk[i] = u[K(i)];
    //  get value of basic function
   bool whatd[last_operatortype];
   for (int i=0;i<last_operatortype;i++) 
     whatd[i]=false;
   whatd[op]=true;
   FB(whatd,K.Vh.Th,K.T,PHat,fb);  
   R r = (fb('.',componante,op),fk);  
   return r;
  }
 
static TypeOfFE_P1Lagrange P1LagrangeP1;
static TypeOfFE_P1Bubble P1BubbleP1;
static TypeOfFE_P2Lagrange P2LagrangeP2;

TypeOfFE  & P2Lagrange(P2LagrangeP2);
TypeOfFE  & P1Bubble(P1BubbleP1);
TypeOfFE  & P1Lagrange(P1LagrangeP1);

static ListOfTFE typefemP1("P1", &P1LagrangeP1);
static ListOfTFE typefemP1b("P1b", &P1BubbleP1);
static ListOfTFE typefemP2("P2", &P2LagrangeP2);
static  ListOfTFE typefemRT("RT0", &RTLagrange);
static  ListOfTFE typefemRTOrtho("RT0Ortho", &RTLagrangeOrtho);
 
 extern  TypeOfFE & RTmodifLagrange, & P1ttdc, & P2ttdc;
 static  ListOfTFE typefemRTmodif("RTmodif", &RTmodifLagrange);
 static ListOfTFE typefemP0("P0", &P0Lagrange);
 static ListOfTFE typefemP1nc("P1nc", &P1ncLagrange);
 static ListOfTFE typefemP1ttdc("P1dc", &P1ttdc);
 static ListOfTFE typefemP2ttdc("P2dc", &P2ttdc);
 
} // fin de namespace Fem2D