matricecreuse_tpl.hpp

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

  f << "# Sparse Matrix (Morse)  " << endl;
  f << "# first line: n m (is symmetic) nbcoef \n";
  f << "# after for each nonzero coefficient:   i j a_ij where (i,j) \\in  {1,...,n}x{1,...,m} \n";  
 
  f << this->n << " " << this->m << " " << symetrique << "  " << nbcoef <<endl;
  int k=lg[0];
  int pold= f.precision();
  for (int i=0;i<this->n;i++)
   { 
    
//    f << i << " : " << lg[i] <<","<< lg[i+1]-1 << " : " ;
    int ke=lg[i+1];
    for (;k<ke;k++)
      f << setw(9) << i+1 << ' ' << setw(9) << cl[k]+1 << ' ' << setprecision( 20) << a[k]<< '\n' ;
     // if (norm(a[k])) f  << cl[k] << " " << a[k]<< ", ";
     // else f  << cl[k] << " 0., " ;
   // f << endl;    
   }
   f.precision(pold);
  return f;
}
template <class R> 
inline R*  MatriceMorse<R>::pij(int i,int j) const 
 {
   if (! (i<this->n && j< this->m)) 
   throwassert(i<this->n && j< this->m);
   int i0=lg[i];
   int i1=lg[i+1]-1;
   while (i0<=i1) // dichotomie
    { 
      int im=(i0+i1)/2;
      if (j<cl[im]) i1=im-1;
      else if (j>cl[im]) i0=im+1;
      else return a+im;      
    }
   return 0;     
 }
template <class R>
template <class FESpace> 
void MatriceMorse<R>::Build(const FESpace & Uh,const FESpace & Vh,bool sym,bool VF)
{
  typedef typename FESpace::Mesh Mesh;
  
  // for galerkine discontinue ....
  // VF : true=> Finite Volume matrices (change the stencil) 
  // VF = false => Finite element 
  // F. Hecht nov 2003
  // -----
  symetrique = sym;
  this->dummy=false;
  a=0;
  lg=0;
  cl=0;
  //  bool same  = &Uh == & Vh;
  ffassert( &Uh.Th == &Vh.Th);  // same Mesh
  const Mesh & Th(Uh.Th);
  //int nbt = Th.nt;
  //int nbv = Th.nv;
  //int nbm = Th.NbMortars;
  int nbe = Uh.NbOfElements;
  int nbn_u = Uh.NbOfNodes;
  int nbn_v = Vh.NbOfNodes;
  
  KN<int> mark(nbn_v);
  KN<int> pe_u(nbn_u+1+Uh.SizeToStoreAllNodeofElement());
  //  les element du node i 
  // sont dans pe_u[k] pour k \in [ pe_u[i] , pe_u[i+1] [
  pe_u=0;
  for (int k=0;k<nbe;k++)
    { 
      int nbne=Uh(k);
      for (int in=0;in<nbne;in++)
        pe_u[(Uh(k,in)+1)]++;
   }
  int kk= nbn_u+1,kkk=kk;
  pe_u[0]=kk;
  for (int in1=1;in1<=nbn_u;in1++)
    { // in1 = in + 1
      kk += pe_u[in1];
      pe_u[in1] = kkk; // store the last of in 
      kkk=kk;
    } 
  if(verbosity>4) 
    cout <<" -- MatriceMorse<R>::Build " << kk << " " << nbn_u << " " << Uh.SizeToStoreAllNodeofElement() 
	 << " " <<  nbn_u+1+Uh.SizeToStoreAllNodeofElement() << endl;
  ffassert(kk== nbn_u+1+Uh.SizeToStoreAllNodeofElement());
  for (int k=0;k<nbe;k++)
    { 
      int nbne=Uh(k);
      for (int in=0;in<nbne;in++)
        pe_u[pe_u[(Uh(k,in)+1)]++] = k;
    }
  
  
  int color=0;
  mark=color++;
  lg = new int [this->n+1];
  ffassert(lg);
  for (int step=0;step<2;step++) 
    { 
      int ilg=0;
      lg[0]=ilg;
      int kij=0;
    for (int in=0;in<nbn_u;in++)
      {
	int nbj=0; // number of j
	int kijs=kij;
	// for all triangle contening node in
	for (int kk= pe_u[in];kk<pe_u[in+1];kk++)
	  {
	    int ke=pe_u[kk];// element of 
	    int tabk[10];
	    int ltab=0;
	    tabk[ltab++]=ke;
	    if( VF) // if Finite volume then add Triangle adj in stencil ...
	      ltab+= Th.GetAllElementAdj(ke,tabk+ltab);
	    tabk[ltab]=-1;
	    for(int ik=0,k=tabk[ik];ik<ltab;k=tabk[++ik])
	      {
		throwassert(k>=0 && k < nbe);
		int njloc = Vh(k);
		for (int jloc=0;jloc<njloc;jloc++)
		  { 
		    int  jn = Vh(k,jloc);
		    if (mark[jn] != color && (!sym ||  jn < in) ) 
		      {
			mark[jn] = color;
			int fdf=Vh.FirstDFOfNode(jn);
			int ldf=Vh.LastDFOfNode(jn);
			if (step)
			  for (int j=fdf;j<ldf;j++)
			    cl[kij++] = j;
			nbj += ldf-fdf;
		      }            
		  }} 
	  }
	int fdf=Uh.FirstDFOfNode(in);
	int ldf=Uh.LastDFOfNode(in);
	int kijl=kij;
	if (step)
	  {
	    HeapSort(cl+kijs,kij-kijs);
	    for (int i=fdf;i<ldf;i++)
	      { 
		if (i!=fdf) //  copy the ligne if not the first 
		  for (int k=kijs;k<kijl;k++)
		    cl[kij++]=cl[k]; 
		if (sym) // add block diag
		  for(int j=fdf;j<=i;j++)
		    cl[kij++]=j;            
		throwassert(kij==lg[i+1]);// verif           
	      }
	  }
	else
	  for (int i=fdf;i<ldf;i++)
	    { 
	      if (sym) ilg += ++nbj; // for the diag block
	      else ilg += nbj;             
	      lg[i+1]=ilg;
	    }
	color++; // change the color
      }
    if (step==0) { // do allocation 
      nbcoef=ilg;
      if (verbosity >3)
        cout << "  -- MorseMatrix: Nb coef !=0 " << nbcoef << endl;
      a = new R[nbcoef];
      cl = new int [nbcoef];}
    ffassert( a && cl);
    for (int i=0;i<nbcoef;i++) 
      a[i]=0;
    
   }
  
}
template<class R> inline void ConjArray( R  *v, int n) 
{
  for (int i=0;i<n;++i)
    v[i] = conj(v[i]);
}
template<> inline void ConjArray<double>(double *v, int n) {}
template<> inline void ConjArray<float>(float *v, int n) {}

template<class R>
 void  MatriceMorse<R>::dotransposition()
 {
   if(symetrique) return; 
   
   ffassert(this->dummy==false);  
   int *llg= new int[nbcoef];
   int *clg= new int[this->m+1];
   
   for (int i=0;i<this->n;i++)
     for (int k=lg[i];k<lg[i+1];k++)
        llg[k]=i;
 
  HeapSort(cl,llg,a,nbcoef);
  for(int k=0;k<this->m;k++)
    clg[k]=-1;

  // build new line end (old column)
  for(int k=0;k<nbcoef;k++)
    clg[cl[k]+1]=k+1;
      
   for(int kk=0, k=0;k<=this->m;k++)
   if (clg[k]==-1)
      clg[k]=kk;
    else kk=clg[k];
    
  clg[this->m]=nbcoef;
  // sort the new column (old line)
  for(int i=0;i<this->m;i++)  
    HeapSort(llg+clg[i],cl+clg[i],a+clg[i],clg[i+1]-clg[i]); 
  delete[] cl;
  delete[] lg;
  Exchange(this->n,this->m);       
  cl=llg;
  lg=clg;
  ConjArray(a,nbcoef);    
 }

template<class R>
 triplet<int,int,bool> BuildCombMat(std::map< pair<int,int>, R> & mij,const list<triplet<R,MatriceCreuse<R> *,bool> >  &lM,bool trans,int ii00,int jj00,bool cnj=false)
  {
    typedef typename list<triplet<R,MatriceCreuse<R> *,bool> >::const_iterator lconst_iterator;
    
    lconst_iterator begin=lM.begin();
    lconst_iterator end=lM.end();
    lconst_iterator i;
    
   // std::map< pair<int,int>, R> mij;
    
    int n=0,m=0;
    bool sym=true;
    for(i=begin;i!=end;i++++)
     {
	if(i->second) // M == 0 => zero matrix 
	{
	    MatriceCreuse<R> & M=*i->second;
	    bool transpose = i->third !=  trans;
	    ffassert( &M);
	    R coef=i->first;
	    if(verbosity>3)
		cout << "                BuildCombMat + " << coef << "*" << &M << " " << sym << "  t = " << transpose << " " <<  i->third << endl;
	    //  change to max FH dec 2007 to hard to satisfy
	   /* if (n==0)*/ { if(transpose) {m=max(m,M.n); n=max(n,M.m);} else{n=max(M.n,n); m=max(M.m,m);}}// Modif mars 2007 FH
	   /* else { if(transpose)  ffassert(n== M.m && m==M.n); else ffassert(n== M.n && m==M.m);}*/
	    sym = M.addMatTo(coef,mij,transpose,ii00,jj00,cnj) && sym;  
	}
     } 
    int nbcoef=mij.size();
    if(sym) nbcoef = (nbcoef+n)/2;

  // return new   MatriceMorse<R>(n,m,mij,sym);   
    return make_triplet(n,m,sym);
  }
  
template<class R>
  MatriceMorse<R> * BuildCombMat(const list<triplet<R,MatriceCreuse<R> *,bool> >  &lM,bool trans,int ii00,int jj00)
  {
   
    std::map< pair<int,int>, R> mij;
    triplet<int,int,bool> nmsym=BuildCombMat(mij,lM,trans,ii00,jj00);

   return new   MatriceMorse<R>(nmsym.first,nmsym.second,mij,nmsym.third);   
     
  }
template<class R>
bool MatriceMorse<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();
  int i,j,k;
  if (symetrique)
   {
     for ( i=0;i<this->n;i++)
       for ( k=lg[i];k<lg[i+1];k++)
         {
           j=cl[k];
           R cij =  coef* ( cnj ? conj(a[k]) : a[k]);
           if(norm(cij)>eps0)
           {
            mij[ij_mat(trans,ii00,jj00,i,j)] += cij ;
           if (i!=j)
             mij[ij_mat(trans,ii00,jj00,j,i)] += cij;
           }
         }
           
   }
  else
   {
     for ( i=0;i<this->n;i++)
       for ( k=lg[i];k<lg[i+1];k++)
         {
           j=cl[k];
           R cij =  coef* ( cnj ? conj(a[k]) : a[k]);

           if(norm(cij)>eps0)
           mij[ij_mat(trans,ii00,jj00,i,j)] += cij;
         }
   }

return symetrique;
}

 
template<class R> 
template<class K>
MatriceMorse<R>::MatriceMorse(int nn,int mm, std::map< pair<int,int>, K> & m, bool sym):
  MatriceCreuse<R>(nn,mm,0),
  nbcoef(m.size()),symetrique(sym),
  a(new R[nbcoef]),
  lg(new int[nn+1]),
  cl(new int[nbcoef]),     
  solver(0)
{
     int k=0;
     bool nosym=!sym;
     typename std::map< pair<int,int>, R>::iterator iter=m.begin(), mend=m.end();
     //  remarque lg est croissant Bug trouver par 
     for(int i=0;i<=nn;i++) lg[i]=0; 
     while(iter!=mend)
      { 
        int i=iter->first.first;
        int j=iter->first.second;
        K & aij=iter->second;
        assert( i < nn && j < mm);
        if(j<=i || nosym)
        {
         cl[k]=j;
         a[k]=aij;
         lg[i+1]=++k;
        }
        ++iter;
       }
    // lg est croissant  on bouche les trou   
   for(int i=1;i<=nn;i++) lg[i]=Max(lg[i-1],lg[i]); 
      
   ffassert(nbcoef==k);  
  }

template<class R>
 void  MatriceMorse<R>::resize(int nn,int mm) 
{
    int nc=0;   
    int *nlg=new int[nn+1],*ncl=0;
    int nm=min(nn,this->n);
    nc =0;
    nlg[0]=nc; 
    if (symetrique)
      {   if( nn != mm) AFAIRE("MatriceMorse<R>::resize symetric  n!=m");
	  for (int i=0;i<nm;i++)
	    {
	      for (int k=lg[i];k<lg[i+1];k++)
		{   int j=cl[k];
		    if( j<this->m && norm(a[k]))		    
			++nc;		   
		}
		nlg[i+1]=nc;
	    }
	  
      }
    else
      {
	  for (int i=0;i<nm;i++)
	    {
	      for (int k=lg[i];k<lg[i+1];k++)
		{
		    int j=cl[k];
		    if(i<this->n && j<this->m && norm(a[k]))
			++nc ;
		}
		nlg[i+1]=nc;
	    }
      }
    for(int i=nm+1;i<=nn;++i)
	nlg[i]=nc;
    ncl = new int[nc];
    R *na=new R[nc];
    nc=0;

    if (symetrique)
      {   if( nn != mm) AFAIRE("MatriceMorse<R>::resize symetric  n!=m");
	  for (int i=0;i<nm;i++)
	      for (int k=lg[i];k<lg[i+1];k++)
		{   int j=cl[k];
		    if( j<this->m && norm(a[k]))		    
		      {na[nc]=a[k];
		       ncl[nc++]=j;}
		}
	  
      }
    else
      {
	  for (int i=0;i<nm;i++)
	      for (int k=lg[i];k<lg[i+1];k++)
		{
		    int j=cl[k];
		    if( j<this->m && norm(a[k]))
		      {na[nc]=a[k];
		       ncl[nc++]=j;}
		}
      }
    
    delete [] cl;
    delete [] lg;
    delete [] a;
    cl=ncl;
    lg=nlg;
    a=na;
    this->n=nn;
    this->m=mm;
    this->N=nn;
    this->M=mm;
    this->nbcoef=nc;
 //   cout << nn << " " << mm << "  " <<  KN_<int>(lg,nn+1) << endl;
    
}
template<class RA>
 template<class RB,class RAB>
 void  MatriceMorse<RA>::prod(const MatriceMorse<RB> & B, MatriceMorse<RAB> & AB)
 {
   //  compute the s
  bool sym=this == & B &&symetrique;
  int *blg=B.lg;
  int *bcl=B.cl;
  ffassert(this->m==B.n); 
  bool delbl= B.symetrique;
  if (delbl)
    {
     int nn=B.n;
      blg = new int[nn+1];
     for (int i=0;i<B.n;i++)
         blg[i]=B.lg[i+1]-B.lg[i];
      blg[nn]=0;   
      
      for (int i=0;i<nn;i++)
        for (int k= B.lg[i];k<B.lg[i+1];k++)
          {  int j=B.cl[k];
              assert(j <= i);
             if (j!=i)  
               blg[j]++;
             }
             
      for (int i=1;i<=nn;i++)
       blg[i]+=blg[i-1];
      int nbnz = blg[nn];
      bcl= new int[nbnz];
      
      for (int i=0;i<B.n;i++)
        for (int k= B.lg[i];k<B.lg[i+1];k++)
          {  int j=B.cl[k];
             assert(j <= i);
             bcl[--blg[i] ]=j;
             if(i !=j)