rnm.hpp

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

// ********** DO NOT REMOVE THIS BANNER **********
// ORIG-DATE:    29 fev 2000  
// -*- Mode : c++ -*-
//
// SUMMARY  : array modelisation 
// USAGE    : LGPL      
// ORG      : LJLL Universite Pierre et Marie Curie, Paris,  FRANCE 
// AUTHOR   : Frederic Hecht
// E-MAIL   : [e-mail]
//

/*
 

 
 Freefem++ is free software; you can redistribute it and/or modify
 it under the terms of the GNU Lesser General Public License as published by
 the Free Software Foundation; either version 2.1 of the License, or
 (at your option) any later version.
 
 Freefem++  is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU Lesser General Public License for more details.
 
 You should have received a copy of the GNU Lesser General Public License
 along with Freefem++; if not, write to the Free Software
 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 

 */
#ifndef KNM_H_
#define KNM_H_
// version sept 2008 FH.
// ----------------------
// une tentative qui ne marche pas 
// de tableau constant 
#include <complex>
#include <iostream>
#include <iomanip>
#include <cmath>
#include <cassert>


using namespace std;
#define const_R   R

#include <cstdlib>
inline void Check_Kn(const char * str,const char * file,int line)
{
 cerr << "CHECK_KN: " << str << " in file: " << file << ", line " << line <<endl;
 assert(0);
 abort();
}

#define K_bigassert(i)  if (!(i)) Check_Kn(#i,__FILE__,__LINE__);

#ifdef CHECK_KN

#define K_throwassert(i)  if (!(i)) Check_Kn(#i,__FILE__,__LINE__);

#else
#define K_throwassert(i) 
#endif
// version du 29 fev 2000  
//  correction   for (... lj++,ui++  qui apelle le produit scalaire
//  petite correction  throwassert 
// ajoute de operateur /= et *= sur des vecteurs
//   suppression de constructeur qui pose de probleme   
//   correction oper +=  ...  dans RNM_op.h ligne 56  change = en oper
// version de 25 nov 99  sans const R
// ajoute de '.' pour extraire une colonne, ou ligne  , ...
//  version du 22 nov 1999   cast to KN_<const R> 
//   version du 21 nov 1999  correction delete 
//  version du 13 nov 1999
//  version du 18 mars 99
//  F. Hecht 
// attention les indexations les indexations peuvent changer
//  puisque que l'on peut prendre la transposer d'une matrice
// tableau 
// mais ils partent de 0 
// version corrigee du 15/11/98
// version avec sous tableau  ---  mars 99 
// -------
//  remarque du 8 mars 99 FH
// class pour prendre des sous-tableau
// attention aux PB de continute dans les tableaux
// on a supposer que les tableaux multi indices pouvait est vue comme 
// un tableau continue ce qui est generalement faux quand l'on en 
// prend un sous tableau
//   exemple: un tableau 3,5 est numerote comme:
//    0  3  6  9 12
//    1  4  7 10 13
//    2  5  8 11 14
//             step
//   indexi  n 1
//   indexj  m n
//   est le sous tableau  3,3  n'est pas numeroter consecutivement 
//  
//    Donc la fonction  IsVector1() nous dit si un tableau 
//    a un 2 ou 3 indices est ou non consecutif en memoire
//    
//  --  ajoute d'une classe VirtualMatrice
// pour modeliser le produit matrice vecteur 
//  x = A*v; via des fonctions virtuelle 
//  ---------------------------------- 
//   version du 6 mars 2001 FH
//   ---  initialisation --
//   --------------------------------
//   version du 9 oct 2001 FH
//   ajoute de constructeur par defaut d'une vecteur
//   +  set , pour definir le vecteur 
//   ou l'affectation (bof bof) 
// ---------------------
//  version sep 2002
//  ajoute  operateur >> pour  KN<R> et KN_<R>
//  --------------------  
//  version  april 2003
//  ajoute un gestion auto de 
//  la fonction InternalError pour les matriceVirtuel
//  --------------------  
//   version jan 2004
//   correction pour go ++ 
//   des operateur  #=  pour les matrices et tenseurs
//  ----------------------
//   version  feb 2004 
//   v(i(:)) =  w   //  i(1:10) 
//   w=u(i(:))  //  
//   version mars 2004 make small correction
//   in  ITAB operator problem if non type R a defi 
//   -------------------
//   Modif pour version avec les Complex   mai 2004                                                                                                 
//   (u,v)  donne le produit complex utiliser dans le produit matrice vecteur
//   (u,conj(v))  donne le produit hermitiene pour le gradient conjugue
//
//   -- de fonction dans le cas real                                                                                                                 
// modif for g++ 4.0 and xlc++   mai 2005
//  adding some this-> 
//   mars 2007
// correction in operator operation:b -1*c 
// aout 2007, 
//  correct y = A*x ; when y is unset 
//  correct y += A*x ; when y is unset 
//  re-correct += sep 2007
//  add size of the matrix in VirtualMatrix class.
// ----------------
inline double  conj(const double & x){return x;}
inline float  conj(const float &x){return x;}
inline long  conj(const long &x){return x;}
inline double  real(const double &x){return x;}
inline float  real(const float &x){return x;}

namespace RNM {
template<class T> inline T Min (const T &a,const T &b){return a < b ? a : b;}
template<class T> inline T Max (const T &a,const T & b){return a > b ? a : b;}
template<class T> inline T Abs (const T &a){return a <0 ? -a : a;}

template<class T> inline void Exchange (T& a,T& b) {T c=a;a=b;b=c;}
template<class T> inline T Max (const T &a,const T & b,const T & c){return Max(Max(a,b),c);}
template<class T> inline T Min (const T &a,const T & b,const T & c){return Min(Min(a,b),c);}
// specialisation cas complex ---
template<class T> 
inline complex<T> Min(const complex<T> &a,complex<T> &b)
{ return complex<T>(min(a.real(),b.real()),min(a.imag(),b.imag()));}
template<class T> 
inline complex<T> Max(const complex<T> &a,const complex<T> &b)
{ return complex<T>(max(a.real(),b.real()),max(a.imag(),b.imag()));}

/*inline complex<double> Min(const complex<double> &a,complex<double> &b)
{ return complex<double>(Min(real(a),real(b)),Min(imag(a),imag(b)));}
inline complex<double> Max(const complex<double> &a,const complex<double> &b)
{ return complex<double>(Max(real(a),real(b)),Max(imag(a),imag(b)));}
*/ }
//  ----                                                                                                                                             

template<class R> class KNMK_ ;
template<class R> class KNM_ ;
template<class R> class KN_ ;
template<class R> class TKN_ ; // KN_ transpose 
template<class R> class notKN_ ; // KN_ not 
template<class R> class notnotKN_ ; // KN_ not not 

template<class R> class KNMK ;
template<class R> class KNM ;
template<class R> class KN ;

template<class R> class conj_KN_ ;
template<class R> class Add_KN_;
template<class R> class Sub_KN_;
template<class R> class Mulc_KN_;
template<class R> class Add_Mulc_KN_;
template<class R> class Mul_KNM_KN_; 
template<class R> class DotStar_KN_;
template<class R> class DotSlash_KN_;

template<class R> class outProduct_KN_;
template<class R> class if_KN_;
template<class R> class if_arth_KN_;
template<class R> class ifnot_KN_;
template<class R,class I> class KN_ITAB; 

template<class R,typename A,typename B> class F_KN_;


#ifndef ffassert
#define ffassert assert
#endif

// gestion des erreur interne --
#ifndef InternalError
typedef void (* TypeofInternalErrorRoutine)(const char *) ;
static TypeofInternalErrorRoutine &InternalErrorRoutinePtr()
{ 
  static TypeofInternalErrorRoutine routine=0;
  return routine;
}

static void InternalError(const char * str) {
  if (InternalErrorRoutinePtr() ) (*InternalErrorRoutinePtr())(str);
  cerr << str; 
  exit(1);
}
inline void  SetInternalErrorRoutine(TypeofInternalErrorRoutine f) 
{
  InternalErrorRoutinePtr()=f;
}
#endif
//  -- 
template<class P,class Q> 
   struct PplusQ { const P & p;const Q & q;
    PplusQ(const P & pp,const Q & qq) : p(pp),q(qq){}
    };

template<class R> 
struct  VirtualMatrice { public:
    int N,M;
    VirtualMatrice(int nn,int mm): N(nn),M(mm) {}
    VirtualMatrice(int nn): N(nn),M(nn) {}
  //  y += A x
  virtual void addMatMul(const KN_<R> &  x, KN_<R> & y) const =0; 
  virtual void addMatTransMul(const KN_<R> &  x, KN_<R> & y) const 
    { InternalError("VirtualMatrice::addMatTransMul not implemented "); }
  virtual bool WithSolver() const {return false;} // by default no solver          
  virtual void Solve( KN_<R> &  x,const KN_<R> & b) const 
    { InternalError("VirtualMatrice::solve not implemented "); } 

#ifdef VersionFreeFempp
  virtual bool ChecknbLine  (int n) const= 0; 
  virtual bool ChecknbColumn  (int m) const =0; 
#else
  virtual bool ChecknbLine  (int n) const {return true;} 
  virtual bool ChecknbColumn  (int m) const {return true;}
#endif
  struct  plusAx { const VirtualMatrice * A; const KN_<R>   x;
   plusAx( const VirtualMatrice * B,const KN_<R> &  y) :A(B),x(y) 
      { ffassert(B->ChecknbColumn(y.N())); }
    };
    
   plusAx operator*(const KN_<R> &  x) const {return plusAx(this,x);}
   
  struct  plusAtx { const VirtualMatrice * A; const KN_<R>   x;
   plusAtx( const VirtualMatrice * B,const KN_<R> &  y) :A(B),x(y) 
    {ffassert(B->ChecknbLine(y.N()));} };
    
  struct  solveAxeqb { const VirtualMatrice * A; const KN_<R>   b;
   solveAxeqb( const VirtualMatrice * B,const KN_<R> &  y) :A(B),b(y) 
    {ffassert(B->ChecknbColumn(y.N()));} };
  
  virtual ~VirtualMatrice(){} 
};

    

//template <class R> class MatriceCreuseMulKN_;
//template <class R> class MatriceCreuseDivKN_;

class ShapeOfArray;

class FromTo{ public:
  long from,to;
  FromTo(long i,long j):from(i),to(j) {K_throwassert(i<j);} 
 };
 
class SubArray{ public:
  const long n,step,start;
//  SubArray(char  nn): n(-1),step(1),start(0) {}  
 explicit SubArray(long nn,long sta=0,long s=1): n(nn),step(s),start(sta) {}
  SubArray(const FromTo& ft) : n(ft.to-ft.from+1),step(1),start(ft.from) {}
  SubArray(const ShapeOfArray & ); // all 
  long end()  const  { return start+ step*n;}
  long last() const  { return start+ step*(n-1);}
  long len1() const  { return step*(n-1);}  
};


class ShapeOfArray{ protected:
  public:

    long n;     //   n  nb of item
    long step;  //   step  nb of between 2 item
    long next;  //  the   next array of same type in matrix for subarray  
              // by default  no next 
  ShapeOfArray(const ShapeOfArray & s,long nn): n(s.n),step(s.n),next(nn) {}              
  ShapeOfArray(long nn): n(nn),step(1),next(-1) {}
  
  ShapeOfArray(long nn,long s): n(nn),step(s),next(-1) {}
  
  ShapeOfArray(long nn,long s,long nextt): n(nn),step(s),next(nextt) {}
  
  ShapeOfArray(const ShapeOfArray &old,const SubArray &sub) 
         : n(sub.n),step(old.step*sub.step),next(old.next)
          { K_throwassert((sub.last())*old.step <= old.last());} // a constructor 
          
  ShapeOfArray(const ShapeOfArray &old,long stepo,long start) 
         : n(old.n-start),step(old.step*stepo),next(old.next) 
          { K_throwassert(n>=0);}        
          
  long end()  const      { return n*step;}
  long last()      const      { return (n-1)*step;}
  long constant()  const { return step==0;}
  long index(long k) const { K_throwassert( (k>=0) && ( (k <n) || !step) );
           return step*k;}
  ShapeOfArray operator*(long stepp) const {return  ShapeOfArray(n,step*stepp,next);}  
  bool SameShape(const ShapeOfArray & a) const 
          { return  !step || !a.step || a.n == n ;} 
  long  N(const ShapeOfArray & a) { return step ? n : a.n;} // size of 2 shape 

  
// protected:
  long operator[](long k) const {  K_throwassert( (k>=0) && ( (k <n) || !step) );
           return step*k;}     
  void init(long nn,long s=1,long nextt=-1) { n=nn; step=s; next=nextt;}         
};

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

inline bool  SameShape(const ShapeOfArray & a,const ShapeOfArray & b) 
           { return  !a.step || !b.step || a.n == b.n ;} 
   
inline long N(const ShapeOfArray & a,const ShapeOfArray & b) 
           { K_throwassert(SameShape(a,b)); return  a.step ? a.n :  b.n ;}

inline SubArray::SubArray(const ShapeOfArray & s) 
           :n(s.n),step(s.step),start(0) {}
           
   

template<class R> 
ostream & operator<<(ostream & f,const KN_<const_R> & v) ;

template<class R> istream & operator>>(istream & f, KN_<R> & v);
template<class R> istream & operator>>(istream & f, KN<R> & v);


template<class R>
class KN_: public  ShapeOfArray {
protected:
  R *v;
public:
  typedef R K; // type of data 
  long N() const {return n;}
  bool unset() const { return !v;}
  void set(R * vv,int nn,int st=1,int nx=-1) {v=vv;n=nn;step=st;next=nx;}
  long size() const{return step?n*step:n;}
  operator R *() const {return v;}
  KN_(const KN_<R> & u) :ShapeOfArray(u),v(u.v){}