lgfem.cpp

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

class OP_MakePtr3 { public:
    class Op : public E_F0mps  { public:
	//  static int GetPeriodic(Expression  bb, Expression & b,Expression & f);
	static const int n_name_param =1;
      static basicAC_F0::name_and_type name_param[] ;
      Expression nargs[n_name_param];
      typedef pfes3 * R;
      typedef pfes3 * A;
      typedef pmesh3 * B;
      typedef TypeOfFE3 * C;
      Expression a,b,c;
      int nbcperiodic ;
      Expression *periodic;
      Op(const basicAC_F0 & args);
      
      AnyType operator()(Stack s) const  {
	A p= GetAny<A>( (*a)(s) );
	B th= GetAny<B>( (*b)(s) );
	C tef= GetAny<C>( (*c)(s) );   
	//  cout << " ----------- " << endl;  
	*p=new pfes3_tef(th,tef,s,nbcperiodic,periodic) ;
	(**p).decrement();
	return  SetAny<R>(p);
      } 
    }; // end Op class 
  
  typedef Op::R Result;
  static  E_F0 * f(const basicAC_F0 & args) { return  new Op(args);}
  static ArrayOfaType  typeargs() {
    return ArrayOfaType(       
			atype<Op::A>(),
			atype<Op::B>(),
			atype<Op::C>(),false ) ;}
};  

void GetPeriodic(Expression perio,    int & nbcperiodic ,    Expression * &periodic)
{
      if ( perio) 
       {
         if( verbosity>1) 
         cout << "  -- Periodical Condition to do" << endl;
         const E_Array * a= dynamic_cast<const  E_Array *>(perio);
         ffassert(a);
         int n = a->size();
        nbcperiodic= n/2;
        if( verbosity>1) 
        cout << "    the number of periodicBC " << n << endl;
        if ( 2*nbcperiodic != n ) CompileError(" Sorry the number of periodicBC must by even"); 
        periodic = new Expression[n*2]; 
        for (int i=0,j=0;i<n;i++,j+=2)
          if (GetPeriodic((*a)[i],periodic[j],periodic[j+1])==0)
            CompileError(" a sub array of periodic BC must be [label, realfunction ]");
        }


}

    
OP_MakePtr2::Op::Op(const basicAC_F0 & args)
  : a(to<A>(args[0])),b(to<B>(args[1])),c(to<C>(args[2]))
     {
      nbcperiodic=0;
      periodic=0;
      args.SetNameParam(n_name_param,name_param,nargs);
      GetPeriodic(nargs[0],nbcperiodic,periodic);
     }
    
OP_MakePtr3::Op::Op(const basicAC_F0 & args)
  : a(to<A>(args[0])),b(to<B>(args[1])),c(to<C>(args[2]))
{
  nbcperiodic=0;
  periodic=0;
  args.SetNameParam(n_name_param,name_param,nargs);
  GetPeriodic(nargs[0],nbcperiodic,periodic);
}

int GetPeriodic(Expression  bb, Expression & b,Expression & f)
    {
      const E_Array * a= dynamic_cast<const E_Array *>(bb);
      if(a && a->size() == 2)
	{
	  b= to<long>((*a)[0]);
	  f= to<double>((*a)[1]);
	  return 1;
	}
      else 
	return 0;
    }

basicAC_F0::name_and_type  OP_MakePtr2::Op::name_param[]= {
  "periodic", &typeid(E_Array) 
};

basicAC_F0::name_and_type  OP_MakePtr3::Op::name_param[]= {
  "periodic", &typeid(E_Array) 
};

inline pfes* MakePtr2(pfes * const &p,pmesh * const &  a){ 
      *p=new pfes_tef(a,&P1Lagrange);
      (**p).decrement();
       return p ;}
       
inline pfes* MakePtr2(pfes * const &p,pfes * const &  a,long const & n){
       *p= new pfes_fes(a,n);
      (**p).decrement();
       return p ;}
       
 long FindTxy(Stack s,pmesh * const &  ppTh,const double & x,const double & y)
{
   R2 P(x,y),PHat;
   bool outside;
     MeshPoint & mp = *MeshPointStack(s);  
   const Mesh * pTh= *ppTh;
   if(pTh == 0) return 0;  
   const Triangle * K=pTh->Find(mp.P.p2(),PHat,outside);
   if (!outside)
     mp.set(*pTh,P,PHat,*K,K->lab);
   else return 0;
   return 1;
}

 


 

template<class K>    
KN<K> * pfer2vect( pair<FEbase<K,v_fes> *,int> p)
 {  
    KN<K> * x=p.first->x();
    if ( !x) {  // defined 
      FESpace * Vh= p.first->newVh();     
      throwassert( Vh);
      *p.first = x = new KN<K>(Vh->NbOfDF);
      *x=K(); 
    }
    return x;}

template<class K>        
long pfer_nbdf(pair<FEbase<K,v_fes> *,int> p)
 {  
   if (!p.first->Vh) p.first->Vh= p.first->newVh();
   throwassert( !!p.first->Vh);
   return p.first->Vh->NbOfDF;
 }
 
double pmesh_area(pmesh * p)
 { throwassert(p && *p) ;  return (**p).area ;}
long pmesh_nt(pmesh * p)
 { throwassert(p && *p) ;  return (**p).nt ;}
long pmesh_nv(pmesh * p)
 { throwassert(p && *p) ;  return (**p).nv ;}
long pVh_ndof(pfes * p)
 { throwassert(p && *p);
   FESpace *fes=**p; ;  return fes->NbOfDF ;}
long pVh_nt(pfes * p)
 { throwassert(p && *p);
   FESpace *fes=**p; ;  return fes->NbOfElements ;}
long pVh_ndofK(pfes * p)
 { throwassert(p && *p);
   FESpace *fes=**p;   return (*fes)[0].NbDoF() ;}

long mp_nuTriangle(MeshPoint * p)
 { throwassert(p  && p->Th && p->T);
   long nu=0;
   if(p->d==2)
     nu=(*p->Th)(p->T);
   else if  (p->d==3)
     nu=(*p->Th3)(p->T3);
   else ffassert(0);
   delete p;
   return nu ;}
   
long mp_region(MeshPoint * p)
 { throwassert(p && p->Th);
   long  nu(p->region);
   delete p;
   return nu ;}


class pVh_ndf : public ternary_function<pfes *,long,long,long> { public:


  class Op : public E_F0mps { public:
      Expression a,b,c;
       Op(Expression aa,Expression bb,Expression cc) : a(aa),b(bb),c(cc) {}       
       AnyType operator()(Stack s)  const 
        { 
           pfes * p(GetAny<pfes *>((*a)(s)));
           long  k(GetAny<long>((*b)(s)));
           long  i(GetAny<long>((*c)(s)));
           throwassert(p && *p);
           FESpace *fes=**p;
           throwassert(fes && k >=0 && k < fes->NbOfElements );
           FElement K=(*fes)[k];
           throwassert(i>=0 && i <K.NbDoF() );
           long ret(K(i));
           return  ret;
         }
   
  };
};
//plus

class Op_CopyArray : public OneOperator { public:
    Op_CopyArray():OneOperator(atype<void>(),atype<E_Array>(),atype<E_Array>()) {}
    E_F0 * code(const basicAC_F0 & args) const ; 
};  
                       
template<class R,int dd>
AnyType pfer2R(Stack s,const AnyType &a)
{
  pair< FEbase<R,v_fes> *  ,int> ppfe=GetAny<pair< FEbase<R,v_fes> *,int> >(a);
  FEbase<R,v_fes> & fe( *ppfe.first);
  int componante=ppfe.second;
  if ( !fe.x()) {
   if ( !fe.x()){
    // CompileError(" Sorry unset fem array ");
     return   SetAny<R>(0.0);
    }
  }
   
  const FESpace & Vh(*fe.Vh);
  const Mesh & Th(Vh.Th);
  assert(Th.ntet==0 && Th.volume==0 && Th.triangles != 0);
  MeshPoint & mp = *MeshPointStack(s);
  const Triangle *K;
  R2 PHat;
  bool outside=false;
  bool qnu=true;
  if ( mp.Th == &Th && mp.T) 
   {
    qnu=false;
    K=mp.T;
    PHat=mp.PHat.p2();
   }
  else if ( mp.other.Th == & Th && mp.other.P.x == mp.P.x && mp.other.P.y == mp.P.y )
   {
    K=mp.other.T;
    PHat=mp.other.PHat.p2();
    outside = mp.other.outside;
   } 
  else {
    if (mp.isUnset()) ExecError("Try to get unset x,y, ...");
    K=Th.Find(mp.P.p2(),PHat,outside);
    mp.other.set(Th,mp.P.p2(),PHat,*K,0,outside);
    }
  // cout << "  ---  " << qnu << "  " << mp.P << " " << mp.outside <<  " " << outside << endl;
  const FElement KK(Vh[Th(K)]);
  if (outside && !KK.tfe->NbDfOnVertex && !KK.tfe->NbDfOnEdge) 
    return   SetAny<R>(0.0); 
/*  if (!outside) 
    {
      if ( Norme2_2( (*K)(PHat) - mp.P ) > 1e-12 )
        cout << "bug ??  " << Norme2_2( (*K)(PHat) - mp.P ) << " " << mp.P << " " << (*K)(PHat) << endl;
    } */
/*  int nbdf=KK.NbDoF();
  
  int N= KK.N;
  KN_<R> U(*fe.x());
  KNMK<R> fb(nbdf,N,3); //  the value for basic fonction
  KN<R> fk(nbdf);
  for (int i=0;i<nbdf;i++) // get the local value
    fk[i] = U[KK(i)];
    //  get value of basic function
  KK.BF(PHat,fb);  
  
  R r = (fb('.',componante,dd),fk);  
*/
  
 const R rr = KK(PHat,*fe.x(),componante,dd);
//  cout << " " << rr << endl;
//  R2 B(mp.P);
/*   if ( r < 0.08001 &&  Norme2_2(mp.P) > 0.05  && Norme2_2(mp.P) < 0.4*0.4  ) 
     {
     int vv=verbosity;
      cout << " f()  triangle  " << Th(K) << " " << mp.P << " " << PHat << " =  " << r << " " <<outside ;
      if (outside) {  verbosity = 200;
         K=Th.Find(mp.P,PHat,outside);
          cout << Th(K) << " " << outside << endl;}
       cout << endl; verbosity=vv;
     } */
//  if ( qnu )   
 //  cout << " f()  triangle       " << Th(K) << " " << mp.P << " " << PHat << " =  " << r <<  endl;
  return SetAny<R>(rr);
}

 
template<class R>
AnyType set_fe (Stack s,Expression ppfe, Expression e)
  { 
   long kkff = Mesh::kfind,  kkth = Mesh::kthrough;
      StackOfPtr2Free * sptr = WhereStackOfPtr2Free(s);

   
    MeshPoint *mps=MeshPointStack(s),mp=*mps;  
    pair<FEbase<R,v_fes> *,int>  pp=GetAny<pair<FEbase<R,v_fes> *,int> >((*ppfe)(s));
    FEbase<R,v_fes> & fe(*pp.first);
    const  FESpace & Vh(*fe.newVh());
    KN<R> gg(Vh.MaximalNbOfDF()); 
    const  Mesh & Th(Vh.Th);
 //   R F[100]; // buffer 
    TabFuncArg tabexp(s,Vh.N);
    tabexp[0]=e;
    
    if(Vh.N!=1)
    {  cerr << " Try to set a  vectorial  FE function  (nb  componant=" <<  Vh.N << ") with one scalar " << endl;
       ExecError(" Error interploation (set)  FE function (vectorial) with a scalar");
    }
    KN<R> * y=new  KN<R>(Vh.NbOfDF);
    KN<R> & yy(*y);
    KN<R> Viso(100);
    // R2 Ptt[3];
    for (int i=0;i<Viso.N();i++)
      Viso[i]=0.01*i; 
      
    FElement::aIPJ ipj(Vh[0].Pi_h_ipj()); 
    FElement::aR2  PtHat(Vh[0].Pi_h_R2()); 
    KN<double>   Aipj(ipj.N());
    KN<R>   Vp(PtHat.N());
    const E_F0 & ff(* (const  E_F0 *) e ) ;

   if (Vh.isFEMesh() )
    {
      
      ffassert(Vh.NbOfDF == Th.nv && Vh.N == 1 );
      for (int iv=0;iv<Th.nv;iv++)
       {
         const Vertex & v(Th(iv));
         int ik=Th.Contening(&v);
         const Triangle & K(Th[ik]);
         int il=-1;
         if  ( &K[0] == &v) il=0;
         if  ( &K[1] == &v) il=1;
         if  ( &K[2] == &v) il=2;
         assert(il>=0);
         mps->set(Th,v,TriangleHat[il],K,v.lab);
         yy[iv] = GetAny<R>( ff(s) );
          sptr->clean(); // modif FH mars 2006  clean Ptr
       }
      
    }
   else     
    
     for (int t=0;t<Th.nt;t++)
      {
         FElement K(Vh[t]);
         int nbdf=K.NbDoF();
         gg=R();   
#ifdef OLDPih    
// old method          
         K.Pi_h(gg,F_Pi_h,F,&tabexp);
#else               
         K.Pi_h(Aipj);
         for (int p=0;p<PtHat.N();p++)
          { 
            mps->set(K.T(PtHat[p]),PtHat[p],K);
            Vp[p]=GetAny<R>( ff(s) );
           }
         for (int i=0;i<Aipj.N();i++)
          { 
           const FElement::IPJ &ipj_i(ipj[i]);
           assert(ipj_i.j==0); // car  Vh.N=0
           gg[ipj_i.i] += Aipj[i]*Vp[ipj_i.p];            
          }
#endif          
         for (int df=0;df<nbdf;df++)         
            (*y)[K(df)] =  gg[df] ;
           sptr->clean(); // modif FH mars 2006  clean Ptr
                       
      }
    *mps=mp;
    fe=y;
     kkff = Mesh::kfind - kkff;
     kkth = Mesh::kthrough -kkth;

    if(verbosity>1)
      ShowBound(*y,cout)  
           << " " << kkth << "/" << kkff << " =  " << double(kkth)/Max<double>(1.,kkff) << endl;
    return SetAny<FEbase<R,v_fes>*>(&fe); 
  }
AnyType set_feoX_1 (Stack s,Expression ppfeX_1, Expression e)
  { // inutile 
    // m恗e chose que  v(X1,X2);
      StackOfPtr2Free * sptr = WhereStackOfPtr2Free(s);
    typedef const interpolate_f_X_1<R>::CODE * code;
    MeshPoint mp= *MeshPointStack(s); 
    code ipp = dynamic_cast<code>(ppfeX_1);
    
    pair<FEbase<R,v_fes> *,int>  pp=GetAny<pair<FEbase<R,v_fes> *,int> >((*ipp->f)(s));
    FEbase<R,v_fes> & fe(*pp.first);
    const  FESpace & Vh(*fe.newVh());
    KN<R> gg(Vh.MaximalNbOfDF()); 
    const  Mesh & Th(Vh.Th);
    R F[100]; // buffer 
    TabFuncArg tabexp(s,Vh.N+2);
    tabexp[0]=e;
    tabexp[1]=ipp->x;
    tabexp[2]=ipp->y;
    
    throwassert(Vh.N==1); 
    KN<R> * y=new  KN<R>(Vh.NbOfDF); 
     for (int t=0;t<Th.nt;t++)
      {
         FElement K(Vh[t]);
         int nbdf=K.NbDoF();
        
         gg=R();
         
         K.Pi_h(gg,FoX_1_Pi_h,F,&tabexp);
         for (int df=0;df<nbdf;df++)
          (*y)[K(df)] =  gg[df] ;
          sptr->clean(); // modif FH mars 2006  clean Ptr          
      }
    *MeshPointStack(s)=mp;