lgmat.cpp

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

  TypeSolveMat    *typemat=&tmat;
  bool initmat=true;
  int umfpackstrategy=0; 
  if (nargs[0]) initmat= ! GetAny<bool>((*nargs[0])(stack));
  if (nargs[1]) typemat= GetAny<TypeSolveMat *>((*nargs[1])(stack));
  if (nargs[2]) eps= GetAny<double>((*nargs[2])(stack));
  // 3 precon 
  if (nargs[4]) NbSpace= GetAny<long>((*nargs[4])(stack));
  if (nargs[6]) tgv= GetAny<double>((*nargs[6])(stack));
  if (nargs[7]) factorize= GetAny<bool>((*nargs[7])(stack));
  
  if (nargs[8]) umfpackstrategy = GetAny<long>((*nargs[8])(stack)); 
  if (nargs[9]) tol_pivot = GetAny<double>((*nargs[9])(stack)); 
  if (nargs[10]) tol_pivot_sym = GetAny<double>((*nargs[10])(stack)); 
  if (nargs[11]) itmax = GetAny<long>((*nargs[11])(stack)); //  frev 2007 OK
   
   if(A->typemat.profile != typemat->profile) 
   {
     cerr << " type of matrix " << A->typemat<<endl;
     cerr << " type of matrix for solver " <<*typemat<<endl;
     
     ExecError(" Set incompatibility between solver and type of matrix");
   }
  if( factorize ) {
    MatriceProfile<R> * pf = dynamic_cast<MatriceProfile<R> *>((MatriceCreuse<R> *) A->A);
    assert(pf);
    switch (typemat->t) {
    case TypeSolveMat::LU: pf->LU(Abs(eps));break;
    case TypeSolveMat::CROUT: pf->crout(Abs(eps));break;
    case TypeSolveMat::CHOLESKY: pf->cholesky(Abs(eps));break;
    default: ExecError("Sorry no factorization for this type for matrix"); 
    }
    
  }    
  SetSolver<R>(stack,*A->A,typemat,VF,eps,NbSpace,itmax,precon,umfpackstrategy,tgv,tol_pivot,tol_pivot_sym);

  return Nothing; 
}



bool SetDefaultSolver()
{

#ifdef HAVE_LIBUMFPACKXXXXXX
    if(verbosity>1)
	cout << " SetDefault sparse solver to UMFPACK" << endl;
    DefSparseSolver<double>::solver  =BuildSolverUMFPack;
    DefSparseSolver<Complex>::solver =BuildSolverUMFPack;    
#else
    if(verbosity>1)
	cout << " SetDefault sparse solver to GMRES (no UMFPACK)" << endl;
    DefSparseSolver<double>::solver  =BuildSolverGMRES<double>;
    DefSparseSolver<Complex>::solver =BuildSolverGMRES<Complex>;    
#endif
    return true;

}


AnyType SetMatrixInterpolation(Stack,Expression ,Expression);

//------

template<class R>
void BuildCombMat(map< pair<int,int>, R> & mij,const KNM_<R> & A, int ii00=0,int jj00=0,R coef=R(1.),bool cnj=false)
{
  double eps0=numeric_limits<double>::min();
  int i,j;
  int n = A.N(),m=A.M();
  for ( i=0;i<n;i++)
   for ( j=0;j<m;j++)
          {
           R cij=coef*A(i,j);
           if (cnj)  cij = conj(cij); 
           if(norm(cij) >eps0)
             mij[ij_mat(false,ii00,jj00,i,j)] += cij;
         
   }

}

void buildInterpolationMatrix(MatriceMorse<R> * m,const FESpace & Uh,const FESpace & Vh,void *data)
{  //  Uh = Vh 

  int op=op_id; //  value of the function
  bool transpose=false;
  bool inside=false;
   int * iU2V=0;  
  if (data)
   {
     int * idata=static_cast<int*>(data);
     transpose=idata[0];
     op=idata[1];
     inside=idata[2];
       iU2V= idata + 5;
     ffassert(op>=0 && op < 4);
   }
  if(verbosity>2) 
    {
      cout << " -- buildInterpolationMatrix   transpose =" << transpose << endl
           << "              value, dx , dy          op = " << op << endl
           << "                            just  inside = " << inside << endl;
    }
  using namespace Fem2D;
  int n=Uh.NbOfDF;
  int mm=Vh.NbOfDF;
  if(transpose) Exchange(n,mm);
  m->symetrique = false;
  m->dummy=false;
  m->a=0;
  m->lg=0;
  m->cl=0;
  m->nbcoef=0;
  m->n=n;
  m->m=mm;
  int n1=n+1;
  const  Mesh & ThU =Uh.Th; // line 
  const  Mesh & ThV =Vh.Th; // colunm
  bool samemesh =  &Uh.Th == &Vh.Th;  // same Mesh
  int thecolor =0;
  //  int nbn_u = Uh.NbOfNodes;
  //int nbn_v = Vh.NbOfNodes;
  
  int nbcoef =0;
  
  KN<int> color(ThV.nt);
  KN<int> mark(n);
  mark=0;
  
  int *lg = new int [n1];
  int * cl = 0;
  double *a=0;
  
  color=thecolor++;
  FElement Uh0 = Uh[0];
  FElement Vh0 = Vh[0];
  
  FElement::aIPJ ipjU(Uh0.Pi_h_ipj()); 
  FElement::aR2  PtHatU(Uh0.Pi_h_R2()); 
  
  //  FElement::aIPJ ipjV(Vh0.Pi_h_ipj()); 
  // FElement::aR2  PtHatV(Vh0.Pi_h_R2()); 
  
  int nbdfVK= Vh0.NbDoF();
  //  int nbdfUK= Uh0.NbDoF();
  int NVh= Vh0.N;
  int NUh= Uh0.N;
  
  //ffassert(NVh==NUh); 
  
  
  int nbp= PtHatU.N(); // 
  //  int nbc= ipjU.N(); // 
  KN<R2> PV(nbp);   //  the PtHat in ThV mesh
  KN<int> itV(nbp); // the Triangle number
  KN<bool> intV(nbp); // ouside or not 
  KN<R>   AipjU(ipjU.N());
  
  KNM<R> aaa(nbp,nbdfVK); 


   const R eps = 1.0e-10;
   const int sfb1=Vh0.N*last_operatortype*Vh0.NbDoF();
   KN<R> kv(sfb1*nbp);
   R * v = kv;
    KN<int> ik(nbp); // the Triangle number
//   KNMK_<R> fb(v+ik[i],VhV0.NbDoF,VhV0.N,1); //  the value for basic fonction

   bool whatd[last_operatortype];
   for (int i=0;i<last_operatortype;i++) 
     whatd[i]=false;
   whatd[op]=true; // the value of function
   KN<bool> fait(Uh.NbOfDF);
   fait=false;
  map< pair<int,int> , double > sij; 
  
  for (int step=0;step<2;step++) 
   {
      
	  for (int it=0;it<ThU.nt;it++)
	    {
	      thecolor++; //  change the current color
	      const Triangle & TU(ThU[it]);
	      FElement KU(Uh[it]);
	      KU.Pi_h(AipjU);
	      int nbkU = 0;
	      if (samemesh)
	        {
	          nbkU = 1; 
	          PV = PtHatU;
	          itV = it;
	        }
	      else 
	       {  
	          const Triangle *ts=0;
	          for (int i=0;i<nbp;i++)
	            {
	              bool outside;
	              ts=ThV.Find(TU(PtHatU[i]),PV[i],outside,ts); 
	              itV[i]= ThV(ts);
	              intV[i]=outside && inside; //  ouside and inside flag 
	            }
	       }
	       
	      
	      for (int p=0;p<nbp;p++)
	         { 

	              KNMK_<R> fb(v+p*sfb1,nbdfVK,NVh,last_operatortype); // valeur de fonction de base de Vh 
	              // ou:   fb(idf,j,0) valeur de la j composante de la fonction idf 
	              Vh0.tfe->FB(whatd,ThV,ThV[itV[p]],PV[p],fb);  
	          }

	      for (int i=0;i<ipjU.N();i++) 
	          { // pour tous le terme 
	           const FElement::IPJ &ipj_i(ipjU[i]);
	          // assert(ipj_i.j==0); // car  Vh.N=0
	           int dfu = KU(ipj_i.i); // le numero de df global 
	           if(fait[dfu]) continue;
	           int jU = ipj_i.j; // la composante dans U
	           int p=ipj_i.p;  //  le points
	           if (intV[p]) continue; //  ouside and inside flag => next 
	           R aipj = AipjU[i];
	           FElement KV(Vh[itV[p]]);
		      int jV=jU;
		      if(iU2V) jV=iU2V[jU];
		    
		    if(jV>=0 && jV<NVh)
			{
			    KNMK_<R> fb(v+p*sfb1,nbdfVK,NVh,last_operatortype); 
			    KN_<R> fbj(fb('.',jV,op)); 
			    
			    for (int idfv=0;idfv<nbdfVK;idfv++) 
				if (Abs(fbj[idfv])>eps) 
				  {
				      int dfv=KV(idfv);
				      int ii=dfu, jj=dfv;
				      if(transpose) Exchange(ii,jj);
				      // le term dfu,dfv existe dans la matrice
				      R c= fbj[idfv]*aipj;
				      //  cout << " Mat inter " << i << " , "<< j << " = " << c << " " <<step << " " << it << " " <<  endl; 
				      if(Abs(c)>eps)
					  //
					  sij[make_pair(ii,jj)] = c;
				      /*   
				       if(step==0)
				       sij.insert(make_pair(i,j));
				       else	                  
				       (*m)(i,j)=c;
				       */
				  }
			}
	                      
	          }
	          
	      for (int df=0;df<KU.NbDoF();df++) 
	          {  
	           int dfu = KU(df); // le numero de df global 
	           fait[dfu]=true;
	           }
	       
	       
	    }
	    if (step==0)
	     {
	         nbcoef = sij.size();
	         cl = new int[nbcoef];
	         a = new double[nbcoef];
	         int k=0;
	         for(int i=0;i<n1;i++)
	           lg[i]=0;
	          
	         for (map<pair<int,int>, double >::iterator kk=sij.begin();kk!=sij.end();++kk)
	          { 
	            int i= kk->first.first;
	            int j= kk->first.second;
	           // cout << " Mat inter " << i << " , "<< j  << endl;
	            cl[k]=j;
	            a[k]= kk->second;	            
	            lg[i+1]=++k;
	          }
	         assert(k==nbcoef);
	         //  on bouche les ligne vide   lg[i]=0; 
	         //  lg est un tableau croissant  =>
	         for(int i=1;i<n1;i++)
	              lg[i]=max(lg[i-1],lg[i]) ;
	         m->lg=lg;
             m->cl=cl;
             m->a=a;
             m->nbcoef=nbcoef;
             fait=false;
	     }
    }
    sij.clear();
  //assert(0); // a faire to do
}

MatriceMorse<R> *  buildInterpolationMatrix1(const FESpace & Uh,const KN_<double> & xx,const KN_<double> & yy ,int *data)
{  //  Uh = Vh 

  int op=op_id; //  value of the function
  int icomp=0;
  bool transpose=false;
  bool inside=false;
  if (data)
   {
     transpose=data[0];
     op=data[1];
     inside=data[2];
     icomp = data[3];
     ffassert(op>=0 && op < 4);
     
   }
  if(verbosity>2) 
    {
      cout << " -- buildInterpolationMatrix   transpose =" << transpose << endl
           << "              value, dx , dy          op = " << op << endl
           << "              composante                 = " << icomp << endl
           << "                            just  inside = " << inside << endl;
    }
  using namespace Fem2D;
  int n=Uh.NbOfDF;
  int mm=xx.N();
  int nbxx= mm;
  if(transpose) Exchange(n,mm);
  const  Mesh & ThU =Uh.Th; // line 
  
   
  FElement Uh0 = Uh[0];
  
   
   
 int nbdfUK= Uh0.NbDoF();
 int NUh= Uh0.N;
 
 ffassert(icomp < NUh && icomp >=0); 
    

   const int sfb1=Uh0.N*last_operatortype*Uh0.NbDoF();
   KN<R> kv(sfb1);
   R * v = kv;
   const R eps = 1.0e-10;

   bool whatd[last_operatortype];
   for (int i=0;i<last_operatortype;i++) 
     whatd[i]=false;
   whatd[op]=true; // the value of function
   KN<bool> fait(Uh.NbOfDF);
   fait=false;
   map< pair<int,int> , double > sij; 
   R2 Phat;
   bool outside;
   
   for(int ii=0;ii<nbxx;ii++)
   {
     const Triangle *ts=ThU.Find(R2(xx[ii],yy[ii]),Phat,outside); 
     if(outside && !inside) continue;
     int it = ThU(ts); 
     FElement KU(Uh[it]);
     KNMK_<R> fb(v,nbdfUK,NUh,last_operatortype);
     Uh0.tfe->FB(whatd,ThU,ThU[it],Phat,fb);  
     KN_<R> Fwi(fb('.',icomp,op));
     for (int idfu=0;idfu<nbdfUK;idfu++) 
       {
        int  j = ii;
        int  i = KU(idfu);
        if(transpose) Exchange(i,j);
        R c = Fwi(idfu);
	if(Abs(c)>eps)
	  sij[make_pair(i,j)] = c;
      }
      }
  
   MatriceMorse<R> * m = new MatriceMorse<R>(n,mm,sij,false);
    sij.clear();
   return m;
  //assert(0); // a faire to do
}


AnyType SetMatrixInterpolation(Stack stack,Expression emat,Expression einter)
{
  using namespace Fem2D;
  
  Matrice_Creuse<R> * sparse_mat =GetAny<Matrice_Creuse<R>* >((*emat)(stack));
  const MatrixInterpolation::Op * mi(dynamic_cast<const MatrixInterpolation::Op *>(einter));
  ffassert(einter);
  int data[ MatrixInterpolation::Op::n_name_param+100];
  data[0]=mi->arg(0,stack,false); // transpose not
  data[1]=mi->arg(1,stack,(long) op_id); ; // get just value
  data[2]=mi->arg(2,stack,false); ; // get just value
  data[3]=mi->arg(3,stack,0L); ; // get just value
  KN<long> U2Vc;
  U2Vc= mi->arg(4,stack,U2Vc); ;
  if( mi->c==0)
  { // old cas 
  pfes * pUh = GetAny< pfes * >((* mi->a)(stack));
  pfes * pVh = GetAny<  pfes * >((* mi->b)(stack));
  FESpace * Uh = **pUh;
  FESpace * Vh = **pVh;
  int NVh =Vh->N;
  int NUh =Uh->N;
  ffassert(NUh< 100-MatrixInterpolation::Op::n_name_param);

      for(int i=0;i<NUh;++i)    
        data[5+i]=i;// 
      for(int i=0;i<min(NUh,(int) U2Vc.size());++i)    
	  data[5+i]= U2Vc[i];//
  if(verbosity>3)
	for(int i=0;i<NUh;++i)
	  {
	    cout << "The Uh componante " << i << " -> " << data[5+i] << "  Componante of Vh  " <<endl;
	  }
	  for(int i=0;i<NUh;++i)    
	if(data[5+i]>=NVh)
	  {
	      cout << "The Uh componante " << i << " -> " << data[5+i] << " >= " << NVh << " number of Vh Componante " <<endl;
	      ExecError("Interpolation incompability beetween componante ");
	  }
      
  ffassert(Vh);