rnm.hpp

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

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


template<class R> 
outProduct_KN_<R> operator*(const KN_<R> &a,const TKN_<R> &b) 
{ return outProduct_KN_<R>(a,b);}

template<class R> 
ifnot_KN_<R> operator*(const KN_<R> &a,const notKN_<R> &b) 
{ return ifnot_KN_<R>(b,a);}

template<class R> 
ifnot_KN_<R> operator*(const KN_<R> &a,const notnotKN_<R> &b) 
{ return if_KN_<R>(b,a);}

template<class R> 
ifnot_KN_<R> operator*(const notKN_<R> &b,const KN_<R> &a) 
{ return ifnot_KN_<R>(b,a);}

template<class R> 
ifnot_KN_<R> operator*(const notnotKN_<R> &b,const KN_<R> & a) 
{ return if_KN_<R>(b,a);}


template<class R> 
R operator*(const TKN_<R> &a,const KN_<R> &b) 
{ return (a,b);}

template<class R>
class KNMK_: public KN_<R> {
  friend class KNMK<R>;
  public:
  ShapeOfArray shapei;
  ShapeOfArray shapej;
  ShapeOfArray shapek;
  public:
  long IsVector1() const {  return (shapei.n*shapej.n*shapek.n) == this->n ;} 
  long N() const {return shapei.n;}
  long M() const {return shapej.n;}
  long K() const {return shapek.n;}
  long size() const { return shapei.n*shapej.n*shapek.n;}
  KNMK_(const ShapeOfArray & s,
        const ShapeOfArray & si,
        const ShapeOfArray & sj,
        const ShapeOfArray & sk,
	    R * u)
    : KN_<R>(u,s),shapei(si),shapej(sj),shapek(sk){} 
    
  KNMK_(R* u,long n,long m,long k)
    : KN_<R>(u, ShapeOfArray(n*m*k)),shapei(n,1,n),shapej(m,n,1),shapek(k,n*m,n*m){};
    
//  KNMK_(const KN_<R> & u,long n,long m,long k)
//   : KN_<R>(ShapeOfArray(n*m*k)),shapei(n,1,n),shapekj(m,n,1),u),
//     shapek(k,n*m,n*m){};

  KNMK_(const KNMK_<R> &U,const SubArray & si,const SubArray & sj,const SubArray & sk)  :
    KN_<R>(U,SubArray(U.ijk(si.len1(),sj.len1(),sk.len1())+1,
                       U.ijk(si.start,sj.start,sk.start))),
                       shapei(U.shapei,si),
                       shapej(U.shapej,sj),
                       shapek(U.shapek,sk){} 

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

    
  long ijk(long i,long j,long k) const 
              { return shapei.index(i)+shapej.index(j)+shapek.index(k);}
  long indexijk(long i,long j,long k) const 
              {return this->index(shapei.index(i)+shapej.index(j)+shapek.index(k));} 
                           
  R & operator()(long i,long j,long k)   const   {return this->v[indexijk(i,j,k)];}
  R & operator()(int i,int j,int k)   const   {return this->v[indexijk(i,j,k)];}
  
//  pas de tableau suivant
 KN_<R>  operator()(const char ,long j,long k)  const  { // le tableau (.,j,k) 
        return KN_<R>(*this,-1,shapei,shapej[j]+shapek[k]);}
 KN_<R>  operator()(long i,const char ,long k)  const  { // le tableau (i,.,k) 
        return KN_<R>(*this,-1,shapej,shapei[i]+shapek[k]);}
 KN_<R>  operator()(long i,long j,const char )  const  { // le tableau (i,j,.) 
        return KN_<R>(*this,-1,shapek,shapei[i]+shapej[j]);}

 KN_<R>  operator()(const char ,int j,int k)  const  { // le tableau (.,j,k) 
        return KN_<R>(*this,-1,shapei,shapej[j]+shapek[k]);}
 KN_<R>  operator()(int i,const char ,int k)  const  { // le tableau (i,.,k) 
        return KN_<R>(*this,-1,shapej,shapei[i]+shapek[k]);}
 KN_<R>  operator()(int i,int j,const char )  const  { // le tableau (i,j,.) 
        return KN_<R>(*this,-1,shapek,shapei[i]+shapej[j]);}
//                                              
 KNM_<R>  operator()(const char ,const char ,long k)  const  { // le tableau (.,.,k) 
        return KNM_<R>(*this,shapei,shapej,shapek[k],shapek.next);} // step = n*m
 //attention les suivants ne marche pas
 KNM_<R>  operator()(const char ,long j,const char )  const  { // le tableau (.,j,.) 
        return KNM_<R>(*this,shapei,shapek,shapej[j],-1/*shapej.next*/);} // step = n
        
 KNM_<R>  operator()(long i,const char ,const char )  const  { // le tableau (i,.,.) 
        return KNM_<R>(*this,shapej,shapek,shapei[i],-1/*shapei.next*/);}  // step = 1

 KNM_<R>  operator()(const char ,const char ,int k)  const  { // le tableau (.,.,k) 
        return KNM_<R>(*this,shapei,shapej,shapek[k],shapek.next);} // step = n*m
 //attention les suivants ne marche pas
 KNM_<R>  operator()(const char ,int j,const char )  const  { // le tableau (.,j,.) 
        return KNM_<R>(*this,shapei,shapek,shapej[j],-1/*shapej.next*/);} // step = n
        
 KNM_<R>  operator()(int i,const char ,const char )  const  { // le tableau (i,.,.) 
        return KNM_<R>(*this,shapej,shapek,shapei[i],-1/*shapei.next*/);}  // step = 1

   KNMK_& operator =(const KNMK_<const_R> & u) ;
   KNMK_& operator+=(const KNMK_<const_R> & u)  ;
   KNMK_& operator-=(const KNMK_<const_R> & u)  ;
   KNMK_& operator/=(const KNMK_<const_R> & u)  ;
   KNMK_& operator*=(const KNMK_<const_R> & u)  ;
   KNMK_& operator =(const_R a)  ; 
   KNMK_& operator+=(const_R a)  ;
   KNMK_& operator-=(const_R a)  ;
   KNMK_& operator/=(const_R a)  ;
   KNMK_& operator*=(const_R a)  ;

  KNMK_  operator()(SubArray si,SubArray sj,SubArray sk) const 
        {return KNMK_(*this,si,sj,sk);}

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

};



template<class R>
class KN :public KN_<R> { public:

  typedef R K;

 // explicit  KN(const R & u):KN_<R>(new R(uu),1,0) {}
  KN() : KN_<R>(0,0) {}
  KN(long nn) : KN_<R>(new R[nn],nn)         {} 
  KN(long nn, R * p) : KN_<R>(new R[nn],nn)  
    { KN_<R>::operator=(KN_<R>(p,nn));}
  KN(long nn,R (*f)(long i) ) : KN_<R>(new R[nn],nn) 
        {for(long i=0;i<this->n;i++) this->v[i]=f(i);}  
  KN(long nn,const  R & a) : KN_<R>(new R[nn],nn) 
        { KN_<R>::operator=(a);} 
  KN(long nn,long s,const  R  a) : KN_<R>(new R[nn],nn,s) 
        { KN_<R>::operator=(a);} 
  template<class S>   KN(const KN_<S> & s):KN_<R>(new R[s.n],s.n) 
        {for (long i=0;i<this->n;i++) this->v[i] = s[i];}
  template<class S>  KN(const KN_<S> & s,R (*f)(S )):KN_<R>(new R[s.n],s.n) 
        {for (long i=0;i<this->n;i++) this->v[i] = f(s[i]);}
  KN(const KN<R> & u):KN_<R>(new R[u.n],u.n)
        { KN_<R>::operator=(u);}
  KN(bool ,KN<R> & u):KN_<R>(u) {u.v=0;u.n=0;}// remove copy for return of local KN. 
    
  //  explicit KN(const KN_<R> & u):KN_<R>(new R[u.n],u.n) 
  //      { KN_<R>::operator=(u);}
        
  ~KN(){delete [] this->v;}
  void CheckSet() { if(!(this->n)) {cerr << "Error RNM set array\n";K_throwassert(0); exit(1);}}
   KN& operator  = (R*  a) { CheckSet(); return operator =(KN_<R>(a,this->n));}
   KN& operator += (R*  a) { CheckSet(); return operator+=(KN_<R>(a,this->n));}  
   KN& operator -= (R*  a) { CheckSet(); return operator-=(KN_<R>(a,this->n));}  
   KN& operator *= (R*  a) { CheckSet(); return operator*=(KN_<R>(a,this->n));}  
   KN& operator /= (R*  a) { CheckSet(); return operator/=(KN_<R>(a,this->n));}  
  
   KN& operator =(const_R a)  
        { if(this->unset()) set(new R[1],1,0,0); KN_<R>::operator= (a);return *this;}
   KN& operator =(const KN_<R>& a)  
        { if(this->unset()) set(new R[a.N()],a.N()); KN_<R>::operator= (a);return *this;}                
   KN& operator =(const KN<R>& a)  
        { if(this->unset()) set(new R[a.N()],a.N()); KN_<R>::operator= (a);return *this;}                
   KN& operator =(const Add_KN_<R> & u)  
        { if(this->unset()) set(new R[u.a.N()],u.a.N());KN_<R>::operator=(u);return *this;}
   KN& operator =(const DotStar_KN_<R> & u)  
        { if(this->unset()) set(new R[u.a.N()],u.a.N());KN_<R>::operator=(u);return *this;}
   KN& operator =(const if_KN_<R> & u)  
        { if(this->unset()) set(new R[u.a.N()],u.a.N());KN_<R>::operator=(u);return *this;}
   KN& operator =(const ifnot_KN_<R> & u)  
        { if(this->unset()) set(new R[u.a.N()],u.a.N());KN_<R>::operator=(u);return *this;}
   KN& operator =(const DotSlash_KN_<R> & u)  
        { if(this->unset()) set(new R[u.a.N()],u.a.N());KN_<R>::operator=(u);return *this;}
   KN& operator =(const Sub_KN_<R> & u)  
        { if(this->unset()) set(new R[u.a.N()],u.a.N());KN_<R>::operator=(u);return *this;}
   KN& operator =(const Mulc_KN_<R> & u)  
        { if(this->unset()) set(new R[u.a.N()],u.a.N());KN_<R>::operator=(u);return *this;}
   KN& operator =(const Add_Mulc_KN_<R> & u)  
        { if(this->unset()) set(new R[u.a.N()],u.a.N());KN_<R>::operator=(u);return *this;}
   KN& operator =(const if_arth_KN_<R> & u)  
        { if(this->unset()) set(new R[u.a.N()],u.a.N());KN_<R>::operator=(u);return *this;}
        
        
   KN& operator =(const Mul_KNM_KN_<R> & u) 
        { if(this->unset()) set(new R[u.b.N()],u.b.N());KN_<R>::operator=(u);return *this;}
//   KN& operator =(const MatriceCreuseMulKN_<R> & Ax) 
//       {if(this->unset()) set(new R[Ax.v.N()],Ax.v.N()); KN_<R>::operator=(Ax);return *this;}
//   KN& operator +=(const MatriceCreuseMulKN_<R> & Ax) 
//       {if(this->unset()) set(new R[Ax.v.N()],Ax.v.N()); KN_<R>::operator+=(Ax);return *this;}
//   KN& operator =(const MatriceCreuseDivKN_<R> & A1x)  
//       { if(this->unset()) set(new R[A1x.v.N()],A1x.v.N());KN_<R>::operator=(A1x);return *this;}
  // correcton aout 2007 FH  add N,M flied in VirtualMatrice
   KN& operator =(const typename VirtualMatrice<R>::plusAx & Ax)  
        { if(this->unset() && Ax.A->N ) set(new R[Ax.A->N],Ax.A->N);KN_<R>::operator=(Ax);return *this;}
   KN& operator =(const typename VirtualMatrice<R>::solveAxeqb & Ab)  
        { if(this->unset()) set(new R[Ab.b.N()],Ab.b.N());KN_<R>::operator=(Ab);return *this;}
   KN& operator +=(const typename  VirtualMatrice<R>::plusAx & Ax)  
  { if(this->unset()  && Ax.A->N) {
        set(new R[Ax.A->N],Ax.A->N);
        KN_<R>::operator=(R());}
    KN_<R>::operator+=(Ax);
    return *this;}
   KN& operator =(const typename VirtualMatrice<R>::plusAtx & Ax)  
        { if(this->unset()&&Ax.A->M) set(new R[Ax.A->M],Ax.A->M);KN_<R>::operator=(Ax);return *this;}
   KN& operator +=(const typename VirtualMatrice<R>::plusAtx & Ax)  
  { if(this->unset()&&Ax.A->M) {
       set(new R[Ax.A->M],Ax.A->M);
      KN_<R>::operator=(R());}
      KN_<R>::operator+=(Ax);
     return *this;}
// end correcton FH
   template<class P,class Q> 
     KN& operator =(const  PplusQ<P,Q> & PQ)  
      { *this=PQ.p; *this+=PQ.q;return *this; } 
   template<class P,class Q> 
     KN& operator +=(const  PplusQ<P,Q> & PQ)  
      { *this+=PQ.p; *this+=PQ.q;return *this; } 
           
   KN& operator -=(const_R a)  
        { KN_<R>::operator-=(a);return *this;}
   KN& operator -=(const KN_<R>& a)  
        { KN_<R>::operator-= (a);return *this;}
   KN& operator -=(const Add_KN_<R> & u)  
        { KN_<R>::operator-=(u);return *this;}
   KN& operator -=(const DotStar_KN_<R> & u)  
        { KN_<R>::operator-=(u);return *this;}
   KN& operator -=(const DotSlash_KN_<R> & u)  
        { KN_<R>::operator-=(u);return *this;}
   KN& operator -=(const Sub_KN_<R> & u)  
        { KN_<R>::operator-=(u);return *this;}
   KN& operator -=(const Mulc_KN_<R> & u)  
        { KN_<R>::operator-=(u);return *this;}
   KN& operator -=(const Add_Mulc_KN_<R> & u)  
        { KN_<R>::operator-=(u);return *this;}
   KN& operator -=(const if_arth_KN_<R> & u)  
        { KN_<R>::operator-=(u);return *this;}
   KN& operator -=(const Mul_KNM_KN_<R> & u) 
        { KN_<R>::operator-=(u);return *this;}
 
   KN& operator +=(const_R a)  
        { KN_<R>::operator += (a);return *this;}
   KN& operator += (const KN_<R>& a)  
        { KN_<R>::operator+= (a);return *this;}
   KN& operator +=(const Add_KN_<R> & u)  
        { KN_<R>::operator+=(u);return *this;}
   KN& operator +=(const DotStar_KN_<R> & u)  
        { KN_<R>::operator+=(u);return *this;}
   KN& operator +=(const DotSlash_KN_<R> & u)  
        { KN_<R>::operator+=(u);return *this;}
   KN& operator +=(const Sub_KN_<R> & u)  
        { KN_<R>::operator+=(u);return *this;}
   KN& operator +=(const Mulc_KN_<R> & u)  
        { KN_<R>::operator+=(u);return *this;}
   KN& operator +=(const Add_Mulc_KN_<R> & u)  
        { KN_<R>::operator+=(u);return *this;}
   KN& operator +=(const if_arth_KN_<R> & u)  
        { KN_<R>::operator+=(u);return *this;}
   KN& operator +=(const Mul_KNM_KN_<R> & u) 
        { KN_<R>::operator+=(u);return *this;}
        

   KN& operator/=(const_R a)  
        { KN_<R>::operator/=(a);return *this;}
   KN& operator /= (const KN_<R>& a)  
        { KN_<R>::operator/= (a);return *this;}
   KN& operator /=(const Add_KN_<R> & u)  
        { KN_<R>::operator/=(u);return *this;}
   KN& operator /=(const Sub_KN_<R> & u)  
        { KN_<R>::operator/=(u);return *this;}
   KN& operator /=(const Mulc_KN_<R> & u)  
        { KN_<R>::operator/=(u);return *this;}
   KN& operator /=(const Add_Mulc_KN_<R> & u)  
        { KN_<R>::operator/=(u);return *this;}
   KN& operator /=(const if_arth_KN_<R> & u)  
        { KN_<R>::operator/=(u);return *this;}
        
   KN& operator /=(const Mul_KNM_KN_<R> & u) 
        { KN_<R>::operator/=(u);return *this;}
        
   KN& operator*=(const_R a)  
        { KN_<R>::operator*=(a);return *this;}
   KN& operator*=(const KN_<const_R>& a)  
        { KN_<R>::operator*= (a);return *this;}
   KN& operator *=(const Add_KN_<R> & u)  
        { KN_<R>::operator*=(u);return *this;}
   KN& operator *=(const Sub_KN_<R> & u)  
        { KN_<R>::operator*=(u);return *this;}
   KN& operator *=(const Mulc_KN_<R> & u)  
        { KN_<R>::operator*=(u);return *this;}
   KN& operator *=(const Add_Mulc_KN_<R> & u)  
        { KN_<R>::operator*=(u);return *this;}
   KN& operator *=(const if_arth_KN_<R> & u)  
        { KN_<R>::operator*=(u);return *this;}
   KN& operator *=(const Mul_KNM_KN_<R> & u) 
        { KN_<R>::operator*=(u);return *this;}
        
  
  template<class I,class T> KN& operator =   (const KN_ITAB<T ,I> & ui)
     {  KN_<R>::operator =(ui);  return *this;}
  template<class I,class T> KN& operator +=   (const KN_ITAB<T ,I> & ui)
     {  KN_<R>::operator +=(ui);  return *this;}
  template<class I,class T> KN& operator -=   (const KN_ITAB<T ,I> & ui)
     {  KN_<R>::operator -=(ui);  return *this;}
  template<class I,class T> KN& operator *=   (const KN_ITAB<T ,I> & ui)
     {  KN_<R>::operator *=(ui);  return *this;}
  template<class I,class T> KN& operator /=   (const KN_ITAB<T ,I> & ui)
     {  KN_<R>::operator /=(ui);  return *this;}
        
        
  //  two opertor to cast to an array of constant      
//    operator KN_<const_R> & ()  
//          { return *  (KN_<const_R>*) this;}
//    operator KN_<const_R> const & ()  const 
//          { return *(const KN_<const_R>*) this;}
//    operator KN<const_R> & () 
//          { return   (KN<const_R> &) *this;}
//    operator KN<const_R> const & ()  const 
//          { return (const KN<const_R>& ) *this;}
  void init(long nn) {this->n=nn;this->step=1;this->next=-1;this->v=new R[nn];}
  void init() {this->n=0;this->step=1;this->next=-1;this->v=0;}
  void init(const KN_<R> & a){init(a.N()); operator=(a);}
  void resize(long nn) {
    if ( nn != this->n) 
     {
       R *vo=this->v;
       long no=std::min(this->n,nn), so=this->step;
       ShapeOfArray::init(nn);
       this->v=new R[this->n];
       // copy
       if(this->v && vo) 
         for(long i=0,j=0;j<no;i++,j+=so) 
           this->v[i]=vo[j]; 
        delete [] vo;} }
  void destroy(){delete [] this->v; this->v=0;this->n=0;}
};

//  Array with 2 indices
//  ---------------------

template<class R>
class KNM: public KNM_<R>{ public:

  KNM(long n,long m) 
        :KNM_<R>(new R[n*m],n,m){}
   KNM(const KNM<R> & u)  // PB si stepi ou stepj nulle
        :KNM_<R>(new R[u.size()],u.N(),u.M()) 
       { KN_<R>::operator=(u);}
  explicit KNM(const KNM_<R> & u)
        :KNM_<R>(new R[u.size()],u.N(),u.M())