matricecreuse_tpl.hpp

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

               bcl[--blg[j]]=i;
          }
    }
   
   set<pair<int,int> > sij;
   double eps0=numeric_limits<double>::min();

     for (int i=0;i<this->n;i++)
       for (int k=lg[i];k<lg[i+1];k++)
         {    
           int j=cl[k];
           if(norm(a[k])<eps0) continue;
           int ii[2],jj[2];
           ii[0]=i;ii[1]=j;
           jj[0]=j;jj[1]=i;
           int kk=1;
           if(symetrique && i != j) kk=2;
           for (int ll=0;ll<kk;ll++)
            {
                int i=ii[ll];
                int j=jj[ll];
		if(j>=B.n) continue; // in case of not equal size A.m != B.n 
                for (int kkb=blg[j];kkb<blg[j+1];kkb++)
                  { 
                   int kz= bcl[kkb];
                   RB bjk;
                   if (B.symetrique && kz > j)
                     bjk=B(kz,j);
                   else
                      bjk=B(j,kz);
                   if( norm(bjk)>eps0 && (!sym || kz<=i))
                     sij.insert(make_pair(i,kz));
                  }
            }
           
         }
    int nn=this->n;
    int mm=B.m;
    int * llg=new int[nn+1];
    int * lcl=new int[sij.size()];  
    RAB * aa = new RAB[sij.size()];
    for(int i=0;i<=nn;i++)
        llg[i]=0;
        
    for (set<pair<int,int> >::iterator iter=sij.begin();iter!=sij.end();++iter)
      { 
        int i=iter->first;
	// int j=iter->second;
        llg[i]++;
       }
     for (int i=1;i<=nn;i++)
       llg[i]+=llg[i-1];
     ffassert(llg[this->n]==(long) sij.size());
     for (set<pair<int,int> >::iterator iter=sij.begin();iter!=sij.end();++iter)
      { 
        int i=iter->first;
        int j=iter->second;
       // cout << i << " , " << j << endl;
        lcl[--llg[i]]=j;
       }
     for(int i=0;i<nn;i++)  
       HeapSort(lcl+llg[i],llg[i+1]-llg[i]); 
       
     AB.n=nn;
     AB.m=mm;
     AB.N=nn;  // add missing jan 2008 FH
     AB.M=mm;  // add missing jan 2008 FH

     AB.lg=llg;
     AB.cl=lcl;
     AB.a=aa;        
     AB.nbcoef=sij.size();
     AB.symetrique=sym;
     AB.dummy=false;
     AB = RAB();
     for (int i=0;i<this->n;i++)
       for (int k=lg[i];k<lg[i+1];k++)
         {    
           int j=cl[k];
           RAB aij = a[k];
           if(norm(aij) <eps0 ) continue;
           int ii[2],jj[2];
           ii[0]=i;ii[1]=j;
           jj[0]=j;jj[1]=i;
           int kk=1;
           if(symetrique && i != j) kk=2;
           for (int ll=0;ll<kk;ll++)
            {
                int i=ii[ll];
                int j=jj[ll];
		if(j>=B.n) continue; // in case of not equal size A.m != B.n 
                for (int kb=blg[j];kb<blg[j+1];kb++)
                  { 
                   int k= bcl[kb];
                   RB bjk;
                   if (B.symetrique && k > j)
                     bjk=B(k,j);
                   else
                      bjk=B(j,k);
                //   cout << i << "," << "," << j << "," << k << " " << aij << " " << bjk << endl;
                   if( norm( bjk)> eps0  && (!sym || k<=i))
                       AB(i,k) += aij*bjk;
                  }
            }
           
         }

    if (delbl) {
      delete [] blg;
      delete [] bcl;
    }
     
     
 }

template<class R>
  void  MatriceMorse<R>::addMatMul(const KN_<R> &  x, KN_<R> & ax) const   
{
  int i,j,k;
  if( ! (this->n==ax.N() && this->m==x.N()))
    {cerr << " Err MatriceMorse<R>:  ax += A x" <<endl;
      cerr << " A.n " << this->n<< " !=  "<< ax.N() << " ax.n \n";
      cerr << " A.m " << this->m<< " != " <<x.N() << " x.n \n" ;
      ffassert(0); 
      abort();
    }
  if (symetrique)
   {
     for (i=0;i<this->n;i++)
       for (k=lg[i];k<lg[i+1];k++)
         {
           j=cl[k];
           ax[i] += a[k]*x[j];
           if (i!=j)
             ax[j] += a[k]*x[i];
         }
           
   }
  else
   {
     for (i=0;i<this->n;i++)
       for (k=lg[i];k<lg[i+1];k++)
         {
           j=cl[k];
           ax[i] += a[k]*x[j];
         }
   }
}

template<class R>
  void  MatriceMorse<R>::addMatTransMul(const KN_<R> &  x, KN_<R> & ax) const   
{
  int i,j,k;
  ffassert(this->m==ax.N());
  ffassert(this->n==x.N());  
  if (symetrique)
   {
     for (i=0;i<this->n;i++)
       for (k=lg[i];k<lg[i+1];k++)
         {
           j=cl[k];
           ax[j] += a[k]*x[i];
           if (i!=j)
             ax[i] += a[k]*x[j];
         }
           
   }
  else
   {
     for (i=0;i<this->n;i++)
       for (k=lg[i];k<lg[i+1];k++)
         {
           j=cl[k];
           ax[j] += a[k]*x[i];
         }
   }
}


template<class R>
MatriceMorse<R>  & MatriceMorse<R>::operator +=(MatriceElementaire<R> & me) {
  int il,jl,i,j;
  int * mi=me.ni, *mj=me.nj;
  if ((this->n==0) && (this->m==0))
   {
   
    //    if(verbosity>3)
    cout << " -- Morse Matrice is empt: let's build it" << endl;
    ffassert(0); 
    /*
    this->n=me.Uh.NbOfDF;
    this->m=me.Vh.NbOfDF;
    switch (me.mtype) {
    case MatriceElementaire<R>::Full : 
      Build(me.Uh,me.Vh,false);    
      break;
    case MatriceElementaire<R>::Symmetric :     
      Build(me.Uh,me.Vh,true);    
      break;
     default:
      cerr << "Big bug type MatriceElementaire is unknown" << (int) me.mtype << endl;
      throw(ErrorExec("exit",1));
      break; }     
    */
   }
  R * al = me.a; 
  R * aij;
  switch (me.mtype) { // modif FH overfloat in array mi and mj => trap on win32
  case MatriceElementaire<R>::Full : ffassert(!symetrique);
    for (il=0; il<me.n; ++il)  { i=mi[il]; 
      for ( jl=0; jl< me.m ; ++jl,++al)  {j=mj[jl];
        aij = pij(i,j);
        throwassert(aij);
	*aij += *al;}}
    break;
     
  case MatriceElementaire<R>::Symmetric : ffassert(symetrique);   
    for (il=0; il<me.n; ++il) {  i=mi[il] ;
      for (jl=0;jl< il+1 ; ++jl) { j=mj[jl];
	 aij =    (j<i) ? pij(i,j) : pij(j,i);
         throwassert(aij);
         *aij += *al++;}}
    break;
  default:
    cerr << "Big bug type MatriceElementaire unknown" << (int) me.mtype << endl;
    exit(1);
    break; 
  }      
  return *this;
} 

template<class R>
  void MatriceMorse<R>::Solve(KN_<R> &x,const KN_<R> &b) const{
    if (solver)    
      solver->Solver(*this,x,b);
    else
  {  cerr << "No Solver defined  for this Morse matrix " << endl;
    throw(ErrorExec("exit",1));}
  }


template<class R>
double MatriceMorse<R>::psor(KN_<R> & x,const  KN_<R> & gmin,const  KN_<R> & gmax , double omega) 
{
  double err=0;
  int n=this->n;
  ffassert(n==this->m);
  ffassert(n==x.N());
  ffassert(n==gmin.N());
  ffassert(n==gmax.N());
  if (symetrique)
   {
     ErrorExec("Error:sorry psor just for no symmetric Morse matrices",1);
   }
  else
   {
     for (int i=0;i<this->n;i++)
      {
       R xnew =x[i];
       R aii=R();
       for (int k=lg[i];k<lg[i+1];k++)
         {
           int j=cl[k];
           if(j!= i) 
             xnew -= a[k]*x[j];
            else aii=a[k];
         }
        if(aii != R())
           xnew /= aii;
         else ErrorExec("Error: psor diagonal coef = 0 ",1);
        R dx  = (xnew - x[i])*omega ;
        R xi = RNM::Min(RNM::Max(x[i]+dx,gmin[i]),gmax[i]);
        dx = x[i]- xi;
        err = Max(err, norm(dx));
        x[i] = xi;
        }
   }  return sqrt(err);
  
}

template<class R>
double MatriceProfile<R>::psor(KN_<R> & x,const  KN_<R> & gmin,const  KN_<R> & gmax , double omega) 
{
  double rr=0;
  ErrorExec("Error:sorry psor just for no symmetric Morse matrices (will do in future FH??? )",2);
  return rr;
  
}

template<class R>
void MatriceProfile<R>::setdiag(const KN_<R> & x) 
{
  ffassert(D);
 ffassert( this->n == x.N());
  KN_<R> d(D,this->n) ;
  d=x;
}
template<class R>
void MatriceProfile<R>::getdiag(KN_<R> & x) const 
{
  ffassert(D);
  ffassert( this->n == x.N());
  KN_<R> d(D,this->n) ;
  x=d;  
}
template<class R>
void MatriceMorse<R>::setdiag(const KN_<R> & x) 
{
 ffassert( this->n == this->m&& this->n == x.N());
 for (int i=0;i<this->n;++i)
    {
      R * p= pij(i,i);
      if(p)     *p = x[i];
      else ffassert( norm(x[i]) < 1e-30);}
}
template<class R>
void MatriceMorse<R>::getdiag(KN_<R> & x) const 
{
 ffassert( this->n == this->m && this->n == x.N());
 for (int i=0;i<this->n;++i)
    {
      R * p= pij(i,i);
      x[i]=  p ?  *p : R() ;
    }
  
}
template<class R>
R MatriceMorse<R>::pscal(const KN_<R> & x,const KN_<R> & y)
{ // (x, Ay)
  R sum=R();
  int i,j,k;
  ffassert(this->n==x.N());
  ffassert(this->m==y.N());  
  if (symetrique)
   {
     for (i=0;i<this->n;i++)
       for (k=lg[i];k<lg[i+1];k++)
         {
           j=cl[k];
           sum += a[k]*x[i]*y[j];
           if (i!=j)
             sum += a[k]*x[j]*y[i];
         }
           
   }
  else
   {
     for (i=0;i<this->n;i++)
       for (k=lg[i];k<lg[i+1];k++)
         {
           j=cl[k];
           sum += a[k]*x[i]*y[j];
         }
   }
  return sum;
}
template<class R>
R MatriceProfile<R>::pscal(const KN_<R> & x,const KN_<R> & y)
{
 if (y.n != this->n || x.n != this->n ) ERREUR(MatriceProfile pscal(xa,x) ," longueur incompatible c (out)") ;
 int i,j,k,kf;
 R sum = R();
 ffassert(this->n == this->m);
 if (D) 
   for (i=0;i<this->n;i++) 
     sum += D[i]*x[i]*y[i];
 else
   for (i=0;i<this->n;i++) // no dia => identyty dai
     sum +=x[i]*y[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++  )
       sum += L[k]*x[i]*y[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++  )
       sum += U[k]*x[i]*y[j],throwassert(i>=0 && i <this->n && j >=0 && j < this->m &&  k>=0 && k < pU[this->n]);
 
 return sum;
}

template<class R>
void MatriceMorse<R>::getcoef(KN_<R> & x) const 
{
 ffassert(x.N()==this->nbcoef);
 x = KN_<R>(this->a,nbcoef);  
}
template<class R>
void MatriceMorse<R>::setcoef(const KN_<R> & x)  
{
 ffassert(x.N()==nbcoef);
  KN_<R>(this->a,nbcoef) = x;
}
template<class R>
int MatriceMorse<R>::NbCoef() const  
{
  return this->nbcoef;
}

template<class R>
void MatriceProfile<R>::getcoef(KN_<R> & x) const 
{
 ffassert(x.N()==this->NbCoef());
 int k=0,kk;
 if (D)
  {  kk=this->n;
     x(SubArray(kk,k))  = KN_<R>(D,kk);
     k += kk; }
 if (L)
  {  kk= pL[this->n];
     x(SubArray(kk,k))  = KN_<R>(L,kk);
     k += kk; }
  if (U && (U != L)) 
  {  kk=  pU[this->n];
     x(SubArray(kk,k))  = KN_<R>(U,kk);
     k += kk; }
   
}
template<class R>
void MatriceProfile<R>::setcoef(const KN_<R> & x)  
{
 ffassert(x.N()==this->NbCoef());
   int k=0,kk;
 if (D)
  {  kk=this->n;
     KN_<R>(D,kk)=x(SubArray(kk,k))   ;
     k += kk; }
 if (L)
  {  kk= pL[this->n];
     KN_<R>(L,kk)=x(SubArray(kk,k))   ;
     k += kk; }
  if (U && (U != L)) 
  {  kk=  pU[this->n];
     KN_<R>(U,kk)=x(SubArray(kk,k)) ;
     k += kk; }

}
template<class R>
int MatriceProfile<R>::NbCoef() const  
{
  int s=0;
  if (D) s += this->n;
  if (L) s += pL[this->n];
  if (U && (U != L)) s += pU[this->n];
  return s;
}
#endif