doperator.hpp

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

 {LinearComb<I,R> r(a);r*=b;return r;}
template<class I,class R>
LinearComb<I,R> operator*(const R & b,const LinearComb<I,R> & a)
 {LinearComb<I,R> r(a);r*=b;return r;}
 

 
class MGauche :public pair<int,operatortype> {public:
  MGauche() {}
  MGauche(int i,operatortype j) {first = i;second= j;}
  MGauche(const pair<int,operatortype> &p) : pair<int,operatortype>(p){}
  bool operator==(const MGauche& a) const {
    return static_cast<bool>(first == a.first && second == a.second);}
  int maxop(int op) const { return Max(op,(int) second);}
    
};

class MDroit :public pair<int,operatortype> {public:
  MDroit(){}
  MDroit(int i,operatortype j) {first = i;second =j;}
  //   first : number of unknow 
  // second  : number of operator ( 0 Id, 1 dx, 2 dy)
  MDroit(const pair<int,operatortype> &p) : pair<int,operatortype>(p){}
  bool operator==(const MDroit& a) const {
    return static_cast<bool>(first == a.first && second == a.second);}
};

inline ostream & operator<<(ostream & f,const MDroit & p)
{ f << p.first <<','<<p.second ;
  return f;}
inline ostream & operator<<(ostream & f,const MGauche & p)
{ f << p.first <<';'<<p.second ;
  return f;}
  
template<typename A,typename B>
inline ostream & operator<<(ostream & f,const pair<A,B> &p)
{ f << p.first <<" "<<p.second ;
  return f;}
  
  
//extern const C_F0 & One, &Zero;
C_F0 operator*(const C_F0 &,const C_F0 &);
C_F0 & operator+=(const C_F0 &,const C_F0 &);
 
typedef LinearComb<MGauche,C_F0> LinearOperatorG;
typedef LinearComb<MDroit ,C_F0> LinearOperatorD;

typedef LinearComb<pair<MGauche,MDroit>,C_F0> BilinearOperator;

inline  int maxop(int op,const MGauche & v) 
    { return Max(op,(int) v.second);}
inline  int maxop(int op,const MDroit & v) 
    { return Max(op,(int) v.second);}
inline  int maxop(int op,const pair<MGauche,MDroit> & b) 
    { return Max(op,(int) b.first.second,(int) b.second.second );}
 
inline BilinearOperator operator*(const LinearOperatorG & a,const LinearOperatorD & b) 
 {
   BilinearOperator r;
   for (LinearOperatorG::const_iterator i=a.v.begin();i!=a.v.end();i++)
     for (LinearOperatorD::const_iterator j=b.v.begin();j!=b.v.end();j++)
       {
	 const LinearOperatorG::K  vi(*i);
	 const LinearOperatorD::K  vj(*j);
	
	 r.add(make_pair(vi.first,vj.first),vi.second*vj.second);
	   
       }
    return r;
 } 
inline BilinearOperator operator*(const LinearOperatorD & b,const LinearOperatorG & a) 
 {
   BilinearOperator r;
   for (LinearOperatorG::const_iterator i=a.v.begin();i!=a.v.end();i++)
     for (LinearOperatorD::const_iterator j=b.v.begin();j!=b.v.end();j++)
       {
	 const LinearOperatorG::K  vi(*i);
	 const LinearOperatorD::K  vj(*j);
	
	 r.add(make_pair(vi.first,vj.first),vi.second*vj.second);
	   
       }
    return r;
 } 
ostream & operator<<(ostream & f,const  BilinearOperator & a); 
 
inline ostream & operator<<(ostream & f,const  BilinearOperator & a)
{
  int k=0;
  for (BilinearOperator::const_iterator i=a.v.begin();i!=a.v.end();i++)
    {
      const BilinearOperator::K  vi(*i);
      const char * www[]={" ","_x ","_y "};
      const  pair<pair<int,int>,pair<int,int> > i1(vi.first);
      const  pair<int,int> ii(i1.first),jj(i1.second);
      f << *(const E_F0 *) vi.second <<  char('u'+ii.first) << www[ii.second] << " " << char('u'+jj.first)<<"'" << www[jj.second] ;
      if ( (++k%5)==0) f << endl ; else f << " ";
    }
  return f;
}
typedef LinearOperatorD LOperaD;
typedef LinearOperatorG LOperaG;
typedef BilinearOperator Opera;
inline LOperaG DiffG(int i,operatortype j) { return LOperaG(make_pair(i,j),*pOne);}
inline LOperaD DiffD(int i,operatortype j) { return LOperaD(make_pair(i,j),*pOne);}

inline LOperaG *newU_(int i) {  LOperaG * r;
  r= new LOperaG(make_pair(i,op_id),*pOne);
 SHOWVERB( cout << "newU_ " << r << endl);
  return r; }
inline LOperaG *newU_x(int i) { return new LOperaG(make_pair(i,op_dx),*pOne);}
inline LOperaG *newU_y(int i) { return new LOperaG(make_pair(i,op_dy),*pOne);}
inline LOperaD *newV_(int i) { return new LOperaD(make_pair(i,op_id),*pOne);}
inline LOperaD *newV_x(int i) { return new LOperaD(make_pair(i,op_dx),*pOne);}
inline LOperaD *newV_y(int i) { return new LOperaD(make_pair(i,op_dy),*pOne);}

template<class L>
 L *  Diff(const L *   u,const  operatortype & d) { 
     throwassert(u);
     L * r= new  L(*u);
    for (typename L::iterator i=r->v.begin();i!=r->v.end();i++)
      {   operatortype & dd=i->first.second;
          
         if (dd != op_id)
           { ffassert(0); i->first.second = op_id;  } // a faire          
         else {
           ffassert(i->second.EvaluableWithOutStack());// a faire derivation des fonctions
           dd = d ; }
    }
    return r;}
    


template<class L>
int  MaxOp(const L *   u) { 
     throwassert(u);
     int op=0;
    for (typename L::const_iterator i=u->v.begin();i!=u->v.end();i++)          
          op = maxop(op,i->first);
    return op;}
    
    
    
template<class L,class A,L* ff(const basicAC_F0 & args)>
class CODE_L1 { public:
  typedef L* Result;
  typedef L* (*func)(const basicAC_F0 & args) ;
  static L* f(const basicAC_F0 & args) { return ff(args);} //  ff :A-> L*
  static ArrayOfaType  typeargs() {return ArrayOfaType(atype<A>);}
};

template<class L>
class CODE_L_Add { public:
  typedef const L* Result;
  static  E_F0 * f(const basicAC_F0 & args) { 
    const L * a(dynamic_cast<const L*>((Expression) args[0]));
    const L * b(dynamic_cast<const L*>((Expression) args[1]));
    throwassert(a  && b );
    return new L(*a+*b);}
   static ArrayOfaType  typeargs() {return ArrayOfaType(atype<const L*>(),atype<const L*>());}
};
template<class L>
class CODE_L_Sub { public:
  typedef const L* Result;
  static  E_F0 * f(const basicAC_F0 & args) { 
    const L * a(dynamic_cast<const L*>((Expression) args[0]));
    const L * b(dynamic_cast<const L*>((Expression) args[1]));
    throwassert(a  && b );
    L * bb = new L(*pminusOne * *b);
    return new L(*a+*bb);}
   static ArrayOfaType  typeargs() {return ArrayOfaType(atype<const L*>(),atype<const L*>());}
};

template<class L>
class CODE_L_Minus { public:
  typedef const L* Result;
  static  E_F0 * f(const basicAC_F0 & args) { 
    const L * a(dynamic_cast<const L*>((Expression) args[0]));
    throwassert(a  && pminusOne );
    return new L(*pminusOne * *a);}
   static ArrayOfaType  typeargs() {return ArrayOfaType(atype<const L*>());}
};
template<class L>
class CODE_L_Plus { public:
  typedef const L* Result;
  static  E_F0 * f(const basicAC_F0 & args) { 
    const L * a(dynamic_cast<const L*>((Expression) args[0]));
    throwassert(a );
    return new L( *a);}
   static ArrayOfaType  typeargs() {return ArrayOfaType(atype<const L*>());}
};

template<class L,class A,class B>
class CODE_L_Mul { public:
  typedef const L* Result;
  static  E_F0 * f(const basicAC_F0 & args) { 
    const A * a(dynamic_cast<const A*>((Expression) args[0]));
    const B * b(dynamic_cast<const B*>((Expression) args[1]));
    SHOWVERB(cout << " CODE_L_Mul " << a << " " << b << endl);
    throwassert(a  && b );
    return new L(*a * *b);}
   static ArrayOfaType  typeargs() {return ArrayOfaType(atype<const A*>(),atype<const B*>());}
};

template<class T,class L>
class CODE_L_MulRL { public:
  typedef const L* Result;
  static  E_F0 * f(const basicAC_F0 & args) { 
    const L * b(dynamic_cast<const L*>((Expression) args[1]));
    if (!b) {
      
      cout << " --- " << ((Expression) args[1]) << typeid((Expression) args[1]).name() << endl;
    throwassert( b );
    }
    return new L(to<T>(args[0]) * *b);}
   static ArrayOfaType  typeargs() {return ArrayOfaType(atype<T>(),atype<const L*>());}
};

template<class L,class T>
class CODE_L_MulLR { public:
  typedef const L* Result;
  static  E_F0 * f(const basicAC_F0 & args) { 
    const L * a(dynamic_cast<const L*>((Expression) args[0]));
    throwassert( a );
    return new L(to<T>(args[1]) * *a);}
   static ArrayOfaType  typeargs() {return ArrayOfaType(atype<Result>(),atype<T>());}
};

template<class L,class T>
class CODE_L_DivLR { public:
  typedef const L* Result;
  static  E_F0 * f(const basicAC_F0 & args) { 
    const L * a(dynamic_cast<const L*>((Expression) args[0]));
    throwassert( a );
    return new L(C_F0(TheOperators,"/",*pOne,args[1]) * *a);}
   static ArrayOfaType  typeargs() {return ArrayOfaType(atype<Result>(),atype<T>());}
};

template<class L,operatortype op>
class CODE_Diff { public:
  typedef const L* Result;
  static  E_F0 * f(const basicAC_F0 & args) { 
   const L * a=dynamic_cast<const L*>((Expression) args[0]) ;
    throwassert( a );
    return  Diff(a,op);}
   static ArrayOfaType  typeargs() {return ArrayOfaType(atype<Result>());}
};

enum TheCode_VF { Code_Jump=1, Code_Mean=2, Code_OtherSide=3};

template<class L,TheCode_VF code>
class Code_VF { public:  
  typedef const L* Result;
  static  E_F0 * f(const basicAC_F0 & args) { 
   const L * u=dynamic_cast<const L*>((Expression) args[0]) ;
    assert( u );
     L * r= new  L(*u);
    for (typename L::iterator i=r->v.begin();i!=r->v.end();i++)
      {   operatortype & dd=i->first.second;
          assert(dd<last_operatortype);
          dd = (operatortype) ((int) dd+last_operatortype*code) ; 
     }    
    return r;}
   static ArrayOfaType  typeargs() {return ArrayOfaType(atype<Result>());}
};