matricecreuse_tpl.hpp

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

#ifndef  MatriceCreuse_tpl_
#ifndef  MatriceCreuse_h_
#include "MatriceCreuse.hpp"
#include <limits>
#include <set>
#include <list>
#include <map>
#endif

#ifndef __MWERKS__
// test blas 
//  on MacOS9 under MWERKS
//  cblas_ddot macos-9 is not 
#ifdef HAVE_CBLAS_H
extern "C" {
#include <cblas.h> 
}
#define WITHBLAS 1
#elif HAVE_VECLIB_CBLAS_H
#include <vecLib/cblas.h> 
#define WITHBLAS 1
#endif  
#endif  
#ifdef WITHBLAS 
template<class R> inline R blas_sdot(const int n,const R *sx,const int incx,const R *sy,const int  incy)
{
  R s=R();
  
  if(incx == 1 && incy == 1)
   for (int k = 0; k< n; k++)
    s += *sx++ * * sy++;
  else
   for (int k = 0; k< n; k++, sx += incx, sy += incy)
    s += *sx * *sy;   
   return  s;
}
template<> inline float blas_sdot(const int n,const  float *sx,const int incx,const float *sy,const int  incy)
{
  return cblas_sdot(n,sx,incx,sy,incy);
}
template<> inline double blas_sdot(const int n,const  double *sx,const int incx,const double *sy,const int  incy)
{
  return cblas_ddot(n,sx,incx,sy,incy);
}
template<> inline  complex<double> blas_sdot(const int n,const  complex<double> *sx,const int incx,const complex<double> *sy,const int  incy)
{
  complex<double> s;
   cblas_zdotu_sub(n,(const void *)sx,incx,(const void *)sy,incy,( void *)&s);
  return s;
}
template<> inline  complex<float> blas_sdot(const int n,const  complex<float> *sx,const int incx,const complex<float> *sy,const int  incy)
{
  complex<float> s;
   cblas_zdotu_sub(n,(const void *)sx,incx,(const void *)sy,incy,( void *)&s);
  return s;
}
#endif
//  end modif FH
using Fem2D::HeapSort;
using std::numeric_limits;

//  -----------
template<class FElement>
inline int  BuildMEK_KK(const int l,int *p,int *pk,int *pkk,const FElement * pKE,const FElement*pKKE)
{
  // routine  build  les array p, pk,pkk 
  // which return number of df int 2 element pKE an pKKE
  // max l size of array p, pk, pkk
  // p[i] is the global number of freedom
  // pk[i] is is the local number in pKE ( -1 if not in pKE element)
  // pkk[i] is is the local number in pKKE ( -1 if not in pKKE element)
  //  remark, if pKKE = 0 => 
     const FElement (*pK[2])={pKE,pKKE};

     int ndf=0; // number of dl
     int * qk=pk, *qkk=pkk;
     for (int k=0;k<2;k++)
      if(pK[k]) 
       {
         if(k) Exchange(qk,qkk);
         const FElement& FEK=*pK[k];
         int nbdf =FEK.NbDoF();
         
         for (int ii=0;ii<nbdf;ii++)
          {
           p[ndf] = FEK(ii); // copy the numbering 
           qk[ndf] = ii;
           qkk[ndf++] = -1;
          } // end for ii
         } 
      ffassert(ndf <=l);
   // compression suppression des doublons
       Fem2D::HeapSort(p,pk,pkk,ndf);
       int k=0;  
        for(int ii=1;ii<ndf;++ii)
          if (p[k]==p[ii]) // doublons 
            { 
              if (pkk[ii]>=0) pkk[k]=pkk[ii];
              if (pk[ii]>=0) pk[k]=pk[ii];
              assert(pk[k] >=0 && pkk[k]>=0);
            }
           else { // copy 
              p[++k] =p[ii];
              pk[k]=pk[ii];
              pkk[k]=pkk[ii];
             }
        ndf=k+1; 
         
   return ndf;
} //  BuildMEK_KK

template<class R,class FES>
void MatriceElementairePleine<R,FES>::call(int k,int ie,int label,void * stack) {
  for (int i=0;i<this->lga;i++) 
     this->a[i]=0;
  if(this->onFace)
    { 
     throwassert(faceelement);
     const Mesh &Th(this->Vh.Th);
     
     int iie=ie,kk=Th.ElementAdj(k,iie);
     if(kk==k|| kk<0) kk=-1;
     if ( &this->Vh == &this->Uh)
      {
       FElement Kv(this->Vh[k]);
       if(kk<0)
        { // return ; // on saute ????  bof bof 
	  this->n=this->m=BuildMEK_KK<FElement>(this->lnk,this->ni,this->nik,this->nikk,&Kv,0);
         int n2 =this->m*this->n; 
         for (int i=0;i<n2;i++) this->a[i]=0;
         faceelement(*this,Kv,Kv,Kv,Kv,this->data,ie,iie,label,stack);
        }
        else
        {
         FElement KKv(this->Vh[kk]);
         this->n=this->m=BuildMEK_KK<FElement>(this->lnk,this->ni,this->nik,this->nikk,&Kv,&KKv);
        
         
         faceelement(*this,Kv,KKv,Kv,KKv,this->data,ie,iie,label,stack);

        }
      }
     else 
      {
        ERREUR("A FAIRE/ TO DO  (see F. hecht) ", 0); 
        ffassert(0); // a faire F. Hecht desole 
       
      }
   }
  else {
  throwassert(element);
  const FElement&Kv(this->Vh[k]);
  int nbdf =Kv.NbDoF();
  for (int i=0;i<nbdf;i++)
     this->ni[i] = Kv(i); // copy the numbering 
  this->m=this->n=nbdf;  

  if(this->ni != this->nj) { // 
    const FElement&Ku(this->Uh[k]);
    int nbdf =Ku.NbDoF();
    for (int i=0;i<nbdf;i++)
      this->nj[i] = Ku(i); // copy the numbering 
     this->m=nbdf;
     int n2 =this->m*this->n; 
     for (int i=0;i<n2;i++) this->a[i]=0;
     element(*this,Ku,Kv,this->data,ie,label,stack);
  }
  else 
    {
     int n2 =this->m*this->n;
     for (int i=0;i<n2;i++) this->a[i]=0;
     element(*this,Kv,Kv,this->data,ie,label,stack);
   // call the elementary mat 
    }  
  }  
}

template<class R,class FES>
void MatriceElementaireSymetrique<R,FES>::call(int k,int ie,int label,void * stack) {
  // mise a zero de la matrice elementaire, plus sur
  for (int i=0;i<this->lga;i++) 
    this->a[i]=0;
  if(this->onFace)
    { 
      ffassert(0); // a faire 
    }
  else {
    
    if (k< this->Uh.Th.nt)
      {
	throwassert(element);
	const FElement K(this->Uh[k]);
	int nbdf =K.NbDoF();
	for (int i=0;i<nbdf;i++)
	  this->ni[i] = K(i); // copy the numbering 
	this->m=this->n = nbdf; 
	
	element(*this,K,this->data,ie,label,stack); 
      }// call the elementary mat 
    else
      {
	ffassert(0); // remove code for the 3d 
	/*
	throwassert(mortar);
	{
	  const FMortar K(&(this->Uh),k);
	  int nbdf = K.NbDoF();
	  for (int i=0;i<nbdf;i++)
	    this->ni[i] = K(i); // copy the numbering 
	  this->m=this->n = nbdf; 
	  // mise a zero de la matrice elementaire, plus sur
	  
	  mortar(*this,K,stack);}
	*/
      }
  }
}
  
template<class R>
MatriceProfile<R>::~MatriceProfile() {
  if(!this->dummy) 
    { //cout << " del mat profile " << endl ;
    if (U && (U !=L))  delete [] U;
    if (D)  delete [] D;
    if (L)  delete [] L;
    if (pU && (pU != pL)) delete [] pU;
    if (pL) delete [] pL;
    //cout << " dl de MatriceProfile " << this << endl;
    }
}
template<class R>
int MatriceProfile<R>::size() const {
  int s = sizeof(MatriceProfile<R>);
  if (D) s += this->n*sizeof(R);
  if (pL) s += this->n*sizeof(int);
  if (pU && (pU != pL)) s += this->n*sizeof(int);
  if (L) s += pL[this->n]*sizeof(int);
  if (U && (U != L)) s += pU[this->n]*sizeof(int);
  return s;
}
/*
template<class R>
int MatriceProfile<R>::MatriceProfile(const MatriceProfile<RR> & A )
  : MatriceCreuse<R>(A.n,A.m,0)
  {
    
    typefac=A.typefac;
    pL=  docpy<int,int>(A.pL,n+1);
    D = docpy<R,RR>(A.D,n);
    if ( A.pL == A.pU ) pU=pL;
    else pU=  docpy<int,int>(A.pU,m+1);
    
      L= docpy<R,RR>(A.L,pL[n]);
      
    if ( A.L == A.U ) U=L;
    else  U= docpy<R,RR>(A.U,pU[m]);
    
  
  }*/
template<class R>
  MatriceMorse<R> *MatriceProfile<R>::toMatriceMorse(bool transpose,bool copy) const 
  {
  // A FAIRE;
    ffassert(0); // TODO
   return 0;
  }
  
 inline pair<int,int> ij_mat(bool trans,int ii00,int jj00,int i,int j) {
  // warning trans sub  matrix and not the block. 
  return trans ? make_pair<int,int>(j+ii00,i+jj00)
                :  make_pair<int,int>(i+ii00,j+jj00) ; }
 
template<class R>
bool MatriceProfile<R>::addMatTo(R coef,std::map< pair<int,int>, R> &mij,bool trans,int ii00,int jj00,bool cnj)
{
   double eps0=numeric_limits<double>::min();
 if( norm(coef)<eps0) return  L == U ;
 int i,j,kf,k;
  if(D)
   for( i=0;i<this->n;i++)
    if( norm(D[i])>eps0)
     mij[ij_mat(trans,ii00,jj00,i,i)] += coef*(cnj? conj(D[i]) : D[i]);
   else
   for(int i=0;i<this->n;i++) // no dia => identity dai
     mij[ij_mat(trans,ii00,jj00,i,i)] += coef;
     
 if (L && pL )    
   for (kf=pL[0],i=0;  i<this->n;   i++  )  
     for ( k=kf,kf=pL[i+1], j=i-kf+k;   k<kf; j++,  k++  )
        if(norm(L[k])>eps0)
        mij[ij_mat(trans,ii00,jj00,i,j)]= coef*(cnj? conj(L[k]) : L[k]);
 if (U && pU)     
   for (kf=pU[0],j=0;  j<this->m;  j++)  
     for (k=kf,kf=pU[j+1], i=j-kf+k;   k<kf; i++,  k++  )
      if(norm(U[k])>eps0)
        mij[ij_mat(trans,ii00,jj00,i,j)]= coef*(cnj? conj(U[k]) : U[k]);
 return L == U ; // symetrique               
}
template<class R>
MatriceProfile<R>::MatriceProfile(const int nn,const R *a)
  :MatriceCreuse<R>(nn,nn,0),typefac(FactorizationNO)
{
  int *pf = new int [this->n+1];
  int i,j,k;
  k=0;
  for (i=0;i<=this->n;k+=i++)
    {
      pf[i]=k;
      //  cout << " pf " << i<< " = " << k  << endl;
    }
  ffassert( pf[this->n]*2 == this->n*(this->n-1));
  pU = pf; // pointeur profile U
  pL = pf; // pointeur profile L
  U = new R[pf[this->n]];
  L = new R[pf[this->n]];
  D = new R[this->n];  
  const R *aij=a;
  for (i=0;i<this->n;i++)
    for (j=0;j<this->n;j++)
      if      (j<i)   L[pL[i+1]-i+j] = *aij++;
      else if (j>i)   U[pU[j+1]-j+i] = *aij++;
      else            D[i] = *aij++;
}

template<class R>
template<class FESpace>
MatriceProfile<R>::MatriceProfile(const FESpace & Vh,bool VF) 
  :MatriceCreuse<R>(Vh.NbOfDF,Vh.NbOfDF,0),typefac(FactorizationNO)
{
   // for galerkine discontinue ....
   // VF : true=> Finite Volume matrices (change the stencil) 
   // VF = false => Finite element 
   // F. Hecht nov 2003
   // -----
  this->dummy=0;
  this->n = this->m = Vh.NbOfDF;
  int i,j,k,ke,ie,mn,jl,iVhk;
  int itab,tabk[5]; 
  int *pf = new int [this->n+1];
  for (i=0;i<this->n;i++)  pf[i]=0;
  for (ke=0;ke<Vh.NbOfElements;ke++)
    { 
      itab=0;
      tabk[itab++]=ke;
      if(VF) itab += Vh.Th.GetAllElementAdj(ke,tabk+itab);
      tabk[itab]=-1;    
      mn = this->n;
      for( k=tabk[ie=0]; ie <itab; k=tabk[++ie])
	{ iVhk=(int) Vh(k);
        for (jl=0;jl<iVhk;jl++) // modif Oct 2008 valgrind
	  { 
	    j=Vh(k,jl) ;
	    mn = Min ( mn , Vh.FirstDFOfNode(j) ) ;}
        }
       //for( k=tabk[ie=0]; ie <itab; k=tabk[++ie])
        { k=ke; // bof bof a verifier finement .... FH
	  iVhk=(int) Vh(k);  
        //for (j=Vh(k,jl=0);jl<(int) Vh(k);j=Vh(k,++jl)) 
	for (jl=0;jl<iVhk;jl++) // modif Oct 2008 valgrind
	     {
		 j=Vh(k,jl);
	      int df1 = Vh.LastDFOfNode(j);
	      for (int df= Vh.FirstDFOfNode(j);  df < df1; df++  )
	       pf[df] = Max(pf[df],df-mn);
	     }
	     }
    }
  int l =0;
  for (i=0;i<this->n;i++)  {int tmp=l;l += pf[i]; pf[i]=tmp;}
  pf[this->n] = l;
  if(verbosity >3) 
    cout << "  -- SizeOfSkyline =" <<l << endl;

  pU = pf; // pointeur profile U
  pL = pf; // pointeur profile L
  D = 0; // diagonal
  U = 0; // upper part
  L = 0; // lower part 
}

template<class R>
void MatriceProfile<R>::addMatMul(const KN_<R> &x,KN_<R> &ax) const 
{if (x.n!= this->n ) ERREUR(MatriceProfile MatMut(xa,x) ," longueur incompatible x (in) ") ;
 if (ax.n!= this->n ) ERREUR(MatriceProfile MatMut(xa,x) ," longueur incompatible ax (out)") ;
 int i,j,k,kf;
 ffassert(this->n == this->m);
 if (D) 
   for (i=0;i<this->n;i++) 
     ax[i] += D[i]*x[i];
 else
   for (i=0;i<this->n;i++) // no dia => identyty dai
     ax[i] +=x[i];
      
 if (L && pL )    
   for (kf=pL[0],i=0;  i<this->n;   i++  )  
     for ( k=kf,kf=pL[i+1], j=i-kf+k;   k<kf; j++,  k++  )
       ax[i] += L[k]*x[j],throwassert(i>=0 && i <this->n && j >=0 && j < this->m && k>=0 && k < pL[this->n]);
       
 if (U && pU)     
   for (kf=pU[0],j=0;  j<this->m;  j++)  
     for (k=kf,kf=pU[j+1], i=j-kf+k;   k<kf; i++,  k++  )
       ax[i] += U[k]*x[j],throwassert(i>=0 && i <this->n && j >=0 && j < this->m &&  k>=0 && k < pU[this->n]);
 

}


template<class R>
void MatriceProfile<R>::operator=(const R & v) {
  if(v!=R())
    { cerr << " Mise a zero d'une matrice MatriceProfile<R>::operator=(R v) uniquement v=" << v << endl;
    throw(ErrorExec("exit",1));
    }
  typefac = FactorizationNO;
  delete [] U;
  delete [] L;
  delete [] D;
  U=L=D=0;
}
template<class R>
MatriceCreuse<R>  & MatriceProfile<R>::operator +=(MatriceElementaire<R> & me) {
  int il,jl,i,j,k;
  int * mi=me.ni, *mj=me.nj;
  if (!D)  // matrice vide 
    { D  = new R[this->n];
    L  = pL[this->n] ? new R[pL[this->n]] :0 ;
    for (i =0;i<this->n;i++) D[i] =0;
    for (k =0;k<pL[this->n];k++) L[k] =0;
    switch (me.mtype) {
    case MatriceElementaire<R>::Full :     
      U  = pU[this->n] ? new R[pU[this->n]] : 0;
      for (k =0;k<pU[this->n];k++) U[k] =0;
      break;
    case MatriceElementaire<R>::Symmetric :     
      U = L; 
      break;
    default:
      cerr << "Big bug type MatriceElementaire unknown" << (int) me.mtype << endl;
      throw(ErrorExec("exit",1));
      break; 
    }
    }
  R * al = me.a; 
  switch (me.mtype) {
  case MatriceElementaire<R>::Full : //throwassert(L !=U);
    for (il=0; il<me.n; ++il) // modif overflow FH win32  oct 2005
     { i=mi[il];
      for ( jl=0; jl< me.m ; ++jl,++al)  // modif overflow FH
        { j=mj[jl] ;
	  if      (j<i)  L[ pL[i+1] - (i-j) ] += *al;
	  else if (j>i)  U[ pU[j+1] - (j-i) ] += *al;
	  else           D[i] += *al;}}
    break;