afunction.hpp

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

template<class R,class A>
 class E_F1_F :public  E_F1 { public:
  typedef  R (*func)(A) ; 
  func f;
  E_F1_F(func ff) : f(ff) {}
  AnyType operator()(Stack s,AnyType & a)  const 
    {return SetAny<R>(f(GetAny<A>(a)));}  
};

template<class R,class A0,class A1>
 class E_F2_F :public  E_F2 { public:
  typedef  R (*func)(const  A0 &,const  A1&) ; 
  func f;
  E_F2_F(func ff) : f(ff) {}
  AnyType operator()(Stack s,AnyType & a0,AnyType & a1)  const 
    {return SetAny<R>(f(GetAny<A0>(a0),GetAny<A1>(a1)));}  
};

template<class R,class TA0,bool RO=true>
 class E_F_F0 :public  E_F0 { public:
   template <class T> struct remove_reference     {typedef T type;};
//   template <class T> struct remove_reference<T&> {typedef T type;};
   template <class T> struct remove_reference<const T&> {typedef T type;};
   typedef typename remove_reference<TA0>::type A0;
   
 
  typedef  R (*func)(  TA0 ) ; 
  func f;
  Expression a;
  E_F_F0(func ff,Expression aa) : f(ff),a(aa) {}
  AnyType operator()(Stack s)  const 
    {return SetAny<R>(f(GetAny<A0>( (*a)(s) )));}  
  bool EvaluableWithOutStack() const 
      {return a->EvaluableWithOutStack();} // 
  bool MeshIndependent() const {return a->MeshIndependent();} // 
  bool ReadOnly() const { return RO  ;}     
  int compare (const E_F0 *t) const { 
     int rr;
    // cout << "cmp " << typeid(*this).name() << " and " << typeid(t).name() << endl;
     const  E_F_F0* tt=dynamic_cast<const E_F_F0 *>(t);
     if (tt && f == tt->f) rr = a->compare(tt->a);
     else rr = E_F0::compare(t);
     return rr;
     } // to give a order in instuction 

   int Optimize(deque<pair<Expression,int> > &l,MapOfE_F0 & m, size_t & n) ;
   virtual ostream & dump(ostream &ff) const  { ff << typeid(*this).name() <<" f= " << ff << " a= "<< *a << ' '  ;return ff; }

};

// modif for xlc++ FH 
template<class R,class TA0,bool RO=true>
class E_F_F0_Opt: public E_F_F0<R,TA0,RO>  { public :
  size_t ia;  
  E_F_F0_Opt(const  E_F_F0<R,TA0,RO>  &t,size_t iaa) 
    : E_F_F0<R,TA0,RO>(t) , ia(iaa) {assert(iaa<2000000 && iaa >0);}
  AnyType operator()(Stack s)  const 
  {
    // A0 x =  *static_cast<A0 *>(static_cast<void*>(static_cast<char *>(s)+ia));
    // cout << " opt f (" << x << " ) = "   << ": " << ia << endl; 
    return SetAny<R>( f( *static_cast<typename E_F_F0<R,TA0>::A0 *>(static_cast<void*>(static_cast<char *>(s)+ia))  ) );}  
  
};   

template<class R,class TA0,bool RO>   
   int  E_F_F0<R,TA0,RO>::Optimize(deque<pair<Expression,int> > &l,MapOfE_F0 & m, size_t & n) 
    {
       int rr = find(m);
       if (rr) return rr;
       return insert(new E_F_F0_Opt<R,TA0,RO>(*this,a->Optimize(l,m,n)),l,m,n);       
    } 
// fin modif xlc++ 

template<class A0>
 class E_VF_F0 :public  E_F0 { public:
  typedef  void (*func)(  A0 ) ; 
  func f;
  Expression a;
  E_VF_F0(func ff,Expression aa) : f(ff),a(aa) {}
  AnyType operator()(Stack s)  const 
    {f(GetAny<A0>( (*a)(s) ));return Nothing;}  
  bool EvaluableWithOutStack() const 
      {return a->EvaluableWithOutStack();} // 
      
 bool MeshIndependent() const { return a->MeshIndependent();  }    

};

inline int clexico(int i,int j) { return i==0 ? j : i;}
inline int clexico(int i,int j,int k) { int ll=clexico(i,j);  return  ll==0 ? k : ll;}

template<class R,class TA0,class TA1>
 class E_F_F0F0 :public  E_F0 { public:
   template <class T> struct remove_reference     {typedef T type;};
   template <class T> struct remove_reference<T&> {typedef T type;};
   typedef typename remove_reference<TA0>::type A0;
   typedef typename remove_reference<TA1>::type A1;
   typedef  R (*func)( A0 , A1 ) ;
    
  func f;
  Expression a0,a1;
  E_F_F0F0(func ff,Expression aa0,Expression aa1) 
    : f(ff),a0(aa0),a1(aa1) {}
  AnyType operator()(Stack s)  const 
    {return SetAny<R>( f( GetAny<A0>((*a0)(s)) , GetAny<A1>((*a1)(s)) ) );}  
   bool EvaluableWithOutStack() const 
      {return a0->EvaluableWithOutStack() && a1->EvaluableWithOutStack();} // 
   bool MeshIndependent() const 
      {return a0->MeshIndependent() && a1->MeshIndependent();} // 
  int compare (const E_F0 *t) const { 
     int rr;
    // cout << "cmp " << typeid(*this).name() << " and " << typeid(t).name() << endl;
     const  E_F_F0F0* tt=dynamic_cast<const E_F_F0F0 *>(t);
     if (tt && f == tt->f) rr= clexico(a0->compare(tt->a0),a1->compare(tt->a1));
     else rr = E_F0::compare(t);
     return rr;
     } // to give a order in instuction 
      
   int Optimize(deque<pair<Expression,int> > &l,MapOfE_F0 & m, size_t & n) ;
    
};

// modif for xlc++
template<class R,class TA0,class TA1>
class E_F_F0F0_Opt: public E_F_F0F0<R,TA0,TA1>  { public :
  size_t ia,ib;  
  E_F_F0F0_Opt(const  E_F_F0F0<R,TA0,TA1> &t,size_t iaa,size_t ibb) 
    : E_F_F0F0<R,TA0,TA1>(t) ,
      ia(iaa),ib(ibb) {}
  AnyType operator()(Stack s)  const 
  {
    //A0 aa =*static_cast<A0 *>(static_cast<void*>(static_cast<char *>(s)+ia));
    //A1 bb=*static_cast<A1 *>(static_cast<void*>(static_cast<char *>(s)+ib)) ;
    //cout << ia << " " << ib <<  "f( " << aa << "," << bb  << " )   = "<< f(aa,bb) << endl;
    return SetAny<R>( f( *static_cast<typename E_F_F0F0<R,TA0,TA1>::A0 *>(static_cast<void*>(static_cast<char *>(s)+ia)) , 
			 *static_cast<typename E_F_F0F0<R,TA0,TA1>::A1 *>(static_cast<void*>(static_cast<char *>(s)+ib)) ) );}  
  
};     
       

template<class R,class TA0,class TA1>
   int E_F_F0F0<R,TA0,TA1>::Optimize(deque<pair<Expression,int> > &l,MapOfE_F0 & m, size_t & n) 
    {

       int rr = find(m);
       if (rr) return rr;

       return insert(new E_F_F0F0_Opt<R,TA0,TA1>(*this,a0->Optimize(l,m,n),a1->Optimize(l,m,n)),l,m,n);
    }
// add modif for xlc++



template<class R,class A0>
 class E_F_F0_ :public  E_F0 { public:
  typedef  R (*func)(const   A0& ) ; 
  func f;
  Expression a;
  E_F_F0_(func ff,Expression aa) : f(ff),a(aa) {}
  AnyType operator()(Stack s)  const 
    {return SetAny<R>(f(GetAny<A0>( (*a)(s) )));}  
   bool EvaluableWithOutStack() const 
      {return a->EvaluableWithOutStack() ;} // 
   bool MeshIndependent() const 
      {return a->MeshIndependent();} // 
    
};
// add FH 07/2008  for pmesh clean 
template<class R,class A0>
class E_F_F0_Add2RC :public  E_F0 { public:
    typedef  R (*func)(const   A0& ) ; 
    func f;
    Expression a;
    E_F_F0_Add2RC(func ff,Expression aa) : f(ff),a(aa) {}
    AnyType operator()(Stack s)  const 
    {return SetAny<R>(Add2StackOfPtr2FreeRC(s,f(GetAny<A0>( (*a)(s) ))));}  
    bool EvaluableWithOutStack() const 
    {return a->EvaluableWithOutStack() ;} // 
    bool MeshIndependent() const 
    {return a->MeshIndependent();} // 
    
};
// end add. 
template<class R,class A0>
 class E_F_F0s_ :public  E_F0mps { public:
  typedef  R (*func)(Stack stack,const   A0& ) ; 
  func f;
  Expression a;
  E_F_F0s_(func ff,Expression aa) : f(ff),a(aa) {}
  AnyType operator()(Stack s)  const 
    {return SetAny<R>(f(s,GetAny<A0>( (*a)(s) )));}  
  bool MeshIndependent() const 
      {return true;} // 

    operator aType () const { return atype<R>();}         
    
};

template<class R,class A0,class A1,class E=E_F0>
 class E_F_F0F0_ :public  E { public:
  typedef  R (*func)(const  A0 &,const  A1 & ) ; 
  func f;
  Expression a0,a1;
  E_F_F0F0_(func ff,Expression aa0,Expression aa1) 
    : f(ff),a0(aa0),a1(aa1) {}
  AnyType operator()(Stack s)  const 
    {return SetAny<R>( f( GetAny<A0>((*a0)(s)) , GetAny<A1>((*a1)(s)) ) );} 
    bool MeshIndependent() const 
      {return a0->MeshIndependent() && a1->MeshIndependent();} // 
 
};
// FH Add 07/2008 
//   class with add 1 to the refcounter  for  mesh . 
template<class R,class A0,class A1,class E=E_F0>
class E_F_F0F0_Add2RC :public  E { public:
typedef  R (*func)(const  A0 &,const  A1 & ) ; 
func f;
Expression a0,a1;
E_F_F0F0_Add2RC(func ff,Expression aa0,Expression aa1) 
: f(ff),a0(aa0),a1(aa1) {}
AnyType operator()(Stack s)  const 
{return SetAny<R>(Add2StackOfPtr2FreeRC(s, f( GetAny<A0>((*a0)(s)) , GetAny<A1>((*a1)(s)) ) ));} 
bool MeshIndependent() const 
{return a0->MeshIndependent() && a1->MeshIndependent();} // 

};
// FH end 07/2008 
template<class R,class A0,class A1,class A2,class E=E_F0>
 class E_F_F0F0F0_ :public  E { public:
  typedef  R (*func)(const  A0 &,const  A1 & , const A2 &) ; 
  func f;
  Expression a0,a1,a2;
  E_F_F0F0F0_(func ff,Expression aa0,Expression aa1,Expression aa2) 
    : f(ff),a0(aa0),a1(aa1),a2(aa2) {}
  AnyType operator()(Stack s)  const 
    {return SetAny<R>( f( GetAny<A0>((*a0)(s)) , GetAny<A1>((*a1)(s)),GetAny<A2>((*a2)(s))  ) );}  
    virtual size_t nbitem() const {return a2->nbitem(); } 
      bool MeshIndependent() const 
      {return a0->MeshIndependent() && a1->MeshIndependent()&& a2->MeshIndependent();} // 

};

template<class R,class A0,class A1,class A2,class E=E_F0>
 class E_F_stackF0F0F0_ :public  E_F0mps { public:
  typedef  R (*func)(Stack, const  A0 &,const  A1 & , const A2 &) ; 
  func f;
  Expression a0,a1,a2;
  E_F_stackF0F0F0_(func ff,Expression aa0,Expression aa1,Expression aa2) 
    : f(ff),a0(aa0),a1(aa1),a2(aa2) {}
  AnyType operator()(Stack s)  const 
    {return SetAny<R>( f(s, GetAny<A0>((*a0)(s)) , GetAny<A1>((*a1)(s)),GetAny<A2>((*a2)(s))  ) );}  
    virtual size_t nbitem() const {return a2->nbitem(); } 
   bool MeshIndependent() const { return true;}
};

template<class R,class A0,class A1>
 class E_F_F0F0_NC :public  E_F0 { public:
  typedef  R (*func)(  A0 &,const  A1 & ) ; 
  func f;
  Expression a0,a1;
  E_F_F0F0_NC(func ff,Expression aa0,Expression aa1) 
    : f(ff),a0(aa0),a1(aa1) {}
  AnyType operator()(Stack s)  const 
    {return SetAny<R>( f( GetAny<A0>((*a0)(s)) , GetAny<A1>((*a1)(s)) ) );}  
   bool MeshIndependent() const 
      {return a0->MeshIndependent() && a1->MeshIndependent() ; } // 

};




 class E_F_StackF0F0 :public  E_F0mps { public:
  typedef   AnyType (*func)(Stack,Expression ,Expression ) ; 
  func f;
  Expression a0,a1;
  E_F_StackF0F0(func ff,Expression aa0,Expression aa1) 
    : f(ff),a0(aa0),a1(aa1) { }
  AnyType operator()(Stack s)  const 
    {return  (*f)(s, a0 , a1) ;}  

};


/*
 class E_F_F0F0_<AnyType,AnyType,AnyType> :public  E_F0 { public:
  typedef  AnyType (*func)(const  AnyType &,const  AnyType & ) ; 
  func f;
  Expression a0,a1;
  E_F_F0F0_(func ff,Expression aa0,Expression aa1) 
    : f(ff),a0(aa0),a1(aa1) {}
  AnyType operator()(Stack s)  const 
    {return  f( (*a0)(s) , (*a1)(s) );}  
   bool MeshIndependent() const 
      {return a0->MeshIndependent() && a1->MeshIndependent() ; } // 

};
*/

class E_F2_func :public  E_F2 { public:
   Function2 f;
   AnyType operator()(Stack s,AnyType & a,AnyType & b)  const {return f(s,a,b);}
   E_F2_func( Function2 ff) : f(ff) {}
};

class E_F0_Func1 :public  E_F0 { public:
   Function1  f;
   E_F0 *a;
   AnyType operator()(Stack s)  const {return f(s,(*a)(s));}
   E_F0_Func1( Function1 f1,E_F0 *aa) : f(f1),a(aa){}
   bool EvaluableWithOutStack() const {return a->EvaluableWithOutStack();} // 
   bool MeshIndependent() const {return a->MeshIndependent();} // 
   int compare (const E_F0 *t) const { 
     int rr;
     const  E_F0_Func1* tt=dynamic_cast<const E_F0_Func1 *>(t);
     if (tt && f == tt->f) rr = a->compare(tt->a);
     else rr = E_F0::compare(t);
     if(tt && 0)
      {
       cout << "\n\t\t\t -------- " << (void *) f << " " << (void *) tt->f << " rr=" << a->compare(tt->a) << endl;
       cout << "\t\t\tcmp E_F0_Func1 " << rr <<" << " << *this << " cmp " << *t << " " << tt << ">>\n";
      }
     return rr;
     } // to give a order in instuction 
  // int Optimize(deque<pair<Expression,int> > &l,MapOfE_F0 & m, size_t & n) const;  // build optimisation

    virtual ostream & dump(ostream &ff) const  { ff << "E_F0_Func1 f= " << f << " a= "<< *a << ' '  ;return ff; }

};
class E_F0_Func2 :public  E_F0 { public:
   Function2  f;
   E_F0 *a,*b;
   AnyType operator()(Stack s)  const {return f(s,(*a)(s),(*b)(s));}
   E_F0_Func2( Function2 f1,E_F0 *aa,E_F0 *bb) : f(f1),a(aa),b(bb){}
   bool EvaluableWithOutStack() const 
     {return a->EvaluableWithOutStack() &&b->EvaluableWithOutStack();} // 
   bool MeshIndependent() const {return a->MeshIndependent() && b->MeshIndependent();} // 

};



//  the variable offset / stack (local variable)
template<class R> class Value1:public E_F0
 { 
  size_t offset;
  public:
  AnyType operator()(Stack s) const { return SetAny<R*>(static_cast<R *>(static_cast<void *>(  static_cast<char *>(s)+offset)));}
  Value1(size_t o):offset(o) {}
};

//  the variable globale
template<class R> class GValue:public E_F0
 { 
  mutable R v;
  public:
  AnyType operator()(Stack ) const { return SetAny<R*>(static_cast<R *>(static_cast<void *>(&v)));}
  GValue(R o):v(o) {}
  bool EvaluableWithOutStack() const {return true;} // 
  
};

//  a constante value 
template<class R> int ccompare(const R & a,const R& b){ return a==b ?  0 :( a<b ? -1 : +1);}
template<class R> int ccompare(const complex<R> & a,const complex<R>& b){ 
  int c=ccompare(a.real(),b.real()); 
  return c==0 ? ccompare(a.imag(),b.imag()): c ;}
  
template<class R> class EConstant:public E_F0
 { 
  const R v;
  public:
  AnyType operator()(Stack ) const { /*cout << " ()" << v << endl*/;return SetAny<R>(v);}
  EConstant(const R & o):v(o) { /*cout << "New constant " << o << endl;*/}
  bool EvaluableWithOutStack() const {return true;} //   
  operator aType () const { return atype<R>();} 
  int compare (const E_F0 *t) const { 
        int rr;
        const  EConstant * tt=dynamic_cast<const EConstant *>(t);
            if (tt) rr = ccompare(v,tt->v);
             else rr = E_F0::compare(t);
         return rr;
       } 
   ostream & dump(ostream &f) const { f << " ((" <<typeid(R).name()  << ") " << v << ") " ;return f;}  
};





//  the variable offset / stack (local variable)

 class LocalVariable:public E_F0
 { 
  size_t offset;
  aType t; //  type of the variable just for check  
  public:
  AnyType operator()(Stack s) const { 
    SHOWVERB( cout << "\n\tget var " << offset << " " <<  t->name() << endl);  
//   return PtrtoAny(static_cast<void *>(static_cast<char *>(s)+offset),t);}
   return PtrtoAny(Stack_offset<void>(s,offset),t);}

  LocalVariable(size_t o,aType tt):offset(o),t(tt) {throwassert(tt);     
     SHOWVERB(cout << "\n--------new var " << offset << " " <<  t->name() << endl);
    }
};


class LocalVariableFES : public LocalVariable { public:
  size_t data;
  LocalVariableFES(size_t o,aType tt,const  size_t & d) 
   : LocalVariable(o,tt),data(d) {}
  size_t nbitem() const { return data;}
};