rnm.hpp

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

  KN_(const KN_<R> & U,const SubArray & sa)  : ShapeOfArray(U,sa),v(U.v + U.index(sa.start)) {}

  KN_ operator()(const SubArray & sa) const { return KN_(*this,sa);} // sub array 
  
  R & operator[](long i) const {return v[index(i)];}
  R & operator()(long i) const {return v[index(i)];}
  R & operator[](int i) const {return v[index(i)];}
  R & operator()(int i) const {return v[index(i)];}

  R operator,(const KN_<const_R> & v) const; // dot  product 

   KN_& operator  =(const KN_<const_R> & u)  ;
   KN_& operator +=(const KN_<const_R> & u)  ;
   KN_& operator -=(const KN_<const_R> & u)  ;
  
   KN_& operator *=(const KN_<const_R> & u)  ;
   KN_& operator /=(const KN_<const_R> & u)  ;
  
  
   KN_& operator = (const_R  a) ;  
   KN_& operator +=(const_R  a) ;
   KN_& operator -=(const_R  a) ;  
   KN_& operator /=(const_R  a) ;
   KN_& operator *=(const_R  a) ;
  
   KN_& operator  = (R*  a) { return operator =(KN_<R>(a,n));}
   KN_& operator += (R*  a) { return operator+=(KN_<R>(a,n));}  
   KN_& operator -= (R*  a) { return operator-=(KN_<R>(a,n));}  
   KN_& operator *= (R*  a) { return operator*=(KN_<R>(a,n));}  
   KN_& operator /= (R*  a) { return operator/=(KN_<R>(a,n));}  


  
  R min() const ;
  R max() const ;
  R sum() const ;
  double norm() const ;
  double l2() const ;
  double l1() const ;
  double linfty() const ;
  double lp(double p) const ;
  
  template<class T> long last(const T &) const;
  template<class T> long first(const T &) const;
  
  void map(R (*f)(R )); // apply the f fonction a all element of the array
  void map(R (*f)(const  R& )); // apply the f fonction a all element of the array

 template<class T>
   void set(R (*f)(const  T& ),KN_<T> & u); // apply the f fonction a all element of the array u
  
   KN_& operator =(const DotStar_KN_<R> & u) ;
   KN_& operator+=(const DotStar_KN_<R> & u) ;
   KN_& operator-=(const DotStar_KN_<R> & u) ;
   KN_& operator*=(const DotStar_KN_<R> & u) ;
   KN_& operator/=(const DotStar_KN_<R> & u) ;

   KN_& operator =(const DotSlash_KN_<R> & u) ;
   KN_& operator+=(const DotSlash_KN_<R> & u) ;
   KN_& operator-=(const DotSlash_KN_<R> & u) ;
   KN_& operator*=(const DotSlash_KN_<R> & u) ;
   KN_& operator/=(const DotSlash_KN_<R> & u) ;

   KN_& operator =(const if_KN_<R> & u) ;
   KN_& operator+=(const if_KN_<R> & u) ;
   KN_& operator-=(const if_KN_<R> & u) ;
   KN_& operator*=(const if_KN_<R> & u) ;
   KN_& operator/=(const if_KN_<R> & u) ;

   KN_& operator =(const ifnot_KN_<R> & u) ;
   KN_& operator+=(const ifnot_KN_<R> & u) ;
   KN_& operator-=(const ifnot_KN_<R> & u) ;
   KN_& operator*=(const ifnot_KN_<R> & u) ;
   KN_& operator/=(const ifnot_KN_<R> & u) ;

   KN_& operator =(const Add_KN_<R> & u) ;
   KN_& operator+=(const Add_KN_<R> & u) ;
   KN_& operator-=(const Add_KN_<R> & u) ;
   KN_& operator*=(const Add_KN_<R> & u) ;
   KN_& operator/=(const Add_KN_<R> & u) ;
  
   template<class I,class T> KN_& operator =  (const KN_ITAB<T,I> & u);
   template<class I,class T> KN_& operator +=  (const KN_ITAB<T,I> & u);
   template<class I,class T> KN_& operator -=  (const KN_ITAB<T,I> & u);
   template<class I,class T> KN_& operator *=  (const KN_ITAB<T,I> & u);
   template<class I,class T> KN_& operator /=  (const KN_ITAB<T,I> & u);


   KN_ITAB< KN_<R>,const KN_<int> >  operator()(const KN_<int> &itab) ;
   KN_ITAB< KN_<R>,const KN_<long> >  operator()(const KN_<long> &itab) ;
   KN_ITAB<const KN_<R>,const KN_<int> > operator()(const KN_<int> &itab) const ;
   KN_ITAB<const KN_<R>,const KN_<long> >  operator()(const KN_<long> &itab) const ;


    

  
   KN_& operator =(const Sub_KN_<R> & u) ;
   KN_& operator-=(const Sub_KN_<R> & u) ;
   KN_& operator+=(const Sub_KN_<R> & u) ;
   KN_& operator*=(const Sub_KN_<R> & u) ;
   KN_& operator/=(const Sub_KN_<R> & u) ;
  
   KN_& operator =(const Mulc_KN_<R> & u) ;
   KN_& operator+=(const Mulc_KN_<R> & u) ;
   KN_& operator-=(const Mulc_KN_<R> & u) ;
   KN_& operator*=(const Mulc_KN_<R> & u) ;
   KN_& operator/=(const Mulc_KN_<R> & u) ;
  
   KN_& operator =(const Add_Mulc_KN_<R> & u) ;
   KN_& operator+=(const Add_Mulc_KN_<R> & u) ;
   KN_& operator-=(const Add_Mulc_KN_<R> & u) ;
   KN_& operator*=(const Add_Mulc_KN_<R> & u) ;
   KN_& operator/=(const Add_Mulc_KN_<R> & u) ;

   KN_& operator =(const if_arth_KN_<R> & u) ;
   KN_& operator+=(const if_arth_KN_<R> & u) ;
   KN_& operator-=(const if_arth_KN_<R> & u) ;
   KN_& operator*=(const if_arth_KN_<R> & u) ;
   KN_& operator/=(const if_arth_KN_<R> & u) ;

  
   KN_& operator =(const Mul_KNM_KN_<R> & u) ; 
   KN_& operator+=(const Mul_KNM_KN_<R> & u) ; 
   KN_& operator-=(const Mul_KNM_KN_<R> & u) ; 
   KN_& operator*=(const Mul_KNM_KN_<R> & u) ; 
   KN_& operator/=(const Mul_KNM_KN_<R> & u) ; 
  
 //  KN_& operator =(const MatriceCreuseMulKN_<R> & ) ;
 //  KN_& operator +=(const MatriceCreuseMulKN_<R> & ) ;
   KN_& operator =(const typename VirtualMatrice<R>::plusAx & Ax)  
    {*this=R(); Ax.A->addMatMul(Ax.x,*this);return *this;}
   KN_& operator =(const typename VirtualMatrice<R>::plusAtx & Ax)  
    {*this=R(); Ax.A->addMatTransMul(Ax.x,*this);return *this;}
   KN_& operator +=(const typename VirtualMatrice<R>::plusAx & Ax)  
    {  Ax.A->addMatMul(Ax.x,*this);return *this;}
   KN_& operator +=(const typename VirtualMatrice<R>::plusAtx & Ax)  
    {  Ax.A->addMatTransMul(Ax.x,*this);return *this;}
   KN_& operator =(const typename VirtualMatrice<R>::solveAxeqb & Ab)  
    {*this=R(); Ab.A->Solve(*this,Ab.b);return *this;}
    
  template<class  A,class B,class C> KN_&  operator =  (const F_KN_<A,B,C>  & u) ;
  template<class  A,class B,class C> KN_&  operator +=  (const F_KN_<A,B,C>  & u) ;
  template<class  A,class B,class C> KN_&  operator -=  (const F_KN_<A,B,C>  & u) ;
  template<class  A,class B,class C> KN_&  operator /=  (const F_KN_<A,B,C>  & u) ;
  template<class  A,class B,class C> KN_&  operator *=  (const F_KN_<A,B,C>  & u) ;
    
   
//   KN_& operator =(const MatriceCreuseDivKN_<R> &)  ;

 friend   ostream & operator<< <R>(ostream & f,const KN_<const_R> & v)  ;
 
  KN_(R *u,const ShapeOfArray & s):ShapeOfArray(s),v(u){}
  KN_(R *u,long nn,long s):ShapeOfArray(nn,s),v(u){}
  KN_(R *u,long nn,long s,long nextt):ShapeOfArray(nn,s,nextt),v(u){}
  KN_(R *u,long nn):ShapeOfArray(nn),v(u){}


  TKN_<R>  t() ; // transpose
  const TKN_<R>  t() const ; // transpose
  notKN_<R>  operator!()  ; //  not
  const  notKN_<R>  operator!() const ; // not

  //  operator KN<R> &();
  // operator const KN<R> &() const;

 private:
 
  KN_&  operator++(){K_throwassert(next>=0);v += next;return *this;} //    ++U
  KN_&  operator--(){K_throwassert(next>=0);v -= next;return *this;} //    --U
  KN_   operator++(int ){K_throwassert(next>=0); KN_ old=*this;v = v +next;return old;} // U++
  KN_   operator--(int ){K_throwassert(next>=0); KN_ old=*this;v = v -next;return old;} // U++
 
  KN_(const KN_<R> & u,long offset) :ShapeOfArray(u),v(&u[offset]){} 
  KN_(const KN_<R> & u,const ShapeOfArray &sh,long startv=0)
         :ShapeOfArray(sh*u.step),v(&u[startv]){}
  KN_(const KN_<R> & u,long nnext,const ShapeOfArray &sh,long startv=0)
         :ShapeOfArray(sh.n,sh.step*u.step,nnext),v(&u[startv]){ }

  
// friend class KN_<R>;   
 friend class KNM_<R>;   
 friend class KNMK_<R>;   
 friend class KN<R>;   
 friend class KNM<R>;   
 friend class KNMK<R>;   
  
};

template<class R>
class KNM_: public KN_<R> {
  public:
  ShapeOfArray shapei;
  ShapeOfArray shapej;
  public:
  long IsVector1() const  {  return (shapei.n*shapej.n) ==  this->n ;} 
  long N() const {return shapei.n;}
  long M() const {return shapej.n;}  
  long size() const { return shapei.n*shapej.n;}
  
  KNM_(R* u,const ShapeOfArray & s,
            const ShapeOfArray & si,
            const ShapeOfArray & sj)
             : KN_<R>(u,s),shapei(si),shapej(sj){} 
  KNM_(R* u,long n,long m)
             : KN_<R>(u,ShapeOfArray(n*m)),shapei(n,1,n),shapej(m,n,1){}
  KNM_(R* u,long n,long m,long s)
             : KN_<R>(u,ShapeOfArray(n*m,s)),shapei(n,1,n),shapej(m,n,1){}                     
  KNM_(KN_<R> u,long n,long m) 
             : KN_<R>(u,ShapeOfArray(m*n)),shapei(n,1,n),shapej(m,n,1){ }
             
  KNM_(const KN_<R> &u,const ShapeOfArray & si,const ShapeOfArray & sj,long offset=0) 
             : KN_<R>(&u[offset],si.last()+sj.last()+1,u.step),shapei(si),shapej(sj) 
             {K_throwassert( offset>=0 && this->n+ (this->v-(R*)u) <= u.n);}
  KNM_(const KN_<R> &u,const ShapeOfArray & si,const ShapeOfArray & sj,long offset,long nnext) 
             : KN_<R>(&u[offset],si.last()+sj.last()+1,u.step,nnext),shapei(si),shapej(sj) 
             {K_throwassert( offset>=0 && this->n+ (this->v-(R*)u) <= u.n);}
  
  KNM_(KNM_<R> U,const SubArray & si,const SubArray & sj)  
             :KN_<R>(U,SubArray(U.ij(si.len1(),sj.len1())+1,U.ij(si.start,sj.start))),
                     shapei(U.shapei,si),shapej(U.shapej,sj){} 

  KNM_(KNM_<R> U,const SubArray & sa,const SubArray & si,const SubArray & sj)  
             :KN_<R>(U,SubArray(sa)),shapei(U.shapei,si),shapej(U.shapej,sj){} 

  KNM_(const KNM_<R> & u) 
             :KN_<R>(u),shapei(u.shapei),shapej(u.shapej) {}

  KNM_ operator()(const SubArray & sa,const SubArray & sb) const 
            { return KNM_(*this,sa,sb);} // sub array 
  
  long ij(long i,long j) const   
            { return shapei.index(i)+shapej.index(j);}
  long indexij(long i,long j)        const   
            { return this->index(shapei.index(i)+shapej.index(j));}
  R & operator()(long i,long j)     const   
            { return this->v[indexij(i,j)];}
  R & operator()(int i,int j)     const   
            { return this->v[indexij(i,j)];}
            
            
  KN_<R> operator()(const char,long j    )  const   // une colonne j  ('.',j)
            { return KN_<R>(&this->v[this->index(shapej.index(j))],shapei*this->step);} 
  KN_<R> operator()(long i    ,const char)  const   // une ligne i  (i,'.')
            { return KN_<R>(&this->v[this->index(shapei.index(i))],shapej*this->step);}  
  KN_<R> operator()(const char,int j    )  const   // une colonne j  ('.',j)
            { return KN_<R>(&this->v[this->index(shapej.index(j))],shapei*this->step);} 
  KN_<R> operator()(int i    ,const char)  const   // une ligne i  (i,'.')
            { return KN_<R>(&this->v[this->index(shapei.index(i))],shapej*this->step);}  
  KN_<R> operator()(const char,const char)  const   // tous 
            { return *this;}
  KNM_<R> t() const
    { return KNM_<R>(this->v,*this,shapej,shapei);} //  before  { return KNM_<R>(*this,shapej,shapei,v);}
  
   KNM_& operator =(const KNM_<const_R> & u) ;
   KNM_& operator =(const_R a)               ;
   KNM_& operator+=(const_R a)               ;
   KNM_& operator-=(const_R a)               ; 
   KNM_& operator/=(const_R a)               ;
   KNM_& operator*=(const_R a)               ; 
   KNM_& operator+=(const KNM_<const_R> & u) ;
   KNM_& operator-=(const KNM_<const_R> & u) ;
   KNM_& operator*=(const KNM_<const_R> & u) ;
   KNM_& operator/=(const KNM_<const_R> & u) ;
   
   KNM_ &operator =(const outProduct_KN_<R> &);
   KNM_ &operator +=(const outProduct_KN_<R> &);
   KNM_ &operator -=(const outProduct_KN_<R> &);
   KNM_ &operator /=(const outProduct_KN_<R> &); // bofbof
   KNM_ &operator *=(const outProduct_KN_<R> &); // bofbof

private:  
  KNM_& operator++() {this->v += this->next;return *this;} // ++U
  KNM_& operator--() {this->v -= this->next;return *this;} // ++U
  KNM_  operator++(int ){KNM_<R> old=*this;this->v = this->v +this->next;return old;} // U++
  KNM_  operator--(int ){KNM_<R> old=*this;this->v = this->v -this->next;return old;} // U--


 friend class KN_<R>;   
// friend class KNM_<R>;   
 friend class KNMK_<R>;   
 friend class KN<R>;   
 friend class KNM<R>;   
 friend class KNMK<R>;   
};

template<class T,class I> 
struct KN_ITAB
{
  KN_ITAB(const T &vv,const I &iindex) : v(vv),index(iindex) {}
  T  v;
  I  index;
  KN_ITAB & operator=(const T & t);
  KN_ITAB & operator+=(const T & t);
  KN_ITAB & operator-=(const T & t);
  KN_ITAB & operator*=(const T & t);
  KN_ITAB & operator/=(const T & t);
  typename T::R & operator[](long i){ return v[index[i]];}
  const typename T::R &  operator[](long i) const { return v[index[i]];}
  long N() const { return index.N();}    
};

template<class R>  KN_ITAB<const KN_<R>,const KN_<int> >  KN_<R>::operator()(const KN_<int>  &itab) const { return KN_ITAB<const KN_<R>,const KN_<int> > (*this,itab);}
template<class R>  KN_ITAB<const KN_<R>,const KN_<long> >  KN_<R>::operator()(const KN_<long>  &itab) const { return KN_ITAB<const KN_<R>,const KN_<long> > (*this,itab);}
template<class R>  KN_ITAB< KN_<R>,const KN_<int> >  KN_<R>::operator()(const KN_<int>  &itab) { return KN_ITAB<KN_<R>,const KN_<int> > (*this,itab);}
template<class R>  KN_ITAB< KN_<R>,const KN_<long> >  KN_<R>::operator()(const KN_<long>  &itab) { return KN_ITAB<KN_<R>,const KN_<long> > (*this,itab);}


template<class R>
struct TKN_:public KN_<R> {
    TKN_(const KN_<R> &x) : KN_<R>(x) {}
};

template<class R>
struct notKN_:public KN_<R> {
    notKN_(const KN_<R> &x) : KN_<R>(x) {}
    notnotKN_<R>  operator!()  ; //  not
    const  notnotKN_<R>  operator!() const ; // not
};

template<class R>
struct notnotKN_:public KN_<R> {
    notnotKN_(const notKN_<R> &x) : KN_<R>(x) {}
    notKN_<R>  operator!()  ; //  notnot
    const  notKN_<R>  operator!() const ; // notnot
};

template<class R>
TKN_<R>  KN_<R>::t() { return *this;} // transpose

template<class R>
const TKN_<R>  KN_<R>::t() const { return *this;} // transpose

template<class R>
notKN_<R>  KN_<R>::operator!() { return *this;} // not

template<class R>
const notKN_<R>  KN_<R>::operator!() const { return *this;} // not

template<class R>
notnotKN_<R>  notKN_<R>::operator!() { return *this;} // not

template<class R>
const notnotKN_<R>  notKN_<R>::operator!() const { return *this;} // not


template<class R>
struct outProduct_KN_ {
    const KN_<R>  a,b;
    R c;
    outProduct_KN_(const KN_<R> & aa, const KN_<R> &bb,R cc=(R)1) : a(aa),b(bb),c(cc) {}
    outProduct_KN_(const KN_<R> * aa, const KN_<R> &bb,R cc=(R)1) : a(*aa),b(bb),c(cc) {}
    outProduct_KN_(const KN_<R> * aa, const KN_<R> *bb,R cc=(R)1) : a(*aa),b(*bb),c(cc) {}
    outProduct_KN_(const Mulc_KN_<R> & aa,const KN_<R> & bb) : a(aa.a),b(bb),c(aa.b) {}    
    outProduct_KN_ operator * (R cc) { return outProduct_KN_(a,b,c*cc);}    
};

template<class R>
struct if_KN_ {
    const KN_<R> & a,&b;
    R c;
    if_KN_(const KN_<R> & aa, const KN_<R> &bb,R cc=1.) : a(aa),b(bb),c(cc) {}
    if_KN_ operator * (R cc) { return if_KN_(a,b,c*cc);}    
};

template<class R>