rnm.hpp

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

        { KNM_<R>::operator=(u);}
        
  ~KNM(){delete [] this->v;}
  
   KNM& operator=(const KNM_<const_R> & u)   
        { KNM_<R>::operator=(u);return *this;}
   KNM& operator=(const_R a)                 
        { KNM_<R>::operator=(a);return *this;}
   KNM& operator+=(const_R  a)               
        { KNM_<R>::operator+=(a);return *this;}
   KNM& operator-=(const_R  a)               
        { KNM_<R>::operator-=(a);return *this;}
   KNM& operator/=(const_R  a)               
        { KNM_<R>::operator/=(a);return *this;}
   KNM& operator*=(const_R  a)               
        { KNM_<R>::operator*=(a);return *this;}
   KNM& operator+=(const KNM_<const_R> & u)  
        { KNM_<R>::operator+=(u);return *this;}
   KNM& operator-=(const KNM_<const_R> & u)  
        { KNM_<R>::operator-=(u);return *this;}

   KNM& operator/=(const KNM_<const_R> & u)  
        { KNM_<R>::operator/=(u);return *this;}
   KNM& operator*=(const KNM_<const_R> & u)  
        { KNM_<R>::operator*=(u);return *this;}


   KNM &operator  =(const outProduct_KN_<R> & u)
        { KNM_<R>::operator =(u);return *this;}
   KNM &operator +=(const outProduct_KN_<R> & u)
        { KNM_<R>::operator+=(u);return *this;}
   KNM &operator -=(const outProduct_KN_<R> & u)
        { KNM_<R>::operator-=(u);return *this;}
   KNM &operator /=(const outProduct_KN_<R> & u) 
        { KNM_<R>::operator/=(u);return *this;}
   KNM &operator *=(const outProduct_KN_<R> & u)
        { KNM_<R>::operator*=(u);return *this;}
  
        
  //  two opertors to cast to un array of constant        
//    operator KNM_<const_R> & ()  
//          { return *  (KNM_<const_R>*) this;}
//    operator KNM_<const_R> const & ()  const 
//          { return *(const KNM_<const_R>*) this;}

//    operator KNM<const_R> & ()  
//          { return *  (KNM<const_R>*) this;}
//    operator KNM<const_R> const & ()  const 
//          { return *(const KNM<const_R>*) this;}
  void init(long nn,long mm) {
    ShapeOfArray::init(nn*mm);
    this->shapei.init(nn,1,nn);
    this->shapej.init(mm,nn,1),
    this->v=new R[nn*mm];}
    
  void resize(long nn,long mm) {     
    long kk=nn*mm;
    long n = this->shapei.n;
    long m = this->shapej.n;
    
    if( n !=nn && m != mm) 
     {
    KNM_ <R> old(*this); 
    long no=std::min(n,nn);
    long mo=std::min(m,mm);
    R *vo=this->v;
    
    // new mat 
    ShapeOfArray::init(kk);
    this->v=new R[this->n];
    this->shapei.init(nn,1,nn);
    this->shapej.init(mm,nn,1);
    
    if(this->v && vo)  // copy
        (*this)(SubArray(no),SubArray(mo)) = old(SubArray(no),SubArray(mo));
      
    delete []vo;
    }
        
   }
    
  void destroy(){delete [] this->v;this->n=0 ;}

};

//  Array with 3 indices
//  ---------------------
template<class R>
class KNMK: public KNMK_<R>{ public:

  KNMK(long n,long m,long k) 
     :KNMK_<R>(new R[n*m*k],n,m,k){}
  explicit KNMK(const KNMK_<R> & u)
     :KNMK_<R>(new R[u.size()],u.N(),u.M(),u.K()) 
     { KNMK_<R>::operator=(u);}
   KNMK(const KNMK<R> & u)
     :KNMK_<R>(new R[u.size()],u.N(),u.M(),u.K()) 
     { KNMK_<R>::operator=(u);}
     
  ~KNMK(){delete [] this->v;}
  
  KNMK& operator=(const KNMK_<const_R> & u)   
     { KNMK_<R>::operator=(u);return *this;}
  KNMK& operator=(const_R a)                  
     { KNMK_<R>::operator=(a);return *this;}
  KNMK& operator+=(const_R  a)                
     { KNMK_<R>::operator+=(a);return *this;}
  KNMK& operator-=(const_R  a)                
     { KNMK_<R>::operator-=(a);return *this;}
  KNMK& operator/=(const_R  a)                
     { KNMK_<R>::operator/=(a);return *this;}
  KNMK& operator*=(const_R  a)                
     { KNMK_<R>::operator*=(a);return *this;}
  KNMK& operator+=(const KNMK_<const_R> & u)  
     { KNMK_<R>::operator+=(u);return *this;}
  // ici jd 
  KNMK& operator-=(const KNMK_<const_R> & u)  
     { KNMK_<R>::operator-=(u);return *this;}
  KNMK& operator*=(const KNMK_<const_R> & u)   
     { KNMK_<R>::operator*=(u);return *this;}
  KNMK& operator/=(const KNMK_<const_R> & u)   
     { KNMK_<R>::operator/=(u);return *this;}
     
//  two opertor to cast to un array of constant          
//    operator KNMK_<const_R> & ()  
//       { return *  (KNMK_<const_R>*) this;}
//    operator KNMK_<const_R> const & ()  const 
//       { return *(const KNMK_<const_R>*) this;}  

//    operator KNMK<const_R> & ()  
//       { return *  (KNMK<const_R>*) this;}
//    operator KNMK<const_R> const & ()  const 
//       { return *(const KNMK<const_R>*) this;}  
};

//  -------------  optimization ---------------------
template<class R> 
class conj_KN_{public:
  const KN_<const_R> & a;
  conj_KN_(const KN_<const_R> & aa) : a(aa){}
};


inline const KN_<long> conj(const KN_<long> &a){ return a;}
inline const KN_<double> conj(const KN_<double> &a){ return a;}
inline const KN_<float> conj(const KN_<float> &a){ return a;}

//template<class R> conj_KN_<R> conj(const KN<R> &a){ return a;}
template<class R> conj_KN_<R> conj(const KN_<R> &a){ return a;}

template<class R> 
class DotStar_KN_{public: 
  const KN_<const_R>  a; const KN_<const_R>  b;
  DotStar_KN_(const KN_<const_R> & aa,const KN_<const_R> & bb) : a(aa),b(bb)  {}
 }; 

 
template<class R> 
class DotSlash_KN_{public: 
  const KN_<const_R>  a; const KN_<const_R>  b;
  DotSlash_KN_(const KN_<const_R> & aa,const KN_<const_R> & bb) : a(aa),b(bb)  {}
 }; 

template<class R> 
class Add_KN_{public: 
  const KN_<const_R>  a; const KN_<const_R>  b;
  Add_KN_(const KN_<const_R> & aa,const KN_<const_R> & bb) 
     : a(aa),b(bb)  { K_throwassert(SameShape(a,b));}
 };  
 
template<class R> 
class Sub_KN_{public: 
  const KN_<const_R>  a; const KN_<const_R>  b;
  Sub_KN_(const KN_<const_R> & aa,const KN_<const_R> & bb) 
    : a(aa),b(bb) { K_throwassert(SameShape(a,b));}
 };
 
template<class R> 
class Mulc_KN_ { public: 
  const KN_<const_R>  a;  const_R  b;
  Mulc_KN_(const KN_<const_R> & aa,const_R  bb) : a(aa),b(bb) {}
  Mulc_KN_(const Mulc_KN_<R> & aa,const_R  bb) : a(aa.a),b(aa.b*bb) {}
  Mulc_KN_ operator-() const {return Mulc_KN_(a,-b);}
  outProduct_KN_<R> operator*(const  TKN_<double> & bb)
{  return outProduct_KN_<R>(a,bb,b);} 

 };  

template<class R> 
class Add_Mulc_KN_ { public:
  const KN_<const_R> a,b;
  const R ca,cb; 
  Add_Mulc_KN_(const Mulc_KN_<R> & aa,const Mulc_KN_<R> & bb)  
        : a(aa.a),b(bb.a),ca(aa.b),cb(bb.b) { K_throwassert(SameShape(a,b));}
  Add_Mulc_KN_(const Mulc_KN_<R> & aa,const KN_<const_R> & bb,const R cbb) 
        : a(aa.a),b(bb),ca(aa.b),cb(cbb)  { K_throwassert(SameShape(a,b));}
  Add_Mulc_KN_(const KN_<const_R> & aa,const R caa,const KN_<const_R> & bb,const R cbb) 
        : a(aa),b(bb),ca(caa),cb(cbb) { K_throwassert(SameShape(a,b));}
 };  

template<class R> 
class if_arth_KN_ { public:
  const KN_<const_R> a,b,c;
  if_arth_KN_(const KN_<R> & aa,const KN_<R> & bb,const KN_<R> & cc)  
        : a(aa),b(bb),c(cc){ K_throwassert(SameShape(a,b)&&SameShape(a,c));}
 };  



template<class R> 
class Mul_KNM_KN_ { public:
  const KNM_<const_R> &A;
  const KN_<const_R> &b;
  Mul_KNM_KN_(const  KNM_<const_R>  &aa,const KN_<const_R>  &bb)  
        : A(aa),b(bb) {K_throwassert(SameShape(A.shapej,b));} 
};


ostream & operator<<(ostream & f,const ShapeOfArray & s);

template<class R> ostream & operator<<(ostream & f,const KN_<const_R>   & v);
template<class R> ostream & operator<<(ostream & f,const KNM_<const_R>  & v);
template<class R> ostream & operator<<(ostream & f,const KNMK_<const_R> & v);
template<class R> inline ostream & operator<<(ostream & f,const KN<const_R>   & v) 
    { return f << (const KN_<const_R> &) v;}
template<class R> inline ostream & operator<<(ostream & f,const KNM<const_R>  & v) 
    { return f << (const KNM_<const_R> &) v;}
template<class R> inline ostream & operator<<(ostream & f,const KNMK<const_R> & v) 
    { return f << (const KNMK_<const_R> &) v;}


template<class R> inline Add_KN_<R> operator+(const KN_<const_R> &a,const KN_<const_R> &b) 
    { return Add_KN_<R>(a,b);}
template<class R> inline Sub_KN_<R> operator-(const KN_<const_R> &a,const KN_<const_R> &b) 
    { return Sub_KN_<R>(a,b);}
template<class R> inline Mulc_KN_<R> operator*(const KN_<const_R> &a,const R &b) 
    { return Mulc_KN_<R>(a,b);}
template<class R> inline Mulc_KN_<R> operator/(const KN_<const_R> &a,const R &b) 
    { return Mulc_KN_<R>(a,1/b);}
template<class R> inline Mulc_KN_<R> operator*(const R &b,const KN_<const_R> &a) 
    { return Mulc_KN_<R>(a,b);}
template<class R> inline Mulc_KN_<R> operator-(const KN_<const_R> &a) 
    { return Mulc_KN_<R>(a,-1);}



template<class R> inline Add_Mulc_KN_<R> operator+(const  Mulc_KN_<R>& a,const Mulc_KN_<R> &b) 
    { return Add_Mulc_KN_<R>(a,b);}
template<class R> inline Add_Mulc_KN_<R> operator-(const  Mulc_KN_<R>& a,const Mulc_KN_<R> &b) 
    { return Add_Mulc_KN_<R>(a,b.a,-b.b);}

template<class R> inline Add_Mulc_KN_<R> operator+(const  Mulc_KN_<R>& a,const KN_<const_R> &b) 
   { return Add_Mulc_KN_<R>(a,b,R(1));}
template<class R> inline Add_Mulc_KN_<R> operator-(const  Mulc_KN_<R>& a,const KN_<const_R> &b) 
   { return Add_Mulc_KN_<R>(a,b,R(-1));}

template<class R> inline Add_Mulc_KN_<R> operator+(const KN_<const_R> & b,const  Mulc_KN_<R>& a) 
   { return Add_Mulc_KN_<R>(a,b,R(1));}

// modif FH mars 2007 
template<class R> inline Add_Mulc_KN_<R> operator-(const KN_<const_R> & a,const  Mulc_KN_<R>& b) 
   { return Add_Mulc_KN_<R>(a,R(1),b.a,-b.b);}// modif FH mars 2007  

template<class R> inline Mul_KNM_KN_<R> operator*(const  KNM_<const_R> & A,const  KN_<const_R> & b) 
    { return Mul_KNM_KN_<R>(A,b);}


template<class R> inline bool  SameShape(const ShapeOfArray & a,const Add_Mulc_KN_<R> & b) 
           { return SameShape(a,b.a) ;} 
template<class R> inline bool  SameShape(const ShapeOfArray & a,const if_arth_KN_<R> & b) 
           { return SameShape(a,b.a) ;} 
template<class R> inline bool  SameShape(const ShapeOfArray & a,const Add_KN_<R> & b) 
           { return SameShape(a,b.a) ;} 
template<class R> inline bool  SameShape(const ShapeOfArray & a,const Sub_KN_<R> & b) 
           { return SameShape(a,b.a) ;} 
template<class R> inline bool  SameShape(const ShapeOfArray & a,const Mulc_KN_<R> & b) 
           { return SameShape(a,b.a) ;} 
template<class R> inline bool  SameShape(const ShapeOfArray & a,const DotStar_KN_<R> & b) 
           { return SameShape(a,b.a) ;} 
template<class R> inline bool  SameShape(const ShapeOfArray & a,const DotSlash_KN_<R> & b) 
           { return SameShape(a,b.a) ;} 
template<class R> inline bool  SameShape(const ShapeOfArray & a,const Mul_KNM_KN_<R> & b) 
           { return a.n==b.A.N() ;} 
 inline bool  SameShape(const ShapeOfArray & a,const VirtualMatrice<double>::plusAx & b) 
           { return true ;} //  pas de test car la matrice peut etre rectangulaire
 inline bool  SameShape(const ShapeOfArray & a,const VirtualMatrice<double>::plusAtx & b) 
           { return true ;} //  pas de test car la matrice peut etre rectangulaire
 inline bool  SameShape(const ShapeOfArray & a,const VirtualMatrice<complex<double> >::plusAx & b) 
           { return true ;} //  pas de test car la matrice peut etre rectangulaire
 inline bool  SameShape(const ShapeOfArray & a,const VirtualMatrice<complex<double> >::plusAtx & b) 
           { return true ;} //  pas de test car la matrice peut etre rectangulaire

 inline bool  SameShape(const ShapeOfArray & a,const double) 
           { return true;} 
 inline bool  SameShape(const ShapeOfArray & a,const complex<double>) 
           { return true;} 
 inline bool  SameShape(const ShapeOfArray & a,const complex<float>) 
           { return true;}            

template<class R>
 inline bool SameShape(KNM<R>& m, const outProduct_KN_<R>& p)
 { return p.a.N()>=m.N() && m.M()>=p.b.N(); } 

template<class R> inline long SameAdress(const KN_<R> &a, const KN_<R> &b) { return &a[0]==&b[0];}
// bof -bof 
//template<class R> inline
//  KN_<R>::operator KN<R> &() { return *(KN<R> *) (void *) this;}
//template<class R> inline
//  KN_<R>::operator const KN<R> &() const { return *(const KN<R> *) ( const void *) this;}

//  operateur y=Ax-b ou y=Ax + b pour le GC
template<class R> 
   PplusQ< typename VirtualMatrice<R>::plusAx, Mulc_KN_<R> >  operator-(const typename VirtualMatrice<R>::plusAx & A,const KN_<R> & B)  
    { return PplusQ< typename VirtualMatrice<R>::plusAx, Mulc_KN_<R> >(A,Mulc_KN_<R>(B,R(-1.)));}

template<class R> 
   PplusQ< typename VirtualMatrice<R>::plusAx, KN_<R> >  operator+(const typename VirtualMatrice<R>::plusAx & A,const KN_<R> & B)  
    { return PplusQ< typename VirtualMatrice<R>::plusAx, KN_<R> >(A,B);}

template<class R> 
   PplusQ< typename VirtualMatrice<R>::plusAx, Mulc_KN_<R> >  operator-(const typename VirtualMatrice<R>::plusAx & A,const KN<R> & B)  
    { return PplusQ< typename VirtualMatrice<R>::plusAx, Mulc_KN_<R> >(A,Mulc_KN_<R>(B,R(-1.)));}

template<class R> 
   PplusQ< typename VirtualMatrice<R>::plusAx, KN_<R> >  operator+(const typename VirtualMatrice<R>::plusAx & A,const KN<R> & B)  
    { return PplusQ< typename VirtualMatrice<R>::plusAx, KN_<R> >(A,B);}
    

template<class R>
KN_<R> diagonal(const KNM<R> & A) { 
  K_throwassert(A.N() == A.M()); 
  return KN_<R>(A,SubArray(A.N(),0,A.N()+1));}

// to def  inv permutation FH mars 2006 
class Inv_KN_long{ public:
  KN_<long>  t;
  Inv_KN_long(const KN_<long> & v)
   : t(v) {}
  Inv_KN_long( KN_<long> const  * & v)
   : t(*v) {}
  operator const KN_<long> & () const {return t;}
};


template<class R,typename A,typename B=R> class  F_KN_ 
{ 
  public: 
  A (*f)(B);
  KN_<R> a;
  long N() const {return a.N();}
  F_KN_( A (*ff)(B),const KN_<R> & aa): f(ff),a(aa) {}
  A operator[](long i) const { return f(a[i]);}
  bool check(long n)  const { return  n <= a.N() || a.constant(); }
  bool constant() const {return a.constant();}
}; 

template<class R,typename A,typename B>
inline bool  SameShape(const ShapeOfArray & a,const F_KN_<R,A,B>  & b) 
           { return  !a.step || b.constant()  || a.n == b.N() ;} 
           
#include "RNM_tpl.hpp"
#ifdef K_throwassert
#undef K_throwassert
#endif
#endif