matricecreuse.hpp

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

#ifndef MatriceCreuse_h_
#define MatriceCreuse_h_

template<class T>
 T Square(const T & r){return r*r;}

 
#ifdef HAVE_LIBUMFPACK_XXXXXXXXXXXXX
extern "C" {
#ifdef HAVE_UMFPACK_H
#include <umfpack.h>
#else
#ifdef HAVE_UMFPACK_UMFPACK_H
#include <umfpack/umfpack.h>
#else
#ifdef HAVE_BIG_UMFPACK_UMFPACK_H
#include <UMFPACK/umfpack.h>
#else
#ifdef HAVE_UFSPARSE_UMFPACK_H
#include <ufsparse/umfpack.h>
#else
#ifdef HAVE_SUITESPARSE_UMFPACK_H
#include <suitesparse/umfpack.h>
#else

  // Defaults to a local version of the UMFPACK headers
#include "../../download/include/umfpack.h"

#endif // HAVE_SUITESPARSE_UMFPACK_H
#endif // HAVE_UFSPARSE_UMFPACK_H
#endif // HAVE_BIG_UMFPACK_UMFPACK_H
#endif // HAVE_UMFPACK_UMFPACK_H
#endif // HAVE_UMFPACK_H
}
#endif

#include "RNM.hpp"
#include "fem.hpp"
#include "FESpace.hpp" 
#include "DOperator.hpp"
#include "QuadratureFormular.hpp" 

using  Fem2D::Mesh;
using  Fem2D::FESpace;
using Fem2D::FElement;
using Fem2D::baseFElement;
using Fem2D::FMortar;
using Fem2D::TypeOfMortar;
using Fem2D::QuadratureFormular;
using Fem2D::QuadratureFormular1d;
using Fem2D::QuadratureFormular_T_5;
using Fem2D::QF_GaussLegendre3;
const double  EPSILON=1e-20;
using Fem2D::onWhatIsEdge;

//#define APROGRAMMER(a) {cerr << "A PROGRAMMER "  #a << endl; exit (1) ;}
#define ERREUR(a,b) {cerr << "ERREUR " #a<< b <<endl; throw(ErrorExec("FATAL ERREUR dans " __FILE__  "\n" #a  " line: ",__LINE__)) ;}

template <class R> class MatriceCreuse; 
template <class R> class MatriceElementaire; 
template <class R,class FESpace> class MatriceElementaireSymetrique;
template <class R,class FESpace> class MatriceElementairePleine;
template <class R> class MatriceMorse;
template <class R> class MatriceProdTensoriel;

//template <class R> R Square(R x){ return x*x;}

template <class T> T* docpyornot(bool nocpy,T* p,int n)
{ 
  T * r=p;
   if( !nocpy) { // do copy 
      r= new T[n]; ffassert(r);
      for(int i=0;i<n;i++) 
        r[i]=p[i];
      }
   return r;
 }
 template <class T,class TT> T* docpy(TT* p,int n)
{ 
   T * r=0;
   if(p && n) { // do copy 
      r= new T[n]; ffassert(r);
      for(int i=0;i<n;i++) 
        r[i]=T(p[i]); // pour cadna ???? FH 
      }
   return r;
 }



template <class R> 
class MatriceElementaire {

public:
  enum TypeOfMatriceElementaire {Full=1,Symmetric=2};


       
  int lga;  // size of array a    
  R* a;          // array  coef --
  int *ni,*nj;   //  list of df   
  // to build matrice on face or edge -----

  
  int n,m;       // n,m number of df
  const TypeOfMatriceElementaire mtype;
  KN<double> data; // to store value of basic function 

  const bool onFace ; //  true if do int on face or edge with jump (VF or GD : Galerkin Discontinus)
  // in with case  add ... 
  const int lnk; // size of the 4 next array
  int *nik,*nikk;  //  number of df in element   k,kk for VF and GD methode 
  int *njk,*njkk;  //  number of df in element   k,kk  for VF and GD methode


  MatriceElementaire(int datasize,int llga
                     ,int *nnj,int * nni,TypeOfMatriceElementaire t=Full)
    
    :   lga(llga),a(new R[lga]),
        ni(nni),nj(nnj),n(0),m(0),mtype(t),data(datasize),
        onFace(false),lnk(0),nik(0),nikk(0),njk(0),njkk(0) 
        {}
       

 //  for discontinous Galerkine method
  MatriceElementaire(int datasize,int llga,int *nni,
                     int lk,
                     TypeOfMatriceElementaire t=Symmetric
                     ) 
    :
    lga(llga),a(new R[lga]),
    ni(nni),nj(nni),n(0),m(0),mtype(t),data(datasize*(lk?2:1)) ,
       onFace(lk!=0),
       lnk(lk),
       nik(lk? new int[lk*2]:0),
       nikk(nik+lk),
       njk(nik),
       njkk(nik+lk)
       { ffassert(lk>=0);}
 
  virtual ~MatriceElementaire() {
    if(ni != nj) 
      delete [] nj;
         
    delete [] ni;
    delete [] a;
    if ( nik) delete[] nik;
     }
      
  virtual R & operator() (int i,int j) =0;
  virtual void call(int ,int ie,int label,void * data) =0;  // 
  const LinearComb<pair<MGauche,MDroit>,C_F0> * bilinearform;
  
  MatriceElementaire & operator()(int k,int ie,int label,void * s=0) {
    call(k,ie,label,s);
    return *this;}
};

template <class FES> 
class MatDataFES { 
public:
  typedef FES FESpace;
  typedef typename FESpace::FElement FElement;

  typedef typename  FESpace::QFElement QFElement; 
  typedef typename  FESpace::QFBorderElement QFBorderElement; 
  CountPointer<const FESpace> cUh,cVh;
  const FESpace &Uh;
  const FESpace &Vh;
  const QFElement & FIT;
  const QFBorderElement & FIE;
  MatDataFES(const FESpace & UUh,const QFElement & fit, const QFBorderElement & fie)
    :Uh(UUh),Vh(UUh),FIT(fit),FIE(fie) {}
  MatDataFES(const FESpace & UUh,const FESpace & VVh,const QFElement & fit, const QFBorderElement & fie)
    :Uh(UUh),Vh(VVh),FIT(fit),FIE(fie) {}
    

};

template <class R,class FES> 
class MatriceElementaireFES :   public MatDataFES<FES> ,   public MatriceElementaire<R> 
{  
public:
  typedef MatriceElementaire<R> MElm ;
  using MElm::Full;
  using MElm::Symmetric;

  typedef typename MElm::TypeOfMatriceElementaire TypeOfMatriceElementaire;
  typedef FES FESpace;

  typedef typename  FESpace::FElement FElement; 
  typedef typename  FESpace::QFElement QFElement; 
  typedef typename  FESpace::QFBorderElement QFBorderElement; 

  MatriceElementaireFES(const FESpace & UUh,const FESpace & VVh,int llga
			,int *nnj,int * nni,TypeOfMatriceElementaire t=Full,
			const QFElement & fit=*QFElement::Default,
			const QFBorderElement & fie =*QFBorderElement::Default) 
                     
    :
    MatDataFES<FES>(UUh,VVh,fit,fie),
    MatriceElementaire<R>(UUh.esize()+VVh.esize(),llga,nnj,nni,t)
  {}
       
  MatriceElementaireFES(const FESpace & UUh,int llga,int *nni,
			TypeOfMatriceElementaire t=Symmetric,
			const QFElement & fit=*QFElement::Default,
			const QFBorderElement & fie =*QFBorderElement::Default)
    :
    MatDataFES<FES>(UUh,UUh,fit,fie),
    MatriceElementaire<R>(UUh.esize(),llga,nni,nni,t)
  {}

  //  for discontinous Galerkine method
  MatriceElementaireFES(const FESpace & UUh,int llga,int *nni,
			int lk,
			TypeOfMatriceElementaire t=Symmetric,
			const QFElement & fit=*QFElement::Default,
			const QFBorderElement & fie =*QFBorderElement::Default) 
    :
    MatDataFES<FES>(UUh,UUh,fit,fie),
    MatriceElementaire<R>(UUh.esize(),llga,nni,lk,t)
  {}
  ~MatriceElementaireFES() {}
  const LinearComb<pair<MGauche,MDroit>,C_F0> * bilinearform;
  
  MatriceElementaireFES & operator()(int k,int ie,int label,void * s=0) {
    this->call(k,ie,label,s);
    return *this;}
};

template <class R,class FES=FESpace> 
class MatriceElementairePleine:public MatriceElementaireFES<R,FES> {

  /* --- stockage --
     //  n = 4 m = 5
     //  0  1  2  3  4
     //  5  6  7  8  9
     // 10 11 12 13 14
     // 15 16 17 18 19
     ------------------*/

public:
  typedef FES FESpace;
  typedef typename  FESpace::Mesh Mesh;
  typedef typename  FESpace::QFElement QFElement;
  typedef typename  FESpace::QFBorderElement QFBorderElement;
  typedef typename  FESpace::FElement FElement; 

  R & operator() (int i,int j) {return this->a[i*this->m+j];}
  // MatPleineElementFunc element;
  void  (* element)(MatriceElementairePleine &,const FElement &,const FElement &, double*,int ie,int label,void *) ; 
   void  (* faceelement)(MatriceElementairePleine &,const FElement &,const FElement &,const FElement &,const FElement &, double*,int ie,int iee, int label,void *) ; 
 void call(int k,int ie,int label,void *);
  
  MatriceElementairePleine & operator()(int k,int ie,int label,void * stack=0)
  {call(k,ie,label,stack);return *this;}
  MatriceElementairePleine(const FESpace & VVh,
                           const QFElement & fit=*QFElement::Default,
                           const QFBorderElement & fie =*QFBorderElement::Default)  
    :MatriceElementaireFES<R,FES>(VVh,
			Square(VVh.MaximalNbOfDF()),
			new int[VVh.MaximalNbOfDF()],this->Full,fit,fie),
    element(0),faceelement(0) {}
 
   //  matrice for VF or Galerkin Discontinus
   MatriceElementairePleine(const FESpace & VVh,bool VF,
                           const QFElement & fit=*QFElement::Default,
                           const QFBorderElement & fie =*QFBorderElement::Default)  
     :MatriceElementaireFES<R,FES>(VVh,
			Square(VVh.MaximalNbOfDF()*2),
			new int[VVh.MaximalNbOfDF()*2],
			VF?VVh.MaximalNbOfDF()*2:0,
			this->Full,fit,fie),			
    element(0),faceelement(0) {}

  MatriceElementairePleine(const FESpace & UUh,const FESpace & VVh,
                               const QFElement & fit=*QFElement::Default,
                               const QFBorderElement & fie =*QFBorderElement::Default) 
    :MatriceElementaireFES<R,FES>(UUh,VVh,
				  UUh.MaximalNbOfDF()*VVh.MaximalNbOfDF(),
				  new int[UUh.MaximalNbOfDF()],
				  new int[VVh.MaximalNbOfDF()],this->Full,fit,fie),
     element(0),faceelement(0) {}

}; 

template <class R,class FES=FESpace> 
class MatriceElementaireSymetrique:public MatriceElementaireFES<R,FES> {



  // --- stockage --
  //   0
  //   1 2
  //   3 4 5
  //   6 7 8 9
  //  10 . . . .
  //

public:
  typedef FES FESpace;
  typedef typename  FESpace::Mesh Mesh;
  typedef typename  FESpace::QFElement QFElement;
  typedef typename  FESpace::QFBorderElement QFBorderElement;
  typedef typename  FESpace::FElement FElement; 

  R & operator()(int i,int j) 
  {return j < i ? this->a[(i*(i+1))/2 + j] : this->a[(j*(j+1))/2 + i] ;}
  void (* element)(MatriceElementaireSymetrique &,const FElement &, double*,int ie,int label,void *) ; 
  void (* mortar)(MatriceElementaireSymetrique &,const FMortar &,void *) ; 
  void call(int k,int ie,int label,void * stack);
  MatriceElementaireSymetrique(const FESpace & VVh,
                               const QFElement & fit=*QFElement::Default,
                               const QFBorderElement & fie =*QFBorderElement::Default) 
    :MatriceElementaireFES<R,FES>(
           VVh,
	   int(VVh.MaximalNbOfDF()*(VVh.MaximalNbOfDF()+1)/2),
	   new int[VVh.MaximalNbOfDF()],this->Symmetric,
       fit,fie),
       element(0),mortar(0) {}
  MatriceElementaireSymetrique & operator()(int k,int ie,int label,void * stack=0) 
  {this->call(k,ie,label,stack);return *this;};
};


template <class R>  class MatriceProfile;

//  classe modele pour matrice creuse
//  ---------------------------------
template <class R> 
class MatriceCreuse : public RefCounter,public VirtualMatrice<R> {
public:
  MatriceCreuse(int NbOfDF,int mm,int ddummy)
         : VirtualMatrice<R>(NbOfDF,mm),n(NbOfDF),m(mm),dummy(ddummy){}
  MatriceCreuse(int NbOfDF)
         : VirtualMatrice<R>(NbOfDF),n(NbOfDF),m(NbOfDF),dummy(1){}
  int n,m,dummy;
  virtual int size() const =0;

  virtual MatriceCreuse & operator +=(MatriceElementaire<R> & )=0;
  virtual void operator=(const R & v) =0; // Mise a zero 
  KN_<R> & MatMul(KN_<R> &ax,const KN_<R> &x) const { 
    ax= R();
    addMatMul(x,ax);
    return ax;}
  virtual ostream& dump (ostream&)  const =0;
  virtual void Solve(KN_<R> & x,const KN_<R> & b) const =0;
  virtual ~MatriceCreuse(){}
  virtual R & diag(int i)=0;
  virtual R & operator()(int i,int j)=0;
  virtual R * pij(int i,int j) const =0; // Add FH 
  virtual  void  resize(int n,int m)  {AFAIRE("MatriceCreuse::resize");}  // a faire dans les classe derive ... // add march 2009  FH 
  virtual MatriceMorse<R> *toMatriceMorse(bool transpose=false,bool copy=false) const {return 0;} // not 
  virtual bool addMatTo(R coef,std::map< pair<int,int>, R> &mij,bool trans=false,int ii00=0,int jj00=0,bool cnj=false)=0;
  // Add FH april 2005
  virtual R pscal(const KN_<R> & x,const KN_<R> & y) =0 ; // produit scalaire  
  virtual double psor(KN_<R> & x,const  KN_<R> & gmin,const  KN_<R> & gmax , double omega) =0;
  virtual void setdiag(const KN_<R> & x)=0 ;
  virtual void getdiag( KN_<R> & x) const =0 ;
  // end add
  virtual int NbCoef() const {return 0;};
  virtual void setcoef(const KN_<R> & x)=0 ;
  virtual void getcoef( KN_<R> & x) const =0 ;
  // Add FH oct 2005
   bool ChecknbLine(int nn) const { return n==nn;}  
   bool ChecknbColumn(int mm) const { return m==mm;} 

  // end ADD

};

template <class R> 
inline ostream& operator <<(ostream& f,const MatriceCreuse<R> & m) 
    {return m.dump(f);}

template <class R>