lgfem.cpp

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

}; 
 class E_P_Stack_inside   :public  E_F0mps { public: 
  AnyType operator()(Stack s)  const { throwassert(* (long *) s);
    return SetAny<double>(MeshPointStack(s)->outside? 0.0 : 1.0 );} 
    operator aType () const { return atype<double>();}         
    
}; 

class E_P_Stack_lenEdge   :public  E_F0mps { public: 
  AnyType operator()(Stack s)  const { throwassert(* (long *) s);
    MeshPoint * mp=MeshPointStack(s);      
    ffassert(mp->T && mp ->e >=0 && mp->d==2);
    double l= mp->T->lenEdge(mp->e);
    return SetAny<double>(l);} 
    operator aType () const { return atype<double>();}         
    
}; 

class E_P_Stack_hTriangle   :public  E_F0mps { public: 
  AnyType operator()(Stack s)  const { throwassert(* (long *) s);
    MeshPoint * mp=MeshPointStack(s);
    assert(mp->T) ;
    double l= mp->T->h();
    return SetAny<double>(l);} 
    operator aType () const  { return atype<double>();}         
    
};

class E_P_Stack_nTonEdge   :public  E_F0mps { public: 
  AnyType operator()(Stack s)  const { throwassert(* (long *) s);
    MeshPoint * mp=MeshPointStack(s);
    assert(mp->T && mp->e > -1 && mp->d==2 ) ;
    long l=mp->Th->nTonEdge(mp->t,mp->e);
    // cout << " nTonEdge " << l << endl;
    return SetAny<long>( l) ;} 
    operator aType () const  { return atype<long>();}         
    
}; 

class E_P_Stack_areaTriangle   :public  E_F0mps { public: 
  AnyType operator()(Stack s)  const { throwassert(* (long *) s);
    MeshPoint * mp=MeshPointStack(s);
    assert(mp->T) ;
      double l=-1; // unset ...
    if(mp->d==2)	
      l= mp->T->area;	
    else if (mp->d==3 && mp->f >=0)
      {	
	  R3 NN = mp->T3->N(mp->f);
	  l= NN.norme()/2.;
      }
    else 
      {
	  cout << "erreur : E_P_Stack_areaTriangle" << mp->d << " " << mp->f << endl;
	  ffassert(0); // undef 
      }
    return SetAny<double>(l);} 
    operator aType () const { return atype<double>();}         
    
}; 

/*
class New_MeshPoint : public E_F0mps {

};*/


 

 
//inline pfes MakePtr(pmesh * const &  a, TypeOfFE * const & tef){ return pfes(new pfes_tef(a,tef)) ;}
//inline pfes MakePtr(pmesh * const &  a){ return pfes(new pfes_tef(a,&P1Lagrange)) ;}
//inline pfes MakePtr(pfes * const &  a,long const & n){ return pfes(new pfes_fes(a,n)) ;}
/*
class E_pfes : public E_F0 {
  const int N;
  Expression *Th,*tef; 
  E_pfes(Expression *TTh,*ttef) : Th(TTh),tef(ttef),N(ttef?ttef->nbitem():0) {}
   virtual AnyType operator()(Stack)  const {
   return AnyType<pfes*>
  }
  
  virtual size_t nbitem() const {return 1;}  
};
  OneOperator_pfes():OneOperator(atype<void>(),atype<E_Array>(),atype<E_Array>()) {}
    E_F0 * code(const basicAC_F0 & args) const ; 
  
};*/

template<class R>
class LinearCG : public OneOperator 
{ public:
  typedef KN<R> Kn;
  typedef KN_<R> Kn_;
  const int cas;

 class MatF_O: VirtualMatrice<R> { public:
   Stack stack;
   mutable  Kn x;
   C_F0 c_x;
   Expression  mat;
   typedef  typename VirtualMatrice<R>::plusAx plusAx;
   MatF_O(int n,Stack stk,const OneOperator * op) 
     : VirtualMatrice<R>(n),stack(stk),
       x(n),c_x(CPValue(x)),
       mat(op->code(basicAC_F0_wa(c_x))) {
         ffassert(atype<Kn >() ==(aType) *op);
         // WhereStackOfPtr2Free(stack)=new StackOfPtr2Free(stack);// FH mars 2005   
         
         }
   ~MatF_O() { 
     // cout << " del MatF_O mat " << endl;
     delete mat;
     // cout << " del MatF_Ocx ..." <<  endl;
      Expression zzz = c_x;
     // cout << " zzz "<< zzz << endl;
     delete zzz;
    // WhereStackOfPtr2Free(stack)->clean(); // FH mars 2005 
     
     }
   void addMatMul(const  Kn_  & xx, Kn_ & Ax) const { 
      ffassert(xx.N()==Ax.N());
      x =xx;
      Ax  += GetAny<Kn>((*mat)(stack));
      WhereStackOfPtr2Free(stack)->clean();
       } 
    plusAx operator*(const Kn &  x) const {return plusAx(this,x);} 
  virtual bool ChecknbLine(int n) const { return true;}  
  virtual bool ChecknbColumn(int m) const { return true;} 
    
};  
 

  class E_LCG: public E_F0mps { public:
   const int cas;// <0 => Nolinear
   static const int n_name_param=4;

   static basicAC_F0::name_and_type name_param[] ;


  Expression nargs[n_name_param];
   
  const OneOperator *A, *C; 
  Expression X,B;

  E_LCG(const basicAC_F0 & args,int cc) :cas(cc)
   {
      args.SetNameParam(n_name_param,name_param,nargs);
      {  const  Polymorphic * op=  dynamic_cast<const  Polymorphic *>(args[0].LeftValue());
         ffassert(op);
         A = op->Find("(",ArrayOfaType(atype<Kn* >(),false)); }
      if (nargs[2]) 
      {  const  Polymorphic * op=  dynamic_cast<const  Polymorphic *>(nargs[2]);
         ffassert(op); 
         C = op->Find("(",ArrayOfaType(atype<Kn* >(),false)); }
       else  C =0;
      X = to<Kn*>(args[1]);
      if (args.size()>2)
        B = to<Kn*>(args[2]);
      else 
        B=0;
   }
     
     virtual AnyType operator()(Stack stack)  const {
       int ret=-1;

      // WhereStackOfPtr2Free(stack)=new StackOfPtr2Free(stack);// FH mars 2005   
      try {
      Kn &x = *GetAny<Kn *>((*X)(stack));
      int n=x.N();
      MatF_O AA(n,stack,A);
      double eps = 1.0e-6;
      int nbitermax=  100;
      if (nargs[0]) eps= GetAny<double>((*nargs[0])(stack));
      if (nargs[1]) nbitermax = GetAny<long>((*nargs[1])(stack));
      if (nargs[3]) eps= *GetAny<double*>((*nargs[3])(stack));
     
      KN<R>  bzero(B?1:n); // const array zero
      bzero=R(); 
      KN<R> *bb=&bzero; 
      if (B) {
        Kn &b = *GetAny<Kn *>((*B)(stack));
        R p = (b,b);
       if (p) 
         {
          // ExecError("Sorry LinearCG work only with nul right hand side, so put the right hand in the function");
          }
         bb = &b;
      }
      if (cas<0) {
       if (C) 
         { MatF_O CC(n,stack,C);
           ret = NLCG(AA,CC,x,nbitermax,eps, 51L-Min(Abs(verbosity),50L) );}
        else 
           ret = NLCG(AA,MatriceIdentite<R>(),x,nbitermax,eps, 51L-Min(Abs(verbosity),50L));
        }
      else 
      if (C) 
       { MatF_O CC(n,stack,C);
         ret = ConjuguedGradient2(AA,CC,x,*bb,nbitermax,eps, 51L-Min(Abs(verbosity),50L) );}
      else 
         ret = ConjuguedGradient2(AA,MatriceIdentite<R>(),x,*bb,nbitermax,eps, 51L-Min(Abs(verbosity),50L));
      if( nargs[3]) *GetAny<double*>((*nargs[3])(stack)) = -(eps);
      }
      catch(...)
      {
       // WhereStackOfPtr2Free(stack)->clean(); // FH mars 2005 
        throw;
      }
     // WhereStackOfPtr2Free(stack)->clean(); // FH mars 2005 
      
      return SetAny<long>(ret);
       
     }  
    operator aType () const { return atype<long>();}         
     
  };
  
  E_F0 * code(const basicAC_F0 & args) const {
    return new E_LCG(args,cas);}
  LinearCG() :   OneOperator(atype<long>(),
                             atype<Polymorphic*>(),
                             atype<KN<R> *>(),atype<KN<R> *>()),cas(2){}
  LinearCG(int cc) :   OneOperator(atype<long>(),
                             atype<Polymorphic*>(),
                             atype<KN<R> *>()),cas(cc){}
   
};


template<class R>
basicAC_F0::name_and_type  LinearCG<R>::E_LCG::name_param[]= {
  {   "eps", &typeid(double)  },
  {   "nbiter",&typeid(long) },
  {   "precon",&typeid(Polymorphic*)},
  {   "veps" ,  &typeid(double*) }
};


template<class R>
class LinearGMRES : public OneOperator 
{ public:
  typedef KN<R> Kn;
  typedef KN_<R> Kn_;
  const int cas;

 class MatF_O: VirtualMatrice<R> { public:
   Stack stack;
   mutable  Kn x;
   C_F0 c_x;
   Expression  mat;
   typedef  typename VirtualMatrice<R>::plusAx plusAx;
   MatF_O(int n,Stack stk,const OneOperator * op) 
     : VirtualMatrice<R>(n),
       stack(stk),
       x(n),c_x(CPValue(x)),
       mat(op->code(basicAC_F0_wa(c_x))) {
       ffassert(atype<Kn >() ==(aType) *op); }
   ~MatF_O() { delete mat;delete c_x.LeftValue();}
   void addMatMul(const  Kn_  & xx, Kn_ & Ax) const { 
      ffassert(xx.N()==Ax.N());
      x =xx;
      Ax  += GetAny<Kn>((*mat)(stack));
      WhereStackOfPtr2Free(stack)->clean(); //  add dec 2008 
   } 
    plusAx operator*(const Kn &  x) const {return plusAx(this,x);} 
  virtual bool ChecknbLine(int n) const { return true;}  
  virtual bool ChecknbColumn(int m) const { return true;} 
    
};  
 

  class E_LGMRES: public E_F0mps { public:
   const int cas;// <0 => Nolinear
   static basicAC_F0::name_and_type name_param[] ;
   static const int n_name_param =5;
   Expression nargs[n_name_param];
  const OneOperator *A, *C; 
  Expression X,B;

  E_LGMRES(const basicAC_F0 & args,int cc) :cas(cc)
   {
      args.SetNameParam(n_name_param,name_param,nargs);
      {  const  Polymorphic * op=  dynamic_cast<const  Polymorphic *>(args[0].LeftValue());
         ffassert(op);
         A = op->Find("(",ArrayOfaType(atype<Kn* >(),false)); }
      if (nargs[2]) 
      {  const  Polymorphic * op=  dynamic_cast<const  Polymorphic *>(nargs[2]);
         ffassert(op); 
         C = op->Find("(",ArrayOfaType(atype<Kn* >(),false)); }
       else  C =0;
      X = to<Kn*>(args[1]);
      if (args.size()>2)
        B = to<Kn*>(args[2]);
      else 
        B=0;
   }
     
     virtual AnyType operator()(Stack stack)  const {
      Kn &x = *GetAny<Kn *>((*X)(stack));
      Kn b(x.n);
     
      if (B)   b = *GetAny<Kn *>((*B)(stack));
      else     b=0.;
      int n=x.N();
      int dKrylov=50;
      MatF_O AA(n,stack,A);
      double eps = 1.0e-6;
      int nbitermax=  100;
      if (nargs[0]) eps= GetAny<double>((*nargs[0])(stack));
      if (nargs[1]) nbitermax = GetAny<long>((*nargs[1])(stack));
      if (nargs[3]) eps= *GetAny<double*>((*nargs[3])(stack));
      if (nargs[4]) dKrylov= GetAny<long>((*nargs[4])(stack));
      
      int ret;
      if(verbosity>4)
        cout << "  ..GMRES: eps= " << eps << " max iter " << nbitermax 
             << " dim of Krylov space " << dKrylov << endl;
        KNM<R> H(dKrylov+1,dKrylov+1);
	int k=dKrylov;//,nn=n;
       double epsr=eps;
      // int res=GMRES(a,(KN<R> &)x, (const KN<R> &)b,*this,H,k,nn,epsr);
      if (cas<0) {
        ErrorExec("NL GMRES:  to do! sorry ",1);
/*       if (C) 
         { MatF_O CC(n,stack,C);
           ret = NLGMRES(AA,CC,x,nbitermax,eps, 51L-Min(Abs(verbosity),50L) );}
        else 
           ret = NLGMRES(AA,MatriceIdentite<R>(),x,nbitermax,eps, 51L-Min(Abs(verbosity),50L));
         ConjuguedGradient  */
        }
      else 
       {
       if (C)
        { MatF_O CC(n,stack,C); 
         ret=GMRES(AA,(KN<R> &)x, (const KN<R> &)b,CC,H,k,nbitermax,epsr);}
       else
         ret=GMRES(AA,(KN<R> &)x, (const KN<R> &)b,MatriceIdentite<R>(),H,k,nbitermax,epsr);       
       }
       /*
      if (C) 
       { MatF_O CC(n,stack,C);
         
         ret = ConjuguedGradient2(AA,CC,x,nbitermax,eps, 51L-Min(Abs(verbosity),50L) );}
      else 
         ret = ConjuguedGradient2(AA,MatriceIdentite<R>(),x,nbitermax,eps, 51L-Min(Abs(verbosity),50L));*/
         
     // if( nargs[3]) *GetAny<double*>((*nargs[3])(stack)) = -(eps);
      return SetAny<long>(ret);
       
     }  
    operator aType () const { return atype<long>();}         
     
  };
  
  E_F0 * code(const basicAC_F0 & args) const {
    return new E_LGMRES(args,cas);}
  LinearGMRES() :   OneOperator(atype<long>(),
                             atype<Polymorphic*>(),
                             atype<KN<R> *>(),atype<KN<R> *>()),cas(2){}
  LinearGMRES(int cc) :   OneOperator(atype<long>(),
                             atype<Polymorphic*>(),
                             atype<KN<R> *>()),cas(cc){}
   
};


template<class R>
basicAC_F0::name_and_type  LinearGMRES<R>::E_LGMRES::name_param[]= {
  {   "eps", &typeid(double)  },
  {   "nbiter",&typeid(long) },
  {   "precon",&typeid(Polymorphic*)},
  {   "veps" ,  &typeid(double*) },
  {   "dimKrylov", &typeid(long) }
};

template<typename int2>
typename map<int,int2>::iterator closeto(map<int,int2> & m, int k)
 {
  typename map<int,int2>::iterator i=  m.find(k);
   if (i==m.end()) 
    {
     i=  m.find(k+1);
     if (i==m.end()) 
      i=  m.find(k-1);
    }
   return i;
 }
 
template<class T,int N>
class Smallvect { public: 
 T v[N];
 T & operator[](int i){return v[i];}
 const T & operator[](int i) const {return v[i];}
};
template<class T,int N>
ostream & operator<<(ostream & f,const Smallvect<T,N> & v)
{
    for(int i=0;i<N;++i)  f << v[i] << ' ';
    return f;
}

template<class T> 
int numeroteclink(KN_<T> & ndfv) 
{
         int nbdfv =0;
         for (int i=0;i<ndfv.N();i++)
           if (ndfv[i]>=i)
           {
             int j=i,ii,kkk=0;
             
             do {
               ii=ndfv[j];
		 ffassert(kkk++<10);
	       assert(nbdfv <= j);