lgmat.cpp

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

		t_Mij[i][j]=6;
	    }
	    else if ( atype<R >()->CastingFrom(rij) )
	    {  		  // frev 2007 
		e_Mij[i][j]=to<R>(c_Mij);
		t_Mij[i][j]=7; //  just un scalaire 
	    }
	    
	    else {  
		
		CompileError(" Block matrix ,  bad type in block matrix");
	    }
	    /*            else if   ( atype<map< pair<int,int>, R> * >()->CastingFrom(rij) ) 
	    {
		e_Mij[i][j]= to<map< pair<int,int>, R> *>(C_F0(eij,rij)).LeftValue();
		t_Mij[i][j]=10;
	    }*/
	}
	
    }
}

template<typename RR>
class  SetRawMatformMat : public OneOperator { 
public:
    typedef Matrice_Creuse<RR> *  A; // Warning  B type of  2 parameter 
    typedef Matrice_Creuse<RR> *  R;
    typedef E_Array B; //   A type of 1 parameter
    
    class CODE : public  E_F0 { public:
	Expression Mat;
	Expression lig;
	Expression col;
	Expression coef;
	bool mi;    
	    CODE(Expression a,const E_Array & tt)  
		: Mat(a),
		 mi(tt.MeshIndependent())
	    {

		    assert(&tt);
		    if(tt.size()!=3) 
			CompileError("Set raw matrix:  [ lg,col, a] = A (size !=3) ");
		    if (    aatypeknlongp->CastingFrom(tt[0].left() ) //// for  compilation error with g++ 3.2.2 (4 times)
			&&  aatypeknlongp->CastingFrom(tt[1].left() )
			&&  atype<KN<RR>* >()->CastingFrom(tt[2].left() ) )
			    {
			      lig = aatypeknlongp->CastTo(tt[0]);
			      col = aatypeknlongp->CastTo(tt[1]);
			      coef = atype<KN<RR>* >()->CastTo(tt[2]);
			    }      
			    else 
				CompileError(" we are waiting for [ lg,col,a] = A");
    }
	    
	    AnyType operator()(Stack stack)  const 
	    {
		A  a=GetAny<A>((*Mat)(stack));
		KN<long> *lg,*cl;
		KN<RR> *cc;
		lg = GetAny<KN<long>*>((*lig)(stack));
		cl = GetAny<KN<long>*>((*col)(stack));
		cc = GetAny<KN<RR>*>((*coef)(stack));
		int n=a->N(),m=a->M();
		map<pair<int,int>,RR> *M=new map<pair<int,int>,RR>;
		if (n >0 && m>0 && a->A) 
		{
		    a->A->addMatTo(RR(1.),*M);
		    // hack 
		    (*M)[make_pair(n-1,m-1)]+=RR();
		}
		int kk = M->size();
		lg->resize(kk);
		cc->resize(kk);
		cl->resize(kk);
		int k=0;
		typename map<pair<int,int>,RR>::const_iterator i;
		//if (!a->sym)
		 for (i=M->begin(); i != M->end();++i,++k)
		  {  
		    (*lg)[k]= i->first.first;
		    (*cl)[k]= i->first.second;
		    (*cc)[k]= i->second;
		  }
		    
		delete M;
		return SetAny<R>(a);
	    } 
	    bool MeshIndependent() const     {return  mi;} // 
	    ~CODE() {}
	    operator aType () const { return atype<R>();}    
    }; // end sub class CODE
    
    
public: // warning hack  A and B 
	E_F0 * code(const basicAC_F0 & args) const 
    { return  new CODE(t[1]->CastTo(args[1]),*dynamic_cast<const E_Array*>( t[0]->CastTo(args[0]).RightValue()));} 
    SetRawMatformMat():   OneOperator(atype<R>(),atype<B>(),atype<A>())  {} // warning with A and B 
    
};


template<typename R>  AnyType RawMatrix<R>::operator()(Stack stack) const
{
    MatriceMorse<R> * amorse =0; 
    KN_<R> cc(GetAny< KN_<R>  >((*coef)(stack)));
    int k= cc.N();
    int n= k;
    int m=n;
    map< pair<int,int>, R> Aij;
    bool sym=false;
    if( lig && col)
    {
	KN_<long> lg(GetAny< KN_<long>  >((*lig)(stack)));
	KN_<long> cl=(GetAny< KN_<long>  >((*col)(stack)));
	n = lg.max()+1;
	m = cl.max()+1;
	ffassert( lg.N()==k && cl.N()==k && lg.min()>=0 && lg.max()>=0);
	sym=false;
	for(int i=0;i<k;++i)
	    Aij[make_pair<int,int>(lg[i],cl[i])]+=cc[i];	
    }
    else
    {
	sym=true;
	for(int i=0;i<n;++i)
	    Aij[make_pair(i,i)]=cc[i];
    }
    //cout <<  " nxm  =" <<n<< "x" << m <<endl; 
    amorse=  new MatriceMorse<R>(n,m,Aij,sym); 
    if(verbosity)
	cout << " -- Raw Matrix    nxm  =" <<n<< "x" << m << " nb  none zero coef. " << amorse->nbcoef << endl;
    
    Matrice_Creuse<R> * sparse_mat =GetAny<Matrice_Creuse<R>* >((*emat)(stack));       
    // sparse_mat->pUh=0;
    // sparse_mat->pVh=0; 
    sparse_mat->A.master(amorse);
    sparse_mat->typemat=(amorse->n == amorse->m) ? TypeSolveMat(TypeSolveMat::GMRES) : TypeSolveMat(TypeSolveMat::NONESQUARE); //  none square matrice (morse)
    
    if(verbosity>3) { cout << "  End Raw Matrix : " << endl;}
    
    return sparse_mat;   
}
template<typename R>  AnyType BlockMatrix<R>::operator()(Stack s) const
{
  typedef list<triplet<R,MatriceCreuse<R> *,bool> > * L;
   KNM<L> Bij(N,M);
   KNM<KNM_<R> * > Fij(N,M); 
   KNM<bool> cnjij(N,M); 
   KNM<R> Rij(N,M); //  to sto
  
   cnjij = false; 
   KN<long> Oi(N+1), Oj(M+1);
   if(verbosity>3) { cout << " Build Block Matrix : " << N << " x " << M << endl;}
   Bij = (L) 0;
   Oi = (long) 0;
   Oj = (long)0;
  for (int i=0;i<N;++i)
   for (int j=0;j<M;++j)
    {
      Fij(i,j)=0;
      Expression eij = e_Mij[i][j];
      int tij=t_Mij[i][j];
      if (eij) 
      {
        cnjij(i,j) = tij%2 == 0; 
        AnyType e=(*eij)(s);
        if (tij==1) Bij(i,j) = to( GetAny< Matrice_Creuse<R>* >( e)) ;
        else if  (tij==2) Bij(i,j) = to( GetAny<Matrice_Creuse_Transpose<R> >(e));
        else if (tij==3)  { KN_<R> x=GetAny< KN_<R>  >( e);  Fij(i,j) = new KNM_<R>(x,x.N(),1);}
        else if (tij==4)  { KN_<R> x=GetAny< Transpose< KN_<R> >   >( e).t ;  Fij(i,j) = new KNM_<R>(x,1,x.N());}
        else if (tij==5)  { KNM<R> * m= GetAny< KNM<R>*  >( e);  Fij(i,j) = new KNM_<R>(*m);}
        else if (tij==6)  { KNM<R> * m= GetAny< Transpose< KNM<R>* >  >( e).t;  Fij(i,j) = new KNM_<R>(m->t()); }
	else if (tij==7)   { Rij(i,j)=GetAny< R  >( e);  Fij(i,j) = new KNM_<R>(&(Rij(i,j)),1,1);}
        
   //     else if  (tij==3) {}
        else {
         cout << " Bug " << tij << endl;
         ExecError(" Type sub matrix block unknown ");
        }
      }
     }
     //  compute size of matrix
     int err=0;
    for (int i=0;i<N;++i)
     for (int j=0;j<M;++j) 
       {
        pair<long,long> nm(0,0);
        
       if (Bij(i,j)) 
         nm = get_NM( *Bij(i,j));
       else if(Fij(i,j)) 
         nm = make_pair<long,long>(Fij(i,j)->N(), Fij(i,j)->M());
          
        if (( nm.first || nm.second)  && verbosity>3)
          cout << " Block [ " << i << "," << j << " ]      =     " << nm.first << " x " << nm.second << " cnj = " << cnjij(i,j) << endl;
        if (nm.first)
          {
          if ( Oi(i+1) ==0 )  Oi(i+1)=nm.first;
          else  if(Oi(i+1) != nm.first)
            { 
                 err++;
                 cerr <<"Error Block Matrix,  size sub matrix" << i << ","<< j << " n (old) "  << Oi(i+1) 
                       << " n (new) " << nm.first << endl;

            }
          }
          if(nm.second) 
          {  
          if   ( Oj(j+1) ==0) Oj(j+1)=nm.second;
          else   if(Oj(j+1) != nm.second) 
            { 
              cerr <<"Error Block Matrix,  size sub matrix" << i << ","<< j << " m (old) "  << Oj(j+1) 
                   << " m (new) " << nm.second << endl;
              err++;}
          }
        }
    if (err)    ExecError("Error Block Matrix,  size sub matrix");
//  cout << "Oi = " <<  Oi << endl;
//  cout << "Oj = " <<  Oj << endl;

    for (int i=0;i<N;++i)
      Oi(i+1) += Oi(i);
    for (int j=0;j<M;++j) // correct 10/01/2007 FH 
      Oj(j+1) += Oi(j);
  long n=Oi(N),m=Oj(M);
  if(verbosity>3)
   {
     cout << "     Oi = " <<  Oi << endl;
     cout << "     Oj = " <<  Oj << endl;
  }
  MatriceMorse<R> * amorse =0; 
{
   map< pair<int,int>, R> Aij;
    for (int i=0;i<N;++i)
     for (int j=0;j<M;++j) 
       if (Bij(i,j)) 
         {
           if(verbosity>3)
             cout << "  Add  Block S " << i << "," << j << " =  at " << Oi(i) << " x " << Oj(j) << " conj = " << cnjij(i,j) << endl;
           BuildCombMat(Aij,*Bij(i,j),false,Oi(i),Oj(j),cnjij(i,j));
         }
       else if (Fij(i,j))
        {
           if(verbosity>3)
             cout << "  Add  Block F " << i << "," << j << " =  at " << Oi(i) << " x " << Oj(j) << endl;
           BuildCombMat(Aij,*Fij(i,j),Oi(i),Oj(j),R(1.),cnjij(i,j));// BuildCombMat
        }
        
           
  amorse=  new MatriceMorse<R>(n,m,Aij,false); 
  }
  if(verbosity)
     cout << " -- Block Matrix NxM = " << N << "x" << M << "    nxm  =" <<n<< "x" << m << " nb  none zero coef. " << amorse->nbcoef << endl;
  
  Matrice_Creuse<R> * sparse_mat =GetAny<Matrice_Creuse<R>* >((*emat)(s));       
  //sparse_mat->pUh=0;
  // sparse_mat->pVh=0; 
  sparse_mat->A.master(amorse);
  sparse_mat->typemat=(amorse->n == amorse->m) ? TypeSolveMat(TypeSolveMat::GMRES) : TypeSolveMat(TypeSolveMat::NONESQUARE); //  none square matrice (morse)
                
     
  // cleanning    
  for (int i=0;i<N;++i)
   for (int j=0;j<M;++j)
    if(Bij(i,j)) delete Bij(i,j);
    else if(Fij(i,j))  delete Fij(i,j);  
   if(verbosity>3) { cout << "  End Build Blok Matrix : " << endl;}
   
 return sparse_mat;  

}

template<class R>
class minusMat { public:
    list<triplet<R,MatriceCreuse<R> *,bool> >  *l;
    minusMat(list<triplet<R,MatriceCreuse<R> *,bool> > *ll):
	l(new list<triplet<R,MatriceCreuse<R> *,bool> >(*ll) )
      {
	    typedef typename list<triplet<R,MatriceCreuse<R> *,bool> >::iterator lci;
	    for (lci i= l->begin();i !=l->end();++i)
		i->first*= R(-1);
      }
};

template<class R>
AnyType mM2L3 (Stack , const AnyType & pp)
{
    minusMat<R> mpp(to(GetAny<Matrice_Creuse<R> *>(pp)));
    return SetAny<minusMat<R> >(mpp);
}
/* template<class R>
AnyType mmM2L3 (Stack , const AnyType & pp)
{
    minusMat<R> &  p(GetAny<minusMat<R> >(pp));
    minusMat<R> mpp(p.l);
    delete  p.l;
    return mpp.l;
}

template<class R>
AnyType mmM2L3c (Stack , const AnyType & pp)
{
    list<triplet<R,MatriceCreuse<R> *,bool> >  *  p(GetAny<minusMat<R> >(pp))
    minusMat<R> mpp(p.l);
    delete  p.l;
    return mpp.l;
}
*/
template <class R>
void AddSparseMat()
{
 aType tkrp = atype<KN<R> *>(); 
 SetMatrix_Op<R>::btype = Dcl_Type<const  SetMatrix_Op<R> * >();
 Dcl_Type<TheDiagMat<R> >();
 Dcl_Type<TheCoefMat<R> >(); // Add FH oct 2005
 Dcl_Type< map< pair<int,int>, R> * >(); // Add FH mars 2005 
 Dcl_Type<  minusMat<R>  >(); // Add FJH mars 2007
 
 basicForEachType * t_MC=atype<  Matrice_Creuse<R>* >();
// basicForEachType * t_MCt=atype<  Matrice_Creuse_Transpose<R> >();
// basicForEachType * t_lM=atype< list<triplet<R,MatriceCreuse<R> *,bool> > * >();
// basicForEachType * t_nM=atype<  minusMat<R> >();
 
 basicForEachType * t_MM=atype<map< pair<int,int>, R> * >();
 
TheOperators->Add("*", 
        new OneBinaryOperator<Op2_mulvirtAv<typename VirtualMatrice<R>::plusAx,Matrice_Creuse<R>*,KN_<R> > >,
        new OneBinaryOperator<Op2_mulvirtAv<typename VirtualMatrice<R>::plusAtx,Matrice_Creuse_Transpose<R>,KN_<R> > >,
        new OneBinaryOperator<Op2_mulvirtAv<typename VirtualMatrice<R>::solveAxeqb,Matrice_Creuse_inv<R>,KN_<R> > >     
        );

TheOperators->Add("*", 
        new OneBinaryOperator<Op2_mulvirtAv<typename VirtualMatrice<R>::plusAx,Matrice_Creuse<R>*,KN_<R> > >( 0  ,tkrp),
        new OneBinaryOperator<Op2_mulvirtAv<typename VirtualMatrice<R>::plusAtx,Matrice_Creuse_Transpose<R>,KN_<R> > >( 0 ,tkrp),
        new OneBinaryOperator<Op2_mulvirtAv<typename VirtualMatrice<R>::solveAxeqb,Matrice_Creuse_inv<R>,KN_<R> > >( 0 ,tkrp)     
        );
        
TheOperators->Add("^", new OneBinaryOperatorA_inv<R>());
  
// matrix new code   FH (Houston 2004)        
 TheOperators->Add("=",
//       new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,const MatrixInterpolation::Op*,E_F_StackF0F0>(SetMatrixInterpolation),
       new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,const Matrix_Prod<R,R>,E_F_StackF0F0>(ProdMat<R,R,R>),
       new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,KN<R> *,E_F_StackF0F0>(DiagMat<R>),       
       new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,Matrice_Creuse_Transpose<R>,E_F_StackF0F0>(CopyTrans<R,R>), 
       new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,Matrice_Creuse<R>*,E_F_StackF0F0>(CopyMat<R,R>) ,
       new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,KNM<R>*,E_F_StackF0F0>(MatFull2Sparse<R>) ,
       new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,map< pair<int,int>, R> * ,E_F_StackF0F0>(MatMap2Sparse<R>) ,
       new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,list<triplet<R,MatriceCreuse<R> *,bool> > *,E_F_StackF0F0>(CombMat<R>) ,
       new OneOperatorCode<BlockMatrix<R> >()
       
       );
       
 TheOperators->Add("<-",
//       new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,const MatrixInterpolation::Op*,E_F_StackF0F0>(SetMatrixInterpolation),
       new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,const Matrix_Prod<R,R>,E_F_StackF0F0>(ProdMat<R,R,R>),
       new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,KN<R> *,E_F_StackF0F0>(DiagMat<R>)  ,
       new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,Matrice_Creuse_Transpose<R>,E_F_StackF0F0>(CopyTrans<R,R>), 
       new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,Matrice_Creuse<R>*,E_F_StackF0F0>(CopyMat<R,R>) ,
       new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,KNM<R>*,E_F_StackF0F0>(MatFull2Sparse<R>) ,
       new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,map< pair<int,int>, R> * ,E_F_StackF0F0>(MatMap2Sparse<R>) ,
       new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,list<triplet<R,MatriceCreuse<R> *,bool> > *,E_F_StackF0F0>(CombMat<R>), 
       new OneOperatorCode<BlockMatrix<R> >()
       
       );
TheOperators->Add("*", 
        new OneBinaryOperator<Op2_pair<Matrix_Prod<R,R>,Matrice_Creuse<R>*,Matrice_Creuse<R>*> >,
        new OneBinaryOperator<Op2_pair<Matrix_Prod<R,R>,Matrice_Creuse_Transpose<R>,Matrice_Creuse<R>* > >, 
        new OneBinaryOperator<Op2_pair<Matrix_Prod<R,R>,Matrice_Creuse_Transpose<R>,Matrice_Creuse_Transpose<R> > >,
        new OneBinaryOperator<Op2_pair<Matrix_Prod<R,R>,Matrice_Creuse<R>*,Matrice_Creuse_Transpose<R> > > ,
        new OneBinaryOperator<Op2_ListCM<R> >  , 
        new OneBinaryOperator<Op2_ListMC<R> >  ,
	new OneBinaryOperator<Op2_ListCMt<R> >  , 
        new OneBinaryOperator<Op2_ListMtC<R> >  
		  
        );
TheOperators->Add("+", 
        new OneBinaryOperator<Op2_ListCMCMadd<R> >,
        new OneBinaryOperator<Op2_ListCMMadd<R> >,
        new OneBinaryOperator<Op2_ListMCMadd<R> >,
//	new OneBinaryOperator<Op2_ListCMCMadd<R> >(t_MCt,t_lM),
//	new OneBinaryOperator<Op2_ListCMCMadd<R> >(t_MC,t_lM),
        new OneBinaryOperator<Op2_ListMMadd<R> >
       
       ); 
 TheOperators->Add("-",  
	 new OneUnaryOperator<Op1_LCMd<R> >
     );
 Add<Matrice_Creuse<R> *>("n",".",new OneOperator1<long,Matrice_Creuse<R> *>(get_mat_n<R>) );
 Add<Matrice_Creuse<R> *>("m",".",new OneOperator1<long,Matrice_Creuse<R> *>(get_mat_m<R>) );
 Add<Matrice_Creuse<R> *>("nbcoef",".",new OneOperator1<long,Matrice_Creuse<R> *>(get_mat_nbcoef<R>) );
 
 
 Add<Matrice_Creuse<R> *>("diag",".",new OneOperator1<TheDiagMat<R> ,Matrice_Creuse<R> *>(thediag<R>) );
 Add<Matrice_Creuse<R> *>("coef",".",new OneOperator1<TheCoefMat<R> ,Matrice_Creuse<R> *>(thecoef<R>) );

// Add<Matrice_Creuse<R> *>("setdiag",".",new OneOperator2<long,Matrice_Creuse<R> *,KN<R> *>(set_diag<R>) );
 TheOperators->Add("=", new OneOperator2<KN<R>*,KN<R>*,TheDiagMat<R> >(get_mat_daig<R>) );
 TheOperators->Add("=", new OneOperator2<TheDiagMat<R>,TheDiagMat<R>,KN<R>*>(set_mat_daig<R>) );
 
// TheOperators->Add("=", new OneOperator2<KN<R>*,KN<R>*,TheDiagMat<R> >(get_mat_daig<R>) );
// TheOperators->Add("=", new OneOperator2<TheDiagMat<R>,TheDiagMat<R>,KN<R>*>(set_mat_daig<R>) );
// ADD oct 2005
 TheOperators->Add("=", new OneOperator2<KN<R>*,KN<R>*,TheCoefMat<R> >(get_mat_coef<R>) );
 TheOperators->Add("=", new OneOperator2<TheCoefMat<R>,TheCoefMat<R>,KN<R>*>(set_mat_coef<R>) );
 
// TheOperators->Add("=", new OneOperator2<KN<R>*,KN<R>*,TheCoefMat<R> >(get_mat_coef<R>) );
// TheOperators->Add("=", new OneOperator2<TheCoefMat<R>,TheCoefMat<R>,KN<R>*>(set_mat_coef<R>) );
 
 Global.Add("set","(",new SetMatrix<R>);
 //Global.Add("psor","(",new  OneOperatorCode<Psor<R> > );
 
 atype<Matrice_Creuse<R> * >()->Add("(","",new OneOperator3_<R*,Matrice_Creuse<R> *,long,long >(get_elementp2mc<R>));
 
 atype<KNM<R>*>()->Add("(","",new OneOperator3_<map< pair<int,int>, R> *,KNM<R>*,Inv_KN_long,Inv_KN_long >(Matrixfull2mapIJ_inv<R>));
 atype<KNM<R>*>()->Add("(","",new OneOperator3_<map< pair<int,int>, R> *,KNM<R>*,KN_<long>,KN_<long> >(Matrixfull2mapIJ<R>));
 
 atype<outProduct_KN_<R>*>()->Add("(","",new OneOperator3_<map< pair<int,int>, R> *,outProduct_KN_<R>*,Inv_KN_long,Inv_KN_long >(Matrixoutp2mapIJ_inv<R>));
 atype<outProduct_KN_<R>*>()->Add("(",