lgmat.cpp

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

  ffassert(Uh);
  
  //  sparse_mat->pUh=pUh;
  //sparse_mat->pVh=pVh;
  sparse_mat->typemat=TypeSolveMat(TypeSolveMat::NONESQUARE); //  none square matrice (morse)
  sparse_mat->A.master(new MatriceMorse<R>(*Uh,*Vh,buildInterpolationMatrix,data));
  //  sparse_mat->A.master(new MatriceMorse<R>(*Uh,*Vh,buildInterpolationMatrix,data));
  }
  else 
  {  // new cas mars 2006
  pfes * pUh = GetAny< pfes * >((* mi->a)(stack));
  KN_<double>  xx = GetAny<  KN_<double>  >((* mi->b)(stack));
  KN_<double>  yy = GetAny<  KN_<double>  >((* mi->c)(stack));
  ffassert( xx.N() == yy.N()); 
  FESpace * Uh = **pUh;
  ffassert(Uh);
  
  //  sparse_mat->pUh=0;
  //  sparse_mat->pVh=0;
  sparse_mat->typemat=TypeSolveMat(TypeSolveMat::NONESQUARE); //  none square matrice (morse)
  sparse_mat->A.master(buildInterpolationMatrix1(*Uh,xx,yy,data));
  }
  return sparse_mat;
}

template<class RA,class RB,class RAB>
AnyType ProdMat(Stack stack,Expression emat,Expression prodmat)
{
  using namespace Fem2D;
  
  Matrice_Creuse<RAB> * sparse_mat =GetAny<Matrice_Creuse<RA>* >((*emat)(stack));
  const Matrix_Prod<RA,RB>  AB = GetAny<Matrix_Prod<RA,RB> >((*prodmat)(stack));
  //  sparse_mat->pUh=AB.A->pUh;
  //sparse_mat->pVh=AB.B->pVh;
  MatriceMorse<RA> *mA= AB.A->A->toMatriceMorse(AB.ta);
  MatriceMorse<RB> *mB= AB.B->A->toMatriceMorse(AB.tb);
  if( !mA && ! mB) ExecError(" Sorry error: in MatProd,  pb trans in MorseMat");
  if( mA->m != mB->n) {
    cerr << " -- Error dim ProdMat A*B : tA =" << AB.ta << " = tB " << AB.tb << endl;
    cerr << " --MatProd " << mA->n<< " "<< mA->m << " x " << mB->n<< " "<< mB->m <<  endl;
    ExecError(" Wrong mat dim in MatProd");
  }
  MatriceMorse<RAB> *mAB=new MatriceMorse<RA>();
  mA->prod(*mB,*mAB);
  
  sparse_mat->typemat=(mA->n == mB->m) ? TypeSolveMat(TypeSolveMat::GMRES) : TypeSolveMat(TypeSolveMat::NONESQUARE); //  none square matrice (morse)
  sparse_mat->A.master(mAB);
  delete mA;
  delete mB;
  return sparse_mat;
}

template<class R>
AnyType CombMat(Stack stack,Expression emat,Expression combMat)
{
  using namespace Fem2D;
  
  Matrice_Creuse<R> * sparse_mat =GetAny<Matrice_Creuse<R>* >((*emat)(stack));
  list<triplet<R,MatriceCreuse<R> *,bool> > *  lcB = GetAny<list<triplet<R,MatriceCreuse<R> *,bool> >*>((*combMat)(stack));
  //  sparse_mat->pUh=0;
  // sparse_mat->pVh=0; 
   MatriceCreuse<R> * AA=BuildCombMat<R>(*lcB,false,0,0);
  sparse_mat->A.master(AA);
  sparse_mat->typemat=(AA->n == AA->m) ? TypeSolveMat(TypeSolveMat::GMRES) : TypeSolveMat(TypeSolveMat::NONESQUARE); //  none square matrice (morse)
  delete lcB;
  return sparse_mat;
}


template<class R>
AnyType MatriceCreuse2map(Stack , const AnyType & mat)
{

  using namespace Fem2D;

  Matrice_Creuse<R> * sparse_mat =GetAny<Matrice_Creuse<R>* >(mat);
  ffassert(sparse_mat);
  int n=sparse_mat->N(),m=sparse_mat->M();
  map<pair<int,int>,R> *M=new map<pair<int,int>,R>;
  if (n >0 && m>0 && sparse_mat->A) 
    {
      sparse_mat->A->addMatTo(R(1.),*M);
      // hack 
      (*M)[make_pair(n-1,m-1)]+=R();
    }
  return M;
}


template<class R>
AnyType DiagMat(Stack stack,Expression emat,Expression edia)
{
  using namespace Fem2D;
  KN<R> * diag=GetAny<KN<R>* >((*edia)(stack));
  Matrice_Creuse<R> * sparse_mat =GetAny<Matrice_Creuse<R>* >((*emat)(stack));
  //  sparse_mat->pUh=0;
  // sparse_mat->pVh=0;
  sparse_mat->typemat=TypeSolveMat(TypeSolveMat::GC); //  none square matrice (morse)
  sparse_mat->A.master(new MatriceMorse<R>((int) diag->N(),(const R*) *diag));
  return sparse_mat;
}



template<class Rin,class Rout>
 struct  ChangeMatriceMorse {
 static  MatriceMorse<Rout> *f(MatriceMorse<Rin> *mr)
 {
    MatriceMorse<Rout>*  mrr=new MatriceMorse<Rout>(*mr);
    delete mr;
    return mrr;
 }
 };

template<class R> 
 struct  ChangeMatriceMorse<R,R> {
 static MatriceMorse<R>* f(MatriceMorse<R>* mr)
 {
   return mr;
 }
 };
template<class R,class RR>
AnyType CopyMat_tt(Stack stack,Expression emat,Expression eA,bool transp)
{
  using namespace Fem2D;
  Matrice_Creuse<R> * Mat;
 
  if(transp)
   {
    Matrice_Creuse_Transpose<R>  tMat=GetAny<Matrice_Creuse_Transpose<R> >((*eA)(stack));
    Mat=tMat; 
   }
  else   Mat =GetAny<Matrice_Creuse<R>*>((*eA)(stack));
  MatriceMorse<R> * mr=Mat->A->toMatriceMorse(transp,false);
  MatriceMorse<RR> * mrr = ChangeMatriceMorse<R,RR>::f(mr);
  
  Matrice_Creuse<RR> * sparse_mat =GetAny<Matrice_Creuse<RR>* >((*emat)(stack));
  //  sparse_mat->pUh=Mat->pUh;
  // sparse_mat->pVh=Mat->pUh;;
  sparse_mat->typemat=TypeSolveMat(TypeSolveMat::GC); //  none square matrice (morse)
  sparse_mat->A.master(mrr);
  return sparse_mat;
}

template<class R,class RR>
AnyType CopyTrans(Stack stack,Expression emat,Expression eA)
{
 return CopyMat_tt<R,RR>(stack,emat,eA,true);
}
template<class R,class RR>
AnyType CopyMat(Stack stack,Expression emat,Expression eA)
{
 return CopyMat_tt<R,RR>(stack,emat,eA,false);
}

template<class R>
AnyType MatFull2Sparse(Stack stack,Expression emat,Expression eA)
{
  KNM<R> * A=GetAny<KNM<R>* >((*eA)(stack));
  Matrice_Creuse<R> * sparse_mat =GetAny<Matrice_Creuse<R>* >((*emat)(stack));
  //  sparse_mat->pUh=0;
  // sparse_mat->pVh=0;
  sparse_mat->typemat=TypeSolveMat(TypeSolveMat::GMRES); //  none square matrice (morse)
  sparse_mat->A.master(new MatriceMorse<R>((KNM_<R> &)*A,0.0));
  
 return sparse_mat;
}

template<class R>
AnyType MatMap2Sparse(Stack stack,Expression emat,Expression eA)
{
   map< pair<int,int>, R> * A=GetAny< map< pair<int,int>, R> * >((*eA)(stack));
   int n=0,m=0;
   // hack:  the last element must exist in the map  to set matrix size 

   typename map< pair<int,int>, R>::const_iterator last= --A->end(); // le last element
   
   if( last != A->end() )
   
        { 
                n = last->first.first+1; 
                m=last->first.second+1;
        } 
        
  Matrice_Creuse<R> * sparse_mat =GetAny<Matrice_Creuse<R>* >((*emat)(stack));
  //  sparse_mat->pUh=0;
  // sparse_mat->pVh=0;
  sparse_mat->typemat=TypeSolveMat(TypeSolveMat::GMRES); //  none square matrice (morse)  
  sparse_mat->A.master(new MatriceMorse<R>(n,m,*A,false));
  delete A; 
 return sparse_mat;
}

template<class RR,class AA=RR,class BB=AA> 
struct Op2_pair: public binary_function<AA,BB,RR> { 
  static RR f(const AA & a,const BB & b)  
  { return RR( a, b);} 
};


template<class R>
long get_mat_n(Matrice_Creuse<R> * p)
 { ffassert(p ) ;  return p->A ?p->A->n: 0  ;}

template<class R>
long get_mat_m(Matrice_Creuse<R> * p)
 { ffassert(p ) ;  return p->A ?p->A->m: 0  ;}

template<class R>
long get_mat_nbcoef(Matrice_Creuse<R> * p)
 { ffassert(p ) ;  return p->A ?p->A->NbCoef(): 0  ;}

template<class R>
pair<long,long> get_NM(const list<triplet<R,MatriceCreuse<R> *,bool> > & lM)
{
      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;
    
    long n=0,m=0;
    for(i=begin;i!=end;i++++)
     {
       ffassert(i->second);
       MatriceCreuse<R> & M=*i->second;
       bool t=i->third;
       int nn= M.n,mm=M.m;
       if (t) swap(nn,mm);
       if ( n==0) n =  nn;
       if ( m==0) m = mm;
       if (n != 0) ffassert(nn == n);
       if (m != 0) ffassert(mm == m);
      }
   return make_pair(n,m);    
}

template<class R>
long get_diag(Matrice_Creuse<R> * p, KN<R> * x)
 { ffassert(p && x ) ;  return p->A ?p->A->getdiag(*x): 0  ;}
template<class R>
long set_diag(Matrice_Creuse<R> * p, KN<R> * x)
 { ffassert(p && x ) ;  return p->A ?p->A->setdiag(*x): 0  ;}

 
template<class R>  
R * get_elementp2mc(Matrice_Creuse<R> * const  & ac,const long & b,const long & c){ 
   MatriceCreuse<R> * a= ac ? ac->A:0 ;
  if(  !a || a->n <= b || c<0 || a->m <= c  ) 
   { cerr << " Out of bound  0 <=" << b << " < "  << a->n << ",  0 <= " << c << " < "  << a->m
           << " Matrix type = " << typeid(ac).name() << endl;
     cerr << ac << " " << a << endl;
     ExecError("Out of bound in operator Matrice_Creuse<R> (,)");}
   R *  p =a->pij(b,c);
   if( !p) { if(verbosity) cerr << "Error: the coef a(" << b << ","   << c << ")  do'nt exist in sparse matrix "
           << " Matrix  type = " << typeid(ac).name() << endl;
       ExecError("Use of unexisting coef in sparse matrix operator a(i,j) ");}
    return  p;}


template<class RR,class AA=RR,class BB=AA> 
struct Op2_mulAv: public binary_function<AA,BB,RR> { 
  static RR f(const AA & a,const BB & b)  
  { return (*a->A * *b );} 
};

template<class RR,class AA=RR,class BB=AA> 
struct Op2_mulvirtAv: public binary_function<AA,BB,RR> { 
  static RR f(const AA & a,const BB & b)  
  { return RR( (*a).A, b );} 
};
 
 
template<class K>
class OneBinaryOperatorA_inv : public OneOperator { public:  
  OneBinaryOperatorA_inv() : OneOperator(atype<Matrice_Creuse_inv<K> >(),atype<Matrice_Creuse<K> *>(),atype<long>()) {}
    E_F0 * code(const basicAC_F0 & args) const 
     { Expression p=args[1];
       if ( ! p->EvaluableWithOutStack() ) 
        { 
          bool bb=p->EvaluableWithOutStack();
          cout << bb << " " <<  * p <<  endl;
          CompileError(" A^p, The p must be a constant == -1, sorry");}
       long pv = GetAny<long>((*p)(0));
        if (pv !=-1)   
         { char buf[100];
           sprintf(buf," A^%ld, The pow must be  == -1, sorry",pv);
           CompileError(buf);}     
       return  new E_F_F0<Matrice_Creuse_inv<K>,Matrice_Creuse<K> *>(Build<Matrice_Creuse_inv<K>,Matrice_Creuse<K> *>,t[0]->CastTo(args[0])); 
    }
};

template<class K>
class Psor :  public E_F0 { public:  
 
   typedef double  Result;
   Expression mat;
   Expression xx,gmn,gmx,oomega;  
   Psor(const basicAC_F0 & args) 
    {   
      args.SetNameParam();
      mat=to<Matrice_Creuse<K> *>(args[0]); 
      gmn=to<KN<K>*>(args[1]); 
      gmx=to<KN<K>*>(args[2]); 
      xx=to<KN<K>*>(args[3]); 
      oomega=to<double>(args[4]); 
      
   }   
    static ArrayOfaType  typeargs() { 
      return  ArrayOfaType( atype<double>(),
                            atype<Matrice_Creuse<K> *>(),
                            atype<KN<K>*>(),
                            atype<KN<K>*>(),
                            atype<KN<K>*>(),
                            atype<double>(),false);}
                            
    static  E_F0 * f(const basicAC_F0 & args){ return new Psor(args);} 
    
    AnyType operator()(Stack s) const {
      Matrice_Creuse<K>* A= GetAny<Matrice_Creuse<K>* >( (*mat)(s) );
      KN<K>* gmin = GetAny<KN<K>* >( (*gmn)(s) );
      KN<K>* gmax = GetAny<KN<K>* >( (*gmx)(s) );
      KN<K>* x = GetAny<KN<K>* >( (*xx)(s) );
      double omega = GetAny<double>((*oomega)(s));
      return A->A->psor(*gmin,*gmax,*x,omega);
    }
  
};
template <class R>
 struct TheDiagMat {
  Matrice_Creuse<R> * A; 
  TheDiagMat(Matrice_Creuse<R> * AA) :A(AA) {ffassert(A);}
  void   get_mat_daig( KN_<R> & x) { ffassert(A && A->A && x.N() == A->A->n  && A->A->n == A->A->m );
     A->A->getdiag(x);}
  void   set_mat_daig(const  KN_<R> & x) { ffassert(A && A->A && x.N() == A->A->n  && A->A->n == A->A->m );
     A->A->setdiag(x);}
 };
 
 template <class R>
 struct TheCoefMat {
  Matrice_Creuse<R> * A; 
  TheCoefMat(Matrice_Creuse<R> * AA) :A(AA) {ffassert(A);}
  void   get_mat_coef( KN_<R> & x) { ffassert(A && A->A && x.N() == A->A->NbCoef()  );
     A->A->getcoef(x);}
  void   set_mat_coef(const  KN_<R> & x) { ffassert(A && A->A && x.N() == A->A->NbCoef() );
     A->A->setcoef(x);}
 };
 
template<class R>
TheDiagMat<R> thediag(Matrice_Creuse<R> * p)
 {  return  TheDiagMat<R>(p);}

template<class R>
TheCoefMat<R> thecoef(Matrice_Creuse<R> * p)
 {  return  TheCoefMat<R>(p);}
 
template<class R>
TheDiagMat<R> set_mat_daig(TheDiagMat<R> dm,KN<R> * x)
{
  dm.set_mat_daig(*x);
  return dm;
}
template<class R>
KN<R> * get_mat_daig(KN<R> * x,TheDiagMat<R> dm)
{
  dm.get_mat_daig(*x);
  return x;
}

template<class R>
TheCoefMat<R> set_mat_coef(TheCoefMat<R> dm,KN<R> * x)
{
  dm.set_mat_coef(*x);
  return dm;
}
template<class R>
KN<R> * get_mat_coef(KN<R> * x,TheCoefMat<R> dm)
{
  dm.get_mat_coef(*x);
  return x;
}

template<class R> 
bool IsRawMat(const basicAC_F0 & args) 
{   
    
    const E_Array & ee= *dynamic_cast<const E_Array*>((Expression) args[1]);
    if (!&ee) return 0;
    
    int N=ee.size();
    if (N==1)
    {
	C_F0 c0(ee[0]);
	return 
	    atype<KN_<R> >()->CastingFrom(ee[0].left());
	    
    }
    else if (N==3)
    {
	C_F0 c0(ee[0]),c1(ee[1]),c2(ee[2]);
	return 
	    atype<KN_<long> >()->CastingFrom(ee[0].left())
	    && 	    atype<KN_<long> >()->CastingFrom(ee[1].left())
	    &&      atype<KN_<R> >()->CastingFrom(ee[2].left());
	    
    }
    return 0;
}


template<typename R>
class RawMatrix :  public E_F0 { public: 
    typedef Matrice_Creuse<R> * Result;
    Expression emat; 
    Expression coef,col,lig;
    RawMatrix(const basicAC_F0 & args) ;
    static ArrayOfaType  typeargs() { return  ArrayOfaType(atype<Matrice_Creuse<R>*>(),atype<E_Array>());}
    AnyType operator()(Stack s) const ;    
};

 template<typename R>
 class BlockMatrix :  public E_F0 { public: