lgfem.cpp

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

    fe=y;
    if(verbosity>1)
    cout << " -- interpole f= g*X^-1, function's bound:  " << y->min() << " " << y->max() << endl; 
    return SetAny<FEbase<R,v_fes>*>(&fe); 
  }

template<class K>
E_set_fev<K>::E_set_fev(const E_Array * a,Expression pp,int ddim) 
  :dim(ddim), aa(*a),ppfe(pp),optimize(true),
   where_in_stack_opt(),optiexp0(),optiexpK() 
{ 
  aa.map(to<K>) ;
  bool kdump=false;
  if(optimize)
    { // new code Optimized  -------
      int n=aa.size();
      deque<pair<Expression,int> > ll;
      MapOfE_F0 m;
      where_in_stack_opt.resize(n);
      size_t top = currentblock->OffSet(0), topbb=top; // FH. bofbof ??? 
      for (int i=0; i<n; i++)
	{
	  Expression ee= aa[i].LeftValue();
	  if (kdump)
	    cout << "Optimize OneOperatorMakePtrFE:  type exp: " << typeid(*ee).name() << " "<<endl;
	  where_in_stack_opt[i]=ee->Optimize(ll, m, top);
	  if (kdump)
	    cout  << "\n\t\t"<< i  << ": " << where_in_stack_opt[i] << endl;
	}
      
      currentblock->OffSet(top-topbb);
      //  
      int k=ll.size(),k0=0,k1=0;
      for (int i=0;i<k;i++)
	if (ll[i].first->MeshIndependent()) k0++;
      deque<pair<Expression,int> > l0(k0),l1(k-k0);
      k0=0,k1=0;
      for (int i=0;i<k;i++)
	if (ll[i].first->MeshIndependent()) 
	  {
	    if (kdump)
	      cout << " mi " << ll[i].second << " " << *(ll[i].first) << endl;
	    l0[k0++]=ll[i];
	  }
	else 
	  {
	    if (kdump)
	      cout << " md " << ll[i].second << " " << *(ll[i].first) << endl;
	    l1[k1++]=ll[i];
	  }
      if (k0)      
	optiexp0 = new E_F0_Optimize(l0,m,0);  // constant part
      if (k1) 
	optiexpK = new E_F0_Optimize(l1,m,0);  // none constant part
      
    }
  
  
}


template<class K>   
AnyType E_set_fev<K>::operator()(Stack s)  const
{
  if(dim== 2)  return  Op2d(s);
  else if(dim == 3)  return Op3d(s);
}


template<class K>   
AnyType E_set_fev<K>::Op3d(Stack s)  const
{
  //  voir E_set_fev3  ( pb de consitance a revoir FH) 
  ffassert(0); // a faire   
}
template<class K>   
AnyType E_set_fev<K>::Op2d(Stack s)  const
{  
  StackOfPtr2Free * sptr = WhereStackOfPtr2Free(s);     
  MeshPoint *mps=MeshPointStack(s), mp=*mps;   
  FEbase<K,v_fes> ** pp=GetAny< FEbase<K,v_fes> **>((*ppfe)(s));
  FEbase<K,v_fes> & fe(**pp);
  const  FESpace & Vh(*fe.newVh());
  KN<K> gg(Vh.MaximalNbOfDF()); 
  
  
  const  Mesh & Th(Vh.Th);
  const int dim=Vh.N;
  K ** copt=0;
  if (optimize)   copt= new K *[dim];
  if(copt) {
    assert((size_t) dim== where_in_stack_opt.size());
    for (int i=0;i<dim;i++)
      {
        int offset=where_in_stack_opt[i];
        assert(offset>10);
        copt[i]= static_cast<K *>(static_cast<void *>((char*)s+offset));
        *(copt[i])=0;
       }
    if (optiexp0) (*optiexp0)(s); // init 
  }
  
  ffassert(dim<100);
  //   R F[100]; // buffer 
  
  TabFuncArg tabexp(s,Vh.N);
  //   const E_Array * aa = dynamic_cast<const E_Array *>(e);
  ffassert( aa.size() == Vh.N);
  for (int i=0;i<dim;i++)
    tabexp[i]=aa[i]; 
  
  KN<K> * y=new  KN<K>(Vh.NbOfDF);
  KN<K> & yy(*y);
  
  FElement::aIPJ ipj(Vh[0].Pi_h_ipj()); 
  FElement::aR2  PtHat(Vh[0].Pi_h_R2()); 
  
  KN<double>   Aipj(ipj.N());
  
  
  //   KNM<K>  Vp(dim,PtHat.N());// bug 
  // g++ error: too many initializers for `const __class_type_info_pseudo
  
  KN<K>  Vp1(dim*PtHat.N());
  
  
  if (Vh.isFEMesh() )
    {
      
      ffassert(Vh.NbOfDF == Th.nv && dim == 1 );
      for (int iv=0;iv<Th.nv;iv++)
	{
	  const E_F0 & ff(* (const  E_F0 *) aa[0]  ) ;
	  const Vertex & v(Th(iv));
	  int ik=Th.Contening(&v);
	  const Triangle & Kt(Th[ik]);
	  int il=-1;
	  if  ( &Kt[0] == &v) il=0;
	  if  ( &Kt[1] == &v) il=1;
	  if  ( &Kt[2] == &v) il=2;
	  assert(il>=0);
	  mps->set(Th,v,TriangleHat[il],Kt,v.lab);
	  if (copt) {
	    if (optiexpK) (*optiexpK)(s); 
            yy[iv] =  *(copt[0]);
          }
          else 
	    yy[iv] = GetAny<K>( ff(s) );
	  sptr->clean(); // modif FH mars 2006  clean Ptr
	}
      
    }
  else
    for (int t=0;t<Th.nt;t++)
      {
	FElement Kt(Vh[t]);
	int nbdf=Kt.NbDoF();
        
	gg=K();
        
#ifdef OLDPih    
	// old method          
	Kt.Pi_h(gg,F_Pi_h,F,&tabexp);
#else               
	Kt.Pi_h(Aipj);
	
	for (int p=0;p<PtHat.N();p++)
          { 
            mps->set(Kt.T(PtHat[p]),PtHat[p],Kt);

      //      KN_<K> Vpp(Vp('.',p));
            KN_<K> Vpp(Vp1,SubArray(dim,p*dim)); // a Change FHHHHHHHH
            if (copt) { // optimize  version 
             if (optiexpK) (*optiexpK)(s);
             for (int j=0;j<dim;j++)
               Vpp[j] = *(copt[j]);
            }
            else  // old version 
            for (int j=0;j<dim;j++)
             if (tabexp[j]) 
               Vpp[j]=GetAny<K>( (*tabexp[j])(s) );
              else Vpp[j]=0;
              
           }
           
         for (int i=0;i<Aipj.N();i++)
          { 
           const FElement::IPJ &ipj_i(ipj[i]);
         //  gg[ipj_i.i] += Aipj[i]*Vp(ipj_i.j,ipj_i.p);           
             gg[ipj_i.i] += Aipj[i]*Vp1(ipj_i.j+ipj_i.p*dim); // index a la main  
            sptr->clean(); // modif FH mars 2006  clean Ptr          
          } 
#endif          

         for (int df=0;df<nbdf;df++)         
           yy[Kt(df)] =  gg[df] ;
      }
   //  MeshPointStack(s)->unset();
    fe=y;
    if (copt) delete [] copt;
    *MeshPointStack(s) = mp;
    if(verbosity>1)
        ShowBound(*y,cout) << endl ;
   //HHHH*/       
    return Nothing;
  }


template<class K>
inline FEbase<K,v_fes> * MakePtrFE(pfes * const &  a){ 
  FEbase<K,v_fes> * p=new FEbase<K,v_fes>(a);
  //cout << "MakePtrFE " << p<< endl; 
  return p ;}
  
template<class K>
inline FEbase<K,v_fes> ** MakePtrFE2(FEbase<K,v_fes> * * const &  p,pfes * const &  a){ 
  *p=new FEbase<K,v_fes>(a);
  //cout << "MakePtrFE2 " << *p<< endl; 
  return p ;}

template<class K>  
inline FEbaseArray<K,v_fes> ** MakePtrFE3(FEbaseArray<K,v_fes> * * const &  p,pfes * const &  a,const long & N){ 
  *p=new FEbaseArray<K,v_fes>(a,N);
  //cout << "MakePtrFE2 " << *p<< endl; 
  return p ;}
  
/*
inline pmesharray*  MakePtr(pmesharray*  const &  p,long   const &  a){ 
  p->first=new pmesh [a];
  p->second=a;
  for (int i=0;i<a;i++) 
     p->first[i]=0; // nuset 
  return p ;}
*/  
template<class K>
class  OneOperatorMakePtrFE : public OneOperator 
{
public:
  // il faut Optimize 
  // typedef double K;
  typedef  FEbase<K,v_fes> ** R;
  typedef pfes* B;
  class CODE : public E_F0mps  
  {
  public:
    Expression fer,fes;
    E_set_fev<K> * e_set_fev;
    const E_Array * v;
    CODE(const basicAC_F0 & args) 
      : 
      fer(to<R>(args[0])),
      fes(to<B>(args[1])),
      e_set_fev(0) 
    {
      if (BCastTo<K>(args[2]) )
	v = new E_Array(basicAC_F0_wa(to<K>(args[2]))); 
      else 
	v = dynamic_cast<const E_Array *>( args[2].LeftValue() );
      if (!v) {
	cout << "Error: type of arg :" << *args[2].left()  << " in " << typeid(K).name() << " case " << endl;
	ErrorCompile(" We wait  a double/complex expression or a array expression",1);
      }
      //v->map(to<K>);
      e_set_fev=  new   E_set_fev<K>(v,fer,v_fes::d);
      
    }
    
    AnyType operator()(Stack stack)  const {
      R  p = GetAny<R>( (*fer)(stack));
      B  a = GetAny<B>( (*fes)(stack)); 
      *p=new FEbase<K,v_fes>(a);
      (*e_set_fev)(stack); 
      //  cout << "MakePtrFE: build p " << p << " " << *p << endl;         
      return SetAny<R>(p);
    }
    operator aType () const { return atype<R>();}         
    
    
  };
  
  E_F0 * code(const basicAC_F0 & args) const 
  { return  new CODE(args);}
  OneOperatorMakePtrFE(aType tt):  // tt= aType<double>() or aType<E_Array>()  
    OneOperator(map_type[typeid(R).name()],map_type[typeid(R).name()],map_type[typeid(B).name()],tt)
  {}
};


// ---  

template<class Result,class A>
class  OneOperator_Ptr_o_R: public OneOperator {
  //  aType r; //  return type 
    typedef Result A::* ptr;
    ptr p;   
    public: 
    E_F0 * code(const basicAC_F0 & args) const 
     { return  new E_F_A_Ptr_o_R<Result,A>(t[0]->CastTo(args[0]),p);} 
    OneOperator_Ptr_o_R(ptr pp):
       OneOperator(atype<Result*>(),atype<A*>()),p(pp) {}     
};

template<class K>  K  *PAddition(const K * a,const K *  b)  {return  new K(*a+*b);}





// 
class fCLD { public:
  typedef pair<int,Label> Key;
  typedef  map<Key,Expression>::iterator iterator;
  map<Key,Expression> *l;
  fCLD (){ l=new  map<Key,Expression>;}
  void operator=(const fCLD & a){ *l=*a.l;}
  void destroy() { delete l;l=0;}
  ~fCLD(){ delete l;l=0;}
  void Add(finconnue *v,int lab,C_F0 f) { 
    const MGauche *pn= v->simple();
     Check(pn,"Def CL Dirichet ");
    
    Label r(lab);
    Key k(make_pair(pn->first,r));
    iterator i=l->find(k);
    Check( i != l->end() ,"Def CL Dirichet already exists");
    l->insert(make_pair(k,CastTo<double>(f)));        
  }
};

//  pour stocker des expression de compilation 


class Convect : public E_F0mps  { public:
    typedef double  Result; // return type 
    Expression u,v,w,ff,dt;
    int d;
    Convect(const basicAC_F0 & args)  : u(0),v(0),w(0),ff(0),dt(0)
    {
      args.SetNameParam(); 
      const E_Array * a = dynamic_cast<const E_Array *>(args[0].LeftValue());
            ffassert(a);
	d= a->size();
       if (d == 3) 
	   w= CastTo<double>((*a)[2]);
       else if (d != 2) 
          { CompileError("convect vector have only 2 or 3 componant");}
       u= CastTo<double>((*a)[0]);
       v= CastTo<double>((*a)[1]);
       
       dt=CastTo<double>(args[1]);
       ff=CastTo<double>(args[2]);
     }
    
    static ArrayOfaType  typeargs() 
      { return  ArrayOfaType(atype<E_Array>(),atype<double>(),atype<double>());}
      
    static  E_F0 * f(const basicAC_F0 & args) { return new Convect(args);} 
    AnyType operator()(Stack s) const ; 
    AnyType eval2(Stack s) const ; 
    AnyType eval3(Stack s) const ; 
    operator aType () const { return atype<Result>();}         
    
};

class Plot :  public E_F0mps { public:
    typedef KN_<R>  tab;
    typedef pferbase sol;
    typedef pf3rbase sol3;
    
    typedef long  Result;
    struct Expression2 {
     long what; // 0 mesh, 1 iso, 2 vector, 3 curve , 4 border , 5  mesh3, 6 iso 3d 
     bool composant;
     Expression e[2];
     Expression2() {e[0]=0;e[1]=0;composant=false;what=0;}
     Expression &operator[](int i){return e[i];}
     sol eval(int i,Stack s,int & cmp) const  {  cmp=-1;
       if (e[i]) {
        if (!composant) {pfer p= GetAny< pfer >((*e[i])(s)); cmp=p.second;return p.first;}
        else {return GetAny< pferbase >((*e[i])(s));}
        }
       else return 0;}
      sol3 eval3(int i,Stack s,int & cmp) const  {  cmp=-1;
	    if (e[i]) {
		if (!composant) {pf3r p= GetAny< pf3r >((*e[i])(s)); cmp=p.second;return p.first;}
		else {return GetAny< pf3rbase >((*e[i])(s));}
	    }
	else return 0;}
	const Mesh & evalm(int i,Stack s) const  { throwassert(e[i]);return  * GetAny< pmesh >((*e[i])(s)) ;}
     const Mesh3 & evalm3(int i,Stack s) const  { throwassert(e[i]);return  * GetAny< pmesh3 >((*e[i])(s)) ;}
     const E_BorderN * evalb(int i,Stack s) const  { throwassert(e[i]);return   GetAny< const E_BorderN *>((*e[i])(s)) ;}
     tab  evalt(int i,Stack s) const  { throwassert(e[i]);return  GetAny<tab>((*e[i])(s)) ;}
    };

   static basicAC_F0::name_and_type name_param[] ;
   static const int n_name_param =20 ;
   Expression bb[4];
    vector<Expression2> l;
    Expression nargs[n_name_param];
    Plot(const basicAC_F0 & args) : l(args.size()) 
    {

      args.SetNameParam(n_name_param,name_param,nargs);
       if ( nargs[8] )
         Box2x2( nargs[8] , bb);   
         
      for (size_t i=0;i<l.size();i++)
       
         if (args[i].left()==atype<E_Array>())
          {
	    cout << "args[i].left()==atype<E_Array>()" << endl;
	    l[i].composant=false;
            const E_Array * a = dynamic_cast<const E_Array *>(args[i].LeftValue());
            ffassert(a);
            if (a->size() >2) { CompileError("plot d'un vecteur