afunction.hpp

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

template <class U>
class LocalVariablePlus : public LocalVariable { public:
  U data;
  LocalVariablePlus(size_t o,aType tt,const  U & d) 
   : LocalVariable(o,tt),data(d) {}
};

//  global variable bof bof 
template<class T> class PValue:public E_F0
 { 
  T * p;
  public:
  AnyType operator()(Stack  ) const { return p;}
  PValue(T * pp):p(pp) {}
};
template<class R> class PPValue:public E_F0
 { 
  R ** p;
  public:
  AnyType operator()(Stack ) const { return SetAny<R*>(*p);}
  PPValue(R ** pp):p(pp) {}
};


template<class R>
Type_Expr CPValue(R & v)
 {
   throwassert(map_type[typeid(R*).name()]);
  return make_pair(map_type[typeid(R*).name()],new PValue<R>(&v));
 }
template<class R>
Type_Expr CPPValue(R *& v)
 {
   throwassert(map_type[typeid(R*).name()]);
  return make_pair(map_type[typeid(R*).name()],new PPValue<R>(&v));
 }
 
template<class R >
Type_Expr CConstant(const R & v)
 {
  throwassert(map_type[typeid(R).name()]);
  return make_pair(map_type[typeid(R).name()],new  EConstant<R>(v));
 }



class CC_F0 {
  Expression f;
  aType r;
public:
  void operator=(const C_F0& c) { f=c.f;r=c.r;;} 
  void operator=(const AC_F0& a) ; //{ f=new E_Array(a); f= atype<E_Array>();};
  void operator=(long ) {f=0;r=0;}
  void operator=(const CListOfInst& c);//{ C_FO cc=c;f=cc.f;r=cc.r}
  operator C_F0 () const {return C_F0(f,r);}
  bool Empty() const {return !f || f->Empty();}
  aType left() const {return r;}
 // operator const C_F0 &() const {return  *this;}
};


class ListOfInst : public  E_F0mps { 
    int n;
    Expression   *   list;
    int   *   linenumber;
    const int nx;
    public:
    ListOfInst():n(0),list(0),linenumber(0),nx(10){}
    ListOfInst(int nn):n(0),list(0),linenumber(0),nx(nn?nn:10){}
    void Add(const C_F0 & ins); 
    AnyType operator()(Stack s) const; 
    operator aType () const { return n ? (aType) * (list[n-1]) : atype<void>();} 
   
   Expression &operator[](int i){return list[i];}
   bool empty() const {return n==0;}
   int size() const {return n;}
   Expression * ptr() const {return list;}
   int * nlines() const {return linenumber;}
 //  void destroy() { if (list) delete [] list; list=0;}
  ~ListOfInst(){ cout << " ----- ~ListOfInst " << endl;
  if(list) delete [] list;list=0;if(linenumber)  delete[] linenumber; linenumber=0;}
};

class CListOfInst  {  private:
  ListOfInst * f;
  const basicForEachType *r;
  public:
   void operator=(const CC_F0 &a){
     f=new ListOfInst();     
       if( !a.Empty() ) {
         f->Add(a);
         r=a.left(); }}
   CListOfInst & operator+=(const CC_F0 & a);//{ if( !a.Empty()){ f->Add(a);r=a.left();};return *this;} 
    operator C_F0 () const  { return C_F0(f,r);}
   void eval(Stack s) {(*f)(s);}
   int size() const {return f->size();}
   Expression * ptr() const {return f->ptr();}
   int * nlines() const { return f->nlines();}
};


AnyType FWhile(Stack ,E_F0 * test,E_F0 * ins);
AnyType FFor(Stack s ,E_F0 * i0,E_F0 * i1,E_F0 * i2,E_F0 * ins);
AnyType FIf(Stack s ,E_F0 * test,E_F0 * i1,E_F0 * i2,E_F0 * notuse);
AnyType TTry(Stack s ,E_F0 * i0,E_F0 * i1,E_F0 * i2,E_F0 * notuse);



extern TableOfIdentifier Global;
void ShowType(ostream & );

template<class T> 
inline C_F0 to(const C_F0 & a) { return map_type[typeid(T).name()]->CastTo(a);}


/*
inline C_F0 toBool(const C_F0 & a)    {return ATYPE(bool)->CastTo(a);}
inline C_F0 toInt(const C_F0 & a)     {return ATYPE(int)->CastTo(a);}
inline C_F0 toLong(const C_F0 & a)    {return ATYPE(long)->CastTo(a);}
inline C_F0 toDouble(const C_F0 & a)  {return ATYPE(double)->CastTo(a);}
inline C_F0 toComplex(const C_F0 & a)  {return ATYPE(Complex)->CastTo(a);}
*/
inline C_F0 While(C_F0 test,C_F0 ins) {return C_F0(new E_F0_CFunc2(FWhile,to<bool>(test),ins),0);}
inline C_F0 For(C_F0 i0,C_F0 i1,C_F0 i2,C_F0 ins) {return C_F0(new E_F0_CFunc4(FFor,i0,to<bool>(i1),i2,ins),0);}
inline C_F0 Try(C_F0 i0,C_F0 i1,C_F0 i2)  {return C_F0(new E_F0_CFunc4(TTry,i0,i1,i2,0),0);}
inline C_F0 FIf(C_F0 i0,C_F0 i1,C_F0 i2) {return C_F0(new E_F0_CFunc4(FIf,to<bool>(i0),i1,i2,0),0);}
inline C_F0 FIf(C_F0 i0,C_F0 i1) {return C_F0(new E_F0_CFunc4(FIf,to<bool>(i0),i1,0,0),0);}
//inline  C_F0 C_F0::PtrValue() const{ 
//   if (!(r && r->un_ptr)) { cerr << "PtrValue: Not a Left value " << *r << endl;CompileError();} 
//   return C_F0(new  E_F0_Func1(r->un_ptr->f,f),r->un_ptr->r);}

 

class basicAC_F0 {
//  version de base d'un tableau d'un parametres  
//   pour les operateurs unaire, binaire, , 
//   pas d'allocation 
  friend class E_Array; // for mapping fonction 
  protected:
  typedef const C_F0 const_C_F0;
 int nb;
 C_F0 *a;
 public:
  typedef map<const char *,C_F0,Keyless> maptype ;
  typedef maptype::iterator iterator;
  typedef maptype::const_iterator const_iterator;
  maptype * named_parameter;
 basicAC_F0 & operator=(int i) {throwassert(i==0);named_parameter=0,nb=0;return *this;} // pas de parametres
 basicAC_F0 & operator=(C_F0 & c) {named_parameter=0;nb=1;a=&c;return *this;}
 basicAC_F0 & operator=(pair<int,C_F0*> p)  {named_parameter=0;nb=p.first;a=p.second;return *this;} 
 const C_F0 & operator [] (int i) const {throwassert(a && i<nb);return a[i];}
 int size() const {return nb;}
 C_F0 * ptr() const  {return a;}
 C_F0  find(const char * k) const  { 
    assert(k);  
   if (named_parameter) { const_iterator i=named_parameter->find(k) ;
    if (i == named_parameter->end() ) return C_F0();
    else  return i->second;}
   else return C_F0();} 
   
 struct name_and_type{ 
  const char * name;
  const type_info * type;
 } ;
 
  void SetNameParam(int n=0,name_and_type *l=0 , Expression * e=0) const ;
};


class AC_F0: public basicAC_F0 { //  a Array of C_F0
//    tableau d'un parametres  max 1024 parametres 
//    avec allocation  
 const static  int MaxSize;
 //  no constructor in this class (this class is in a union )
  public:
  AC_F0 & operator=(pair<const char *,const C_F0> p) {
    named_parameter=0; a=new C_F0[MaxSize]; nb=0;Add(p.first,p.second);return *this;}
  AC_F0 & operator+=(pair<const char *,const C_F0> p) {Add(p.first,p.second);return *this;} 

  AC_F0 & operator=(long k) {throwassert(k==0);named_parameter=0;a=new C_F0[MaxSize]; nb=0;return *this;}     
  AC_F0 & operator=(const C_F0& c) {named_parameter=0; a=new C_F0[MaxSize]; nb=0;a[nb++]=c;return *this;} 
  AC_F0 & operator+=(const C_F0& c) { 
       if ( ! (a&& nb<MaxSize))
         CompileError("Sorry  number of parameters > 1024");
       a[nb++]=c;return *this;} 
  AC_F0 & Add(const char * nm,const C_F0 &c)  {
     if (!named_parameter) named_parameter=new maptype();
    iterator i=named_parameter->find(nm);
    if(i==named_parameter->end()) named_parameter->insert(make_pair(nm,c));
    else {cerr << " the named in the list already exists "<< nm <<endl; CompileError();}
    return *this;}
  int size() const {return nb;}
  const C_F0 & operator [] (int i) const {throwassert(a && i<nb);return a[i];}
  void destroy() {
      
        nb=0;       
        if(named_parameter) 
          delete named_parameter;
        if (a)
           delete []  a;
        a=0;named_parameter=0;}
  
}; 

class  basicAC_F0_wa : public basicAC_F0 { public:
 basicAC_F0_wa(const C_F0 & e) {
   named_parameter=0;
   nb=1;
   a= new C_F0[nb];
   a[0]=e;
 }
 basicAC_F0_wa(const C_F0 & e,const C_F0 & ee) {
   named_parameter=0;
   nb=2;
   a= new C_F0[nb];
   a[0]=e;
   a[1]=ee;
 }
 

 basicAC_F0_wa(C_F0 e,const basicAC_F0 & b) { 
   named_parameter=0;
   if (b.named_parameter) named_parameter = new maptype(*b.named_parameter);
   nb=1+b.size();
   a= new C_F0[nb];
   a[0]=e;
   for (int i=1;i<nb;i++) a[i]=b[i-1];}
 ~basicAC_F0_wa(){delete [] a;} 

 basicAC_F0_wa(const basicAC_F0 & b) { 
   named_parameter=0;
   if (b.named_parameter) named_parameter = new maptype(*b.named_parameter);
   nb=b.size();
   a= new C_F0[nb];
   for (int i=0;i<nb;i++) a[i]=b[i];}
   
  private: 
   void operator=(const basicAC_F0 & b);
};



class E_Array  :public E_F0 {  public: 
  basicAC_F0_wa *v;// the value
  E_Array(const basicAC_F0 & aa) : v(new basicAC_F0_wa(aa))  {throwassert(v);}
  AnyType operator()(Stack)  const {
     cerr << " No evaluation of an E_array" << endl;
     throwassert(0);
     return  Nothing;}
 const C_F0 & operator [] (int i) const {throwassert(v );return (*v)[i];}
 int size() const {return v->size();}
 size_t nbitem() const {return v->size();}
 void map(C_F0 (*mapping)(const C_F0 & )) const 
   { for (int i=0;i<v->size();i++) 
      v->a[i]=(*mapping)(v->a[i]);}
  virtual bool MeshIndependent() const {
    for (int i=0;i<v->size();i++) 
      if (v->a[i].MeshIndependent()) return false;
     return false;
   
  } // 
  operator aType () const { return atype<void>();} 
 
  };
class PlotStream;  
class E_Border ;
class E_BorderN :public E_F0mps { public: 
   const E_Border * b;
   Expression  n;
   const E_BorderN * next;
   E_BorderN(const E_Border *  bb, C_F0  nn,const E_BorderN * nx=0) ; 
    E_BorderN(const E_BorderN & bb,const E_BorderN * nx)
    : b(bb.b),n(bb.n),next(nx)
    {
      int kk=1;
	if(bb.next) {// modif FH. 13/02/2008
	       
	    const E_BorderN ** pnext = &next;
	    E_BorderN  *pp;
	    next = bb.next;  // copy bb;       
	    for(int step=0;step<2;++step)
	    {
	      while (*pnext)
	       {
		 kk++;
		 pp = new E_BorderN(**pnext); // copy 
		 *pnext = pp;
		 pnext = & pp->next;
	       }
	      if(step==0)
	       *pnext= nx;  // copy de nx
	    }
	   // cout << "  BorderN : nb item : " << kk << " == " << size()<< endl;  
	}
    }   
    AnyType operator()(Stack)  const {
     return  SetAny<const  E_BorderN *>(this);}  
  operator aType () const { return atype<const  E_BorderN *>();}         
    int size() const { int k=0;for(const E_BorderN  *pp=this;pp; pp=pp->next) k++; return k;}  
   E_BorderN * operator+( const E_BorderN & bb)  const 
   { throwassert(bb.next==0);
     return new E_BorderN(bb,this);}
  long Nbseg(Stack stack) const { return GetAny<long>((*n)(stack));}
  double from(Stack stack) const ;//{ return GetAny<double>((*n)(stack));}
  double to(Stack stack) const ;//{ return GetAny<double>((*b)(stack));}
  double * var(Stack stack) const ;//{ return GetAny<double*>((*n)(stack));}
  void code(Stack stack) const ;
  long label()const  ;
  void Plot(Stack stack) const ;
  void SavePlot(Stack stack,PlotStream & plot ) const;
    
  void BoundingBox(Stack stack,double  &xmin,double & xmax, double & ymin,double & ymax) const ;
};

class AddBorderOperator: public  OneOperator{
  typedef const E_BorderN * A;
    public: 
    E_F0 * code(const basicAC_F0 & args) const 
     {
       A a0=dynamic_cast<A>((Expression) args[0]);
       A a1=dynamic_cast<A>((Expression) args[1]);
       ffassert( a0 && a1);
       //ffassert(a1->next==0);  // change FH 13/02/2008
       return  new E_BorderN(*a1,a0);} 
    AddBorderOperator(): 
      OneOperator(map_type[typeid(A).name()],map_type[typeid(A).name()],map_type[typeid(A).name()])
      {pref = 0;}

};


class  OneOperator_borderN : public OneOperator {public:
  const  E_Border * theborder;
    E_F0 * code(const basicAC_F0 & a) const 
     { return  new E_BorderN(theborder,a[0]);} 
    OneOperator_borderN(const  E_Border * b)
      : OneOperator(atype<const E_BorderN *>(),atype<long>()),
      theborder(b){}
};

class E_Border  :public Polymorphic  {  public: 
  static basicAC_F0::name_and_type name_param[] ;
  static const int n_name_param =0;
  static long Count;
  Expression xvar,xfrom,xto,xcode;
  basicAC_F0_wa * tab;
  long label;
  E_Border(const E_Array * a) : 
    xvar(0),xfrom(0),xto(0),xcode(0), tab(a? a->v:0) ,label(++Count) 
  {
    assert(tab); 
    Add("(",new OneOperator_borderN(this));
  }
  
  E_Border(const basicAC_F0 & aa) :    
    xvar(to<double*>(aa[0])),
    xfrom(to<double>(aa[1])),
    xto(to<double>(aa[2])),
    xcode(aa[3].LeftValue()),
    tab(0),
    label(++Count)
  {    Add("(",new OneOperator_borderN(this));}

  AnyType operator()(Stack)  const {
    return  SetAny<const  E_Border *>(this);}
  double length(Stack ) const { ffassert(0);return 0.0; /* a faire */ }
};
  
inline  E_BorderN::E_BorderN(const E_Border *bb, C_F0  nn,const E_BorderN * nx)  
     :b(bb),n(::to<long>(nn)),next(nx) { throwassert(b);}

inline  double E_BorderN::from(Stack stack) const { return b->xfrom ? GetAny<double>((*b->xfrom)(stack)): double(0.0);}
inline  double  E_BorderN::to(Stack stack) const { return b->xto? GetAny<double>((*b->xto)(stack)): b->length(stack) ;}
inline  double *  E_BorderN::var(Stack stack) const { return b->xvar ? GetAny<double*>((*b->xvar)(stack)): (double*) 0 ;}
inline  void  E_BorderN::code(Stack stack)const { (*b->xcode)(stack);}
inline  long  E_BorderN::label()const { return b->label;}

inline ArrayOfaType::ArrayOfaType(const basicAC_F0 & aa) : n(aa.size()),t(n ? (n<=4 ? tt : new aType[n]):0),ellipse(false) { 
   for (int i=0;i<n;i++) t[i]=aa[i].left();}
   
inline ArrayOfaType::ArrayOfaType(const ArrayOfaType & aa) : n(aa.n),t(n<=4?tt:new aType[n]),ellipse(aa.ellipse) { 
   for (int i=0;i<n;i++) t[i]=aa.t[i];}   


inline C_F0 TableOfIdentifier::Find(const char * name) const  {
    const_iterator i=m.find(name); 
    if ( i == m.end()) { return C_F0();}
    else return C_F0(i->second);}

inline C_F0 TableOfIdentifier::Find(const char * name,const basicAC_F0 & args) const  {
    const_iterator i=m.find(name); 
    if ( i == m.end()) {cerr<<"No operator " << name<<endl;
    cerr <<*this << endl;CompileError("TableOfIdentifier::Find");return C_F0();}
    else {return C_F0(C_F0(i->second),"(",args);}}
//  Attention il y a moralement un bug
//  les initialisation   x = y   ( passe par l'operateur binaire <-  dans TheOperators
//   les initialisation   x(y)   ( passe par l'operateur unaire <-  du type de x
//   -------

inline size_t align8(size_t &off) 
{ 
 register size_t o= off %8 ;
  off += o ? 8-o : 0;
 return off;
}


template<class T>
inline Type_Expr  NewVariable(aType t,size_t &off) 
{ 
   size_t o= align8(off);//  align    
 //  off += t->un_ptr_type->size;
 // bug    off += t->size;
   off += t->un_ptr_type->size; // correction 16/09/2003 merci