lgmat.cpp

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

   typedef Matrice_Creuse<R> * Result;
   int N,M; 
   Expression emat; 
   Expression ** e_Mij;
   int ** t_Mij;
   BlockMatrix(const basicAC_F0 & args) ;
   ~BlockMatrix() ;
      
    static ArrayOfaType  typeargs() { return  ArrayOfaType(atype<Matrice_Creuse<R>*>(),atype<E_Array>());}
    static  E_F0 * f(const basicAC_F0 & args){
	if(IsRawMat<R>(args)) return new RawMatrix<R>(args);
	else return new BlockMatrix(args);
    }  
    AnyType operator()(Stack s) const ;
    
};

template<typename R>  
map< pair<int,int>, R> *Matrixfull2mapIJ_inv (KNM<R>   * const & pa,const Inv_KN_long & iii,const Inv_KN_long & jjj)
{
   const  KN_<long> &ii(iii), &jj(jjj);
   const KNM<R> & a(*pa);
   int N=a.N(),M=a.M();
   long n = ii(SubArray(N)).max()+1;
   long m= jj(SubArray(M)).max()+1;
   
/*
  long minn = ii(SubArray(N)).min()+1;
  long minm= jj(SubArray(M)).min()+1;
  if ( !(0 <= minn && 0 <=  minm) ) 
  {
  cerr << " Out of Bound  in A(I^-1,J^1) :  "<< minn << " " << minm <<" =>  negative value!! " << endl;
  ExecError("Out of Bound Error");
  }
*/
   
  // cout << "  ### n m " << n << " " << m << endl; 
   map< pair<int,int>, R> *pA= new map< pair<int,int>, R>;
   map< pair<int,int>, R> & A(*pA);
   A[make_pair(n-1,m-1)] = R(); // Hack to be sure that the last term existe 
  
   for (long i=0;i<N;++i)
    for (long j=0;j<M;++j)
     { R aij=a(i,j);
       //cout << i << " " << j << " :: " << ii[i] << " " << jj[j] << " = " << aij << endl;
       if(ii[i]>=0 && jj[j]>=0 && norm(aij)>1e-40) 
         A[make_pair(ii[i],jj[j])] += aij;
     }
      
  return pA;
}

template<typename R>  
map< pair<int,int>, R> *Matrixfull2mapIJ (KNM<R>   * const & pa,const KN_<long> & ii,const  KN_<long> & jj)
{
   const KNM<R> & a(*pa);
   int N=a.N(),M=a.M();
   long n = ii.N();
   long m= jj.N();
  // cout << "  ### n m " << n << " " << m << endl; 
   map< pair<int,int>, R> *pA= new map< pair<int,int>, R>;
   map< pair<int,int>, R> & A(*pA);
   A[make_pair(n-1,m-1)] += R(); // Hack to be sure that the last term existe 
  
   for (long il=0;il<N;++il)
    for (long jl=0;jl<M;++jl)
     { 
       long i = ii[il];
       long j = jj[jl];
       if( i>=0 && j >=0) {
          if ( !(0 <= i && i < N && 0 <= j && j < M ) )
            {
              cerr << " Out of Bound  in A(I,J) : " << i << " " << j << " not in " << "[0,"<<N<<"[x[0," << M << "[ \n";
              ExecError("Out of Bound Error");
             }
       
          R aij=a(i,j);
       //cout << i << " " << j << " :: " << ii[i] << " " << jj[j] << " = " << aij << endl;
         if (norm(aij)>1e-40) 
           A[make_pair(il,jl)] += aij;
       }
     }
      
  return pA;
}

template<class R>
AnyType Matrixfull2map (Stack , const AnyType & pp)
{
   const KNM<R> & a(*GetAny<KNM<R> *>(pp));
   int N=a.N(),M=a.M();
   int n = N;
   int m= M;
   map< pair<int,int>, R> *pA= new map< pair<int,int>, R>;
   map< pair<int,int>, R> & A(*pA);
   A[make_pair(n-1,m-1)] = R(); // Hack to be sure that the last term existe 
  
   for (int i=0;i<N;++i)
    for (int j=0;j<M;++j)
     { R aij=a(i,j);
     if (norm(aij)>1e-40) 
      A[make_pair(i,j)] += aij;
     }
      
  return pA;
}


template<class R>
map< pair<int,int>, R> *Matrixoutp2mapIJ_inv (outProduct_KN_<R>   * const & pop,const Inv_KN_long & iii,const Inv_KN_long & jjj)
{
   const KN_<long> &ii(iii), &jj(jjj);
   const outProduct_KN_<R> & op(*pop);
   long  N=op.a.N(),M=op.b.N();
   long  n = ii(SubArray(N)).max()+1;
   long m= jj(SubArray(M)).max()+1;
/*
   long minn = ii(SubArray(N)).min()+1;
   long minm= jj(SubArray(M)).min()+1;
     if ( !(0 <= minn && 0 <=  minm) ) 
        {
            cerr << " Out of Bound  in A(I^-1,J^1) :  "<< minn << " " << minm <<" =>  negative value!! " << endl;
            ExecError("Out of Bound Error");
        }
 */
   map< pair<int,int>, R> *pA= new map< pair<int,int>, R>;
   map< pair<int,int>, R> & A(*pA);
   A[make_pair(n-1,m-1)] = R(); // Hack to be sure that the last term existe 
  
   for (int i=0;i<N;++i)
    for (int j=0;j<M;++j)
     { 
       R aij=op.a[i]*conj(op.b[j]);
       if(ii[i]>=0 && jj[j]>=0 && norm(aij)>1e-40) 
//       if (norm(aij)>1e-40 &) 
          A[make_pair(ii[i],jj[j])] += aij;
     }   
  delete pop;
    
  return pA;
}


template<class R>
map< pair<int,int>, R> *Matrixmapp2mapIJ1 (map< pair<int,int>, R> *const &  B,const Inv_KN_long & iii,const Inv_KN_long  & jjj)
{
    const KN_<long> &ii(iii), &jj(jjj);  
    typedef typename map< pair<int,int>, R>::const_iterator It;
    
    map< pair<int,int>, R> *pA= new map< pair<int,int>, R>;
    map< pair<int,int>, R> & A(*pA);
    int n=0,m=0;
    // hack:  the last element must exist in the map  to set matrix size 
/*    
    It lastb= --B->end(); // le last element
    
    if( lastb != B->end() )	
    { 
	nb = last->first.first+1; 
	mb=last->first.second+1;
    } */
    int N=ii.N(),M=jj.N();
   // A[make_pair(n-1,m-1)] = R(); // Hack to be sure that the last term existe 
    
    for (It k=B->begin();k!=B->end();++k)
    {
	int il =  k->first.first;
	int jl =  k->first.second;
	if ( !( 0 <= il && il < N && 0 <= jl && jl < M )  )
	{
	    cerr << " Out of Bound  in (Map)(I,J) : " << il << " " << jl << " not in " << "[0,"<<N<<"[x[0," << M << "[ \n";
	    ExecError("Out of Bound Error");
	}
	int i=ii(il);
	int j=jj(jl);
	n=max(i,n);
	m=max(j,m);
	R aij =k->second;
	if(i >=0 && j>=0) 
	  A[make_pair(i,j)] += aij;
    }
     A[make_pair(n,m)] += R(); // Hack to be sure that the last term existe 
     delete B;
    
    return pA;
}

template<class R>
map< pair<int,int>, R> *Matrixmapp2mapIJ (map< pair<int,int>, R> *const &  B,const KN_<long> & ii,const KN_<long>  & jj)
{
    
    typedef typename map< pair<int,int>, R>::const_iterator It;
    typedef typename multimap< int,int>::iterator  MI;
    
    map< pair<int,int>, R> *pA= new map< pair<int,int>, R>;
    map< pair<int,int>, R> & A(*pA);
    multimap< int,int > I,J;
    int N=ii.N(),M=jj.N();
    for (int i=0;i<N;++i)
	if(ii[i]>=0)
	  I.insert(make_pair(ii[i],i));
    for (int j=0;j<M;++j)
	if(jj[j]>=0)
	    J.insert(make_pair(jj[j],j));
    int n=0,m=0;

    for (It k=B->begin();k!=B->end();++k)
    {
	int il =  k->first.first;
	int jl =  k->first.second;
	R aij =k->second;
	pair<MI,MI> PPi=I.equal_range(il);
	pair<MI,MI> PPj=J.equal_range(jl);
	for(MI pi=PPi.first ; pi !=  PPi.second; ++pi)
	{
	    int i=pi->second;
	    for(MI pj=PPj.first ; pj !=  PPj.second; ++pj)
	    { 
		int j=pj->second;
		n=max(i,n);
	        m=max(j,m);	
	       if(i >=0 && j>=0) 
	         A[make_pair(i,j)] += aij;
	    }
	}   
    }
    A[make_pair(n,m)] =+ R(); // Hack to be sure that the last term existe 
    delete B;
    
    return pA;
}

template<class R>
map< pair<int,int>, R> *Matrixoutp2mapIJ (outProduct_KN_<R>   * const & pop,const KN_<long> & ii,const KN_<long>  & jj)
{
   const outProduct_KN_<R> & op(*pop);
   long N=op.a.N(),M=op.b.N();
   long n=ii.N(),m=jj.N();
   
   map< pair<int,int>, R> *pA= new map< pair<int,int>, R>;
   map< pair<int,int>, R> & A(*pA);
   A[make_pair(n-1,m-1)] = R(); // Hack to be sure that the last term existe 
   
   for (long il=0;il<n;++il)
    for (long jl=0;jl<m;++jl)
     { 
       long i = ii[il];
       long j = jj[jl];
       if(i>=0 && j >=0)
        {
               if ( !( 0 <= i && i < N && 0 <= j && j < M )  )
                {
                    cerr << " Out of Bound  in (a*b')(I,J) : " << i << " " << j << " not in " << "[0,"<<N<<"[x[0," << M << "[ \n";
                    ExecError("Out of Bound Error");
                }
               R aij=op.a[i]*conj(op.b[j]);
               if (norm(aij)>1e-40) 
                  A[make_pair(il,jl)] += aij;
               }
     }   
  delete pop;
    
  return pA;
}


template<class R>
AnyType Matrixoutp2map (Stack , const AnyType & pp)
{
   const outProduct_KN_<R> & op(*GetAny<outProduct_KN_<R> *>(pp));
   long N=op.a.N(),M=op.b.N();
   long n = N;
   long m= M;
   map< pair<int,int>, R> *pA= new map< pair<int,int>, R>;
   map< pair<int,int>, R> & A(*pA);
   A[make_pair(n-1,m-1)] = R(); // Hack to be sure that the last term existe 
  
   for (long i=0;i<N;++i)
    for (long j=0;j<M;++j)
     { 
      R aij=op.a[i]*conj(op.b[j]);
      if (norm(aij)>1e-40) 
        A[make_pair(i,j)] += aij;
     } 
  delete &op;        
  return pA;
}


template<typename R>  BlockMatrix<R>::~BlockMatrix() 
{  
    if (e_Mij)
    {  cout << " del Block matrix "<< this << " " << e_Mij <<" N = " << N << " M = " << M << endl;
	for (int i=0;i<N;i++)
	{ delete [] e_Mij[i];
	    delete [] t_Mij[i];
	}
	delete [] e_Mij;
	delete [] t_Mij;
	N=0;
	M=0;
	e_Mij=0;
	t_Mij=0; }
}   

template<typename R>  RawMatrix<R>::RawMatrix(const basicAC_F0 & args) 
{
    args.SetNameParam();
    emat = args[0];
    
    const E_Array & ee= *dynamic_cast<const E_Array*>((Expression) args[1]);
    
    int N=ee.size();
    if (N==1)
    {
	C_F0 c0(ee[0]);
	coef=to<KN_<R> >(ee[0]);
	lig=0;
	col=0;
    }
    
    else if (N==3)
    {
	C_F0 c0(ee[0]),c1(ee[1]),c2(ee[2]);
	coef=to<KN_<R> >(ee[2]);
	lig=to<KN_<long> >(ee[0]);
	col=to<KN_<long> >(ee[1]);
	
    }
    
    
}
template<typename R>  BlockMatrix<R>::BlockMatrix(const basicAC_F0 & args) 
{   
    N=0;
    M=0;
    args.SetNameParam();
    emat = args[0];
    const E_Array & eMij= *dynamic_cast<const E_Array*>((Expression) args[1]);
    N=eMij.size();
    int err =0;
    for (int i=0;i<N;i++)
    {
        const E_Array* emi= dynamic_cast<const E_Array*>((Expression)  eMij[i]);
        if (!emi) err++;
        else
        { 
	    if ( i==0) 
		M = emi->size();
	    else
		if(M != emi->size()) err++;
        }
    }
    if (err) {
	CompileError(" Block matrix : [[ a, b, c], [ a,b,c ]] or Raw Matrix [a] or [ l, c, a ] ");
    }
    assert(N && M);
    e_Mij = new Expression * [N];
    t_Mij = new int * [N];
    for (int i=0;i<N;i++)
    {
	const E_Array li= *dynamic_cast<const E_Array*>((Expression)  eMij[i]);
	
	e_Mij[i] =  new Expression [M];
	t_Mij[i] = new int [M]; 
	for (int j=0; j<M;j++)  
	{
	    C_F0 c_Mij(li[j]);
	    Expression eij=c_Mij.LeftValue();
	    aType rij = c_Mij.left();
	    if ( rij == atype<long>() &&  eij->EvaluableWithOutStack() )
	    {
		long contm = GetAny<long>((*eij)(0));
		/*  prev  version
		if(contm !=0) 
		CompileError(" Block matrix , Just 0 matrix");
		e_Mij[i][j]=0;
		t_Mij[i][j]=0;*/
		if(contm==0)
		{
		    e_Mij[i][j]=0;
		    t_Mij[i][j]=0;
		}
		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 , Just 0 matrix");
	    }
	    else if ( rij ==  atype<Matrice_Creuse<R> *>()) 
	    {
		e_Mij[i][j]=eij;
		t_Mij[i][j]=1;
	    } 
	    else if ( rij ==  atype<Matrice_Creuse_Transpose<R> >()) 
	    {
		e_Mij[i][j]=eij;
		t_Mij[i][j]=2;
	    } 
	    else if ( atype<KN_<R> >()->CastingFrom(rij) )
	    {  
		e_Mij[i][j]=to<KN_<R> >(c_Mij);
		t_Mij[i][j]=3;
		
	    } 
	    else if ( atype<Transpose<KN_<R> > >()->CastingFrom(rij) )
	    {  
		
		e_Mij[i][j]=to<Transpose<KN_<R> > >(c_Mij);
		t_Mij[i][j]=4;
	    } 
	    else if ( atype<KNM<R> *  >()->CastingFrom(rij) )
	    {  
		
		e_Mij[i][j]=to<KNM<R> * >(c_Mij);
		t_Mij[i][j]=5;
	    }
	    else if ( atype<Transpose< KNM<R> * > >()->CastingFrom(rij) )
	    {  
		
		e_Mij[i][j]=to<Transpose<KNM<R> *> >(c_Mij);