array_tlp.hpp

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

	         M = lli.size(); ffassert( M>0 );
	        for (int i=0;i<N;i++) tab[i] = new Expression [M];
	       }
	     else {  
	        if ( M != li->size() ) { 
	        cout << " line " << i << " the size of the column change " << M << " to " << li->size() << endl;
	        CompileError(" Is not a matrix,  M is not constant" ); } } 
	        
	    for (int j=0;j<M;j++)
              tab[i][j]=atype<RR>()->CastTo(  lli[j]);
	   }
	 else  // li == 0 
	  CompileError(" we are waiting for  vector of scalar [  , , ,  ] ");
	 }
	 
    }
    
    AnyType operator()(Stack stack)  const 
    {
      A  a=GetAny<A>((*a0)(stack));
      if (isinit) 
        a->init(N,M);
      else
	a->resize(N,M);
      
       for (int i =0;i<N;++i)
       for (int j =0;j<M;++j)
          (*a)(i,j)=   GetAny< RR >( (*(tab[i][j]))(stack)) ; 
      return SetAny<R>(a);
    } 
    bool MeshIndependent() const     {return  mi;} // 
    ~CODE() { for (int i=0;i<N;i++) delete [] tab[i]; delete [] tab; }
    operator aType () const { return atype<R>();}    
  }; // end sub class CODE
  
  
    public: 
    E_F0 * code(const basicAC_F0 & args) const 
     { return  new CODE(t[0]->CastTo(args[0]),*dynamic_cast<const E_Array*>( t[1]->CastTo(args[1]).LeftValue()));} 
    InitMatfromAArray():   OneOperator(atype<R>(),atype<A>(),atype<B>())  {}
  
};

template<typename RR>
class  SetArrayofKNfromKN : public OneOperator { 
public:
    typedef KN_<RR>  A; // Warning  B type of  1 parameter 
    typedef KN_<RR>  R;
    typedef E_Array B; //   A type of 2 parameter
    
    class CODE : public  E_F0 { public:
       Expression a0;
       int N;
       Expression * tab;
       int * what;//  0  RR, 1 KN<RR>, 
       const  bool mi;

    CODE(Expression a,const E_Array & tt)  
      : a0(a),N(tt.size()),
	tab(new Expression [N]),
	what(new int[N])  ,
	mi(tt.MeshIndependent())
      {
        assert(&tt);
	//      int err=0;
        for (int i=0;i<N;i++)
	if(atype<RR*>()->CastingFrom(tt[i].left() ) ) 
	  {
          tab[i]=atype<RR*>()->CastTo(tt[i]);
	    what[i]=0;
	  }
	else if(atype<KN_<RR> >()->CastingFrom(tt[i].right() ) ) 
	  {
	    tab[i]=atype<KN_<RR> >()->CastTo(tt[i].RightExp());
	    what[i]=1;
	  }      
	else 
	  CompileError(" we are waiting for scalar or vector of scalar");
    }
    
    AnyType operator()(Stack stack)  const 
    {
      A  a=GetAny<A>((*a0)(stack));
      KN<AnyType> v(N);
      KN<int>  nn(N+1);
      for (int i=0;i<N;i++)
        v[i]= (*(tab[i]))(stack);
      
      int n=0; 
      for (int i=0;i<N;i++)
	{
	  if (what[i]==0) nn[i]=1;
	  else if (what[i]==1) nn[i]=GetAny<KN_<RR> >(v[i]).size();
          n += nn[i];
	}
      ffassert(n == a.size()); 
      for (int i=0,j=0 ;i<N; j += nn[i++])
	
        if (what[i]==0)
         * GetAny<RR*>(v[i]) = a[j];
        else if (what[i]==1) { // hack FH 
           KN_<RR> tab(GetAny<KN_<RR> >(v[i])); 
           tab  =a(SubArray(nn[i],j));
           }
      return SetAny<R>(a);
    } 
    bool MeshIndependent() const     {return  mi;} // 
    ~CODE() { delete [] tab; delete[] what;}
    operator aType () const { return atype<R>();}    
  }; // end sub class CODE
  
  
    public: // warning hack  A and B 
    E_F0 * code(const basicAC_F0 & args) const 
     { return  new CODE(t[1]->CastTo(args[1]),*dynamic_cast<const E_Array*>( t[0]->CastTo(args[0]).RightValue()));} 
    SetArrayofKNfromKN():   OneOperator(atype<R>(),atype<B>(),atype<A>())  {} // warning with A and B 
  
};
   


template<class K> long get_n(KN<K> * p){ return p->N();}
template<class K> long get_n(KNM<K> * p){ return p->N();}

template<class K> long get_m(KNM<K> * p){ return p->M();}
template<class K> K get_max(KN<K> * p){ return p->max();}
template<class K> K get_min(KN<K> * p){ return p->min();}
template<class K> K get_sum(KN<K> * p){ return p->sum();}
template<class K> double get_l2(KN<K> * p){ return p->l2();}
template<class K> double get_l1(KN<K> * p){ return p->l1();}
template<class K> double get_linfty(KN<K> * p){ return p->linfty();}

template<class K,class T > K get_sum0(const T & p){ return p.sum();}
template<class K,class T > K get_max0(const T &p){ return p.max();}
template<class K,class T > K get_min0(const T &p){ return p.min();}
template<class K,class T> K  get_l2_0(const T &p){ return p.l2();}
template<class K,class T> K  get_l1_0(const T &p){ return p.l1();}
template<class K,class T> K  get_linfty_0(const T &p){ return p.linfty();}

 
 class ostream_precis { public:
 ostream_precis(ostream * ff) :f(ff) {}
  ostream * f;
   operator long () const {return f->precision();}
 };

template<class A,class B> B castto(const A & a){ return a;}
 
template<class K>
void ArrayDCL()
{
    Dcl_TypeandPtr<KN<K> >(0,0,0,::Destroy<KN<K> >);
  //  Dcl_Type<KN<Complex> *>(0,::Destroy<KN<Complex> >);
   // Dcl_Type<KN<K> *>(0,::Destroy<KN<K> >); // Modif 17102005 
   // attention un exp KN<> * right est un KN<> et non un KN<> *

    Dcl_Type<KNM<K> *>(0,::Destroy<KNM<K> >);
   //  Dcl_Type< Transpose<KN<K> *> > ();  remove mars 2006 FH 
    Dcl_Type< outProduct_KN_<K>* >();
    Dcl_Type< Transpose<KN_<K> > > ();
    Dcl_Type< Transpose< KNM<K> *> >();
    //Dcl_Type< Transpose<KN<Complex> > > ();
    Dcl_TypeandPtr<KN_<K> >(0,0,0,0);
    //Dcl_TypeandPtr<KN_<Complex> >(0,0,0,0);

    Dcl_Type<Add_KN_<K> >();
    
    Dcl_Type<DotStar_KN_<K> >();
    Dcl_Type<DotSlash_KN_<K> >();
    Dcl_Type<Sub_KN_<K> >();
    Dcl_Type<Mulc_KN_<K> >();
    Dcl_Type<Mul_KNM_KN_<K> >();
    Dcl_Type<Add_Mulc_KN_<K> *>();
    Dcl_Type<if_arth_KN_<K> *>();
    // for    B(I) and B(I^-1)
    Dcl_Type<pair<KN_<K>,Inv_KN_long> *>();
    Dcl_Type<pair<KN_<K>,KN_<long> > *>();
    

     map_type[typeid(KN_<K> ).name()]->AddCast(
       new E_F1_funcT<KN_<K>,KN_<K>*>(UnRef<KN_<K> >),
     //  new E_F1_funcT<KN_<K>,KN<K>*>(UnRef<KN_<K>,KN<K>* >), inutil cas KN<K> est right expression de KN<K>* 
       new E_F1_funcT<KN_<K>,KN<K> >(Cast<KN_<K>,KN<K> >)
       
       );
    //   ,new E_F1_funcT<KN_<K>,K>(ValueToKN_<K>),
    //   new E_F1_funcT<KN_<K>,K*>(PtrToKN_<K>)       
       
     // Ajoute FH   
     map_type[typeid(KN<K> ).name()]->AddCast(
       new E_F1_funcT<KN<K>,KN<K>*>(UnRef<KN<K> >)
    //   ,new E_F1_funcT<KN_<K>,K>(ValueToKN_<K>),
    //   new E_F1_funcT<KN_<K>,K*>(PtrToKN_<K>)       
       ); 
    map_type_of_map[make_pair(atype<long>(),atype<K>())]=atype<KN<K>*>(); // vector
    map_pair_of_type[make_pair(atype<long>(),atype<long>())] =atype<pair<long,long> >();   
    map_type_of_map[make_pair(atype<pair<long,long> >(),atype<K>())]=atype<KNM<K>*>(); // matrix                                               
}



template<class A,class B> pair<A,B> * pBuild(const A & a,const B & b)
  { return new pair<A,B>(a,b);}

// add mars 2006
template<class K,class L,class OP>
struct set_A_BI: public binary_function<KN_<K>,pair<KN_<K>, KN_<L> > *,KN_<K> > {
  static KN_<K> f(const KN_<K>   & a, pair<KN_<K>, KN_<L> > * const & b)  {
    KN_<K> x(a);
    OP op;
     const KN_<K> & y(b->first);
    const KN_<L> & I(b->second);
    L  N = x.N();
    L n = y.N();
    
    L maxI=I(SubArray(N)).max() ;
    L minI=I(SubArray(N)).min() ;
    
    if( maxI >= n || I.N()  < N) 
       { cerr << " Out of Bound x=y(I)  :  0 <= " << minI << " < "<< maxI << "< " << n  << endl;
         cerr << " or I.N() " << I.N() << " > " << N << endl;
         ExecError("Out of Bound error");
       }
       
    for(int i=0;i<N;i++)
      if(I[i]>=0) 
      op(x(i),y(I[i]));
    delete b;
    return a;
  
  }
};  

template<class K,class L,class OP>
struct set_AI_B: public binary_function<pair<KN_<K>, KN_<L> > * ,KN_<K>, NothingType > {
  static NothingType  f( pair<KN_<K>, KN_<L> > * const & b,const KN_<K>   & a)  {
    KN_<K> x(a);
    OP op;
     const KN_<K> & y(b->first);
    const KN_<L> & I(b->second);
    L  N = I.N();
    L n = y.N();
    
    L maxI=I(SubArray(N)).max() ;
    L minI=I(SubArray(N)).min() ;
    
    if(  maxI >= n || x.N()  < N ) 
       { cerr << " Out of Bound x(I)=y  :  0 <= " << minI << " < "<< maxI << "< " << n  << endl;
         cerr << " or x.N() " << I.N() << " > " << N << endl;
         ExecError("Out of Bound error");
       }
       
    for(int i=0;i<N;i++)
      if(I[i] >=0) 
      op(y(I[i]),x[i]);
    delete b;   
    return  NothingType();
  
  }
};  

template<class K> 
struct Op3_paac: public ternary_function<KN_<K>,KN_<K>,K,if_arth_KN_<K>*> { 
static if_arth_KN_<K>* f(Stack s,const KN_<K> & a,const KN_<K> & b,const  K & c )  {
    //K cc(c);
    KN_<K> kc(new(NewAllocTmp(s,sizeof(c))) K(c),1,0);
  return new if_arth_KN_<K>(a,b,kc);}
};   
template<class K> 
struct Op3_paca: public ternary_function<KN_<K>,K,KN_<K>,if_arth_KN_<K>*> { 
    static if_arth_KN_<K>* f(Stack s,const KN_<K> & a,const  K & b,const KN_<K> & c )  {
	//K bb(b);
	KN_<K> kb(new(NewAllocTmp(s,sizeof(b))) K(b),1,0);
    return new if_arth_KN_<K>(a,kb,c);}
};   

template<class K> 
struct Op3_pacc: public ternary_function<KN_<K>,K,K,if_arth_KN_<K>*> { 
    static if_arth_KN_<K>* f(Stack s,const KN_<K> & a,const K & b,const  K & c )  {
	K cc(c),bb(b);
	KN_<K> kc(new(NewAllocTmp(s,sizeof(c))) K(c),1,0),
	       kb(new(NewAllocTmp(s,sizeof(b))) K(b),1,0);    
    return new if_arth_KN_<K>(a,kb,kc);}
};   

  
extern aType aaaa_knlp;
template<class K,class Z>
void ArrayOperator()
{
     Dcl_Type< Resize<KN<K> > > ();
     Dcl_Type< Resize<KNM<K> > > ();
     aType knrp = atype<KN<K> *>();
     aType knr_ = atype<KN_<K> >();
     typedef KN<Z> ZN;
      
     aType knlp=  aaaa_knlp ;
     //atype<KN<Z> *>();
     // aType knl_ = atype<KN_<Z> >();
    
     atype<KN<K>* >()->Add("[","",new OneOperator2_<K*,KN<K>*,Z >(get_elementp_<K,KN<K>*,Z>));
     atype<KN<K>* >()->Add("(","",new OneOperator2_<K*,KN<K>*,Z >(get_elementp_<K,KN<K>*,Z>));
     atype<KN_<K> >()->Add("(","",new OneOperator2_<KN_<K>,KN_<K>,char >(fSubArrayc<KN_<K> >));
     atype<KN_<K> >()->Add("(","",new OneOperator2_<KN_<K>,KN_<K>,SubArray>(fSubArray<K> ));
     atype<KN<K>*>()->Add("(","",new OneOperator2_<KN_<K>,KN<K>*,SubArray>(fSubArrayp<K> ));
     atype<KN<K>* >()->Add("(","",new OneOperator2_<KN<K>*,KN<K>*,char >(fSubArrayc<KN<K>* >));

     atype<KNM<K>* >()->Add("(","",new OneOperator3_<KN_<K>,KNM<K>*,Z,SubArray >(get_element_is<KN_<K>,KNM<K>*,Z,SubArray>));
     atype<KNM<K>* >()->Add("(","",new OneOperator3_<KN_<K>,KNM<K>*,SubArray,Z >(get_element_si<KN_<K>,KNM<K>*,SubArray,Z>));
     atype<KNM<K>* >()->Add("(","",new OneOperator3_<KN_<K>,KNM<K>*,Z,char >(get_element_lineorcol<KN_<K>,KNM<K>*,Z,char>));
     atype<KNM<K>* >()->Add("(","",new OneOperator3_<KN_<K>,KNM<K>*,char,Z >(get_element_lineorcol<KN_<K>,KNM<K>*,char,Z>));

     atype<KNM<K>* >()->Add("(","",new OneOperator3_<K*,KNM<K>*,Z,Z >(get_elementp2_<K,KNM<K>*,Z,Z>));

     Add<KN<K> *>("sum",".",new OneOperator1<K,KN<K> *>(get_sum));
     Add<KN<K> *>("min",".",new OneOperator1<K,KN<K> *>(get_min));
     Add<KN<K> *>("max",".",new OneOperator1<K,KN<K> *>(get_max));
     Add<KN<K> *>("l2",".",new OneOperator1<double,KN<K> *>(get_l2));
     Add<KN<K> *>("l1",".",new OneOperator1<double,KN<K> *>(get_l1));
     Add<KN<K> *>("linfty",".",new OneOperator1<double,KN<K> *>(get_linfty));
     
     
     Add<KN_<K> >("sum",".",new OneOperator1_<K,KN_<K> >(get_sum0<K,KN_<K> >));
     Add<KN_<K> >("min",".",new OneOperator1_<K,KN_<K> >(get_min0<K,KN_<K> >));
     Add<KN_<K> >("max",".",new OneOperator1_<K,KN_<K> >(get_max0<K,KN_<K> >));
     Add<KN_<K> >("l2",".",new OneOperator1_<double,KN_<K> >(get_l2_0<double,KN_<K> >));
     Add<KN_<K> >("l1",".",new OneOperator1_<double,KN_<K> >(get_l1_0<double,KN_<K> >));
     Add<KN_<K> >("linfty",".",new OneOperator1_<double,KN_<K> >(get_linfty_0<double,KN_<K> >));
    
     Add<KN<K> >("sum",".",   new OneOperator1_<K,KN<K> >(get_sum0<K,KN<K> >));
     Add<KN<K> >("min",".",   new OneOperator1_<K,KN<K> >(get_min0<K,KN<K> >));
     Add<KN<K> >("max",".",   new OneOperator1_<K,KN<K> >(get_max0<K,KN<K> >));
     Add<KN<K> >("l2",".",    new OneOperator1_<double,KN<K> >(get_l2_0<double,KN<K> >));
     Add<KN<K> >("l1",".",    new OneOperator1_<double,KN<K> >(get_l1_0<double,KN<K> >));
     Add<KN<K> >("linfty",".",new OneOperator1_<double,KN<K> >(get_linfty_0<double,KN<K> >));
     

     Add<KN<K> *>("resize",".",new OneOperator1< Resize<KN<K> >,KN<K> *>(to_Resize));
     Add<KNM<K> *>("resize",".",new OneOperator1< Resize<KNM<K> >,KNM<K> *>(to_Resize));
     
     Add<Resize<KN<K> > >("(","",new OneOperator2_<KN<K> *,Resize<KN<K> > , Z   >(resize1));
     Add<Resize<KNM<K> > >("(","",new OneOperator3_<KNM<K> *,Resize<KNM<K> > , Z, Z  >(resize2));