lgfem.cpp

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

             //  assert(ii>=nbdfv);
               ndfv[j]=nbdfv ;
               j=ii;
               }             
             while (j!=nbdfv);
             if (verbosity > 100) 
             cout << "    ndf: " <<  j << " " <<  ii  << " <- " << nbdfv << " " <<  kkk <<  endl;
              nbdfv++;
           }
       return nbdfv;
}

bool  InCircularList(const int *p,int i,int k)
//  find k in circular list:  i , p[i], p[p[i]], ... 
{ 
    int j=i,l=0;    
    do {
	if (j==k) return true;
	ffassert(l++<10);
        j=p[j];
    } while (j!=i);    
    return false;
}

bool BuildPeriodic( 
  int nbcperiodic,
  Expression *periodic,
  const Mesh &Th,Stack stack,
  int & nbdfv, KN<int> & ndfv,int & nbdfe, KN<int> & ndfe) { 
  
/*
  build numbering of vertex form 0 to nbdfv-1
  and build numbering  of  edge form 0 to nbdfe-1
  we removing common vextex or common edge    
  --  we suppose one df by vertex 
      nbdfv number of df on vertex 
      ndfv[i]  given the numero of the df of the vertex 
  -- we suppose 1 df
*/  
   typedef Smallvect<int,2> int2;    
    if (nbcperiodic ) {
      
   //    KN<int> ndfv(Th.nv);
    //   KN<int> ndfe(Th.neb);
       ffassert(ndfv.N()==Th.nv);
       ffassert(ndfe.N()==Th.neb);
        
       MeshPoint *mp=MeshPointStack(stack),smp=*mp;   
       int n= nbcperiodic;
       if (verbosity >2)
         cout << " Nb of pair of periodic conditions: = " << n <<  endl;
       int * link1=0;
       int * link2=0;
       KN<int*> plk1(n),plk2(n);
       KN<int> nlk1(n),nlk2(n);
       KN<int> lab1(n),lab2(n);
#ifndef  HUGE_VAL      
       const double infty= numeric_limits<double>::infinity();
#else
       const double infty= HUGE_VAL;
#endif       
       int nblink1, nblink2;
       int *plink1 , *plink2;
        for (int step=0;step<2;step++)
         {
           nblink1=0,     nblink2=0;
           plink1=link1,  plink2=link2;
           for (int ip=0, k=0;ip<n;ip++,k+=4)
            {
               int label1=GetAny<long>((*periodic[k+0])(stack));
               int label2=GetAny<long>((*periodic[k+2])(stack));
               lab1[ip]=label1;
               lab2[ip]=label2;
               
               int l1=nblink1;
               int l2=nblink2;
               plk1[ip]= plink1;
               plk2[ip]= plink2;
               
               for (int ke=0;ke<Th.neb;ke++)
                {
                 if (Th.bedges[ke].lab==label1)
                  {
                    if (plink1) *plink1++=ke;
                    nblink1++;
                  }
                 else if (Th.bedges[ke].lab==label2)
                  {
                    if (plink2) *plink2++=ke;
                    nblink2++;
                   }
                 }
               nlk1[ip]= nblink1-l1;
               nlk2[ip]= nblink2-l2;              
            }
            if(step) break; // no reallocl 
            if (verbosity >3)
            cout << "  Periodic = " << nblink1 << " " << nblink2 << " step=" << step << endl;
            link1 = new int[nblink1];
            link2 = new int[nblink2];
            if(nblink1 != nblink2)
            {
             ExecError("Periodic:  the both number of edges is not the same ");
            }
         }
        if ( nblink1 >0) 
        {
        for (int ip=0, k=0;ip<n;ip++,k+=4)
          {
            map<int,int2> m;
            const int kk1=1,kk2=3;
            int label1=lab1[ip],label2=lab2[ip];
            int n1=nlk1[ip],n2=nlk2[ip];
            int *pke1=plk1[ip], *pke2=plk2[ip];
            double xmn=infty,xmx=-infty,hmn=infty;
            if (verbosity >1)
            cout << "  --Update: periodic  couple label1= " << label1 << ", n edges= " << n1 << "; "
                                          << ", label2= " << label2<<  ", n edges= " << n2 <<endl; 
            if (n1 != n2) ExecError("periodic BC:  the number of edges is not the same");
            for (int i1=0;i1<n1;i1++)
             {
              const BoundaryEdge & e =Th.bedges[pke1[i1]];
              if (e.lab==label1) 
                {           
                 mp->set(e[0].x,e[0].y);
                 double x0=GetAny<double>((*periodic[k+kk1])(stack));
                 mp->set(e[1].x,e[1].y);
                 double x1=GetAny<double>((*periodic[k+kk1])(stack));
		 if(verbosity>5)
			cout << "lab1:  e[" << pke1[i1] << "]  v0:   " <<  e[0].x << " " << e[0].y << "  s = " << x0
			<< "\t v1 " <<  e[0].x << " " << e[0].y << "  s = " << x1 << endl;
		    xmn=Min(x1,x0,xmn);
                 xmx=Max(x1,x0,xmx); 
                 hmn=Min(hmn,Abs(x1-x0));
               }                                
             }
            ffassert(hmn>1.0e-20);
            double coef = 8/hmn;
            double x0 = xmn;
            if (verbosity > 2)
            cout << "  --Update: periodic " << xmn << " " << xmx << " " << " h=" << hmn << endl;
            ffassert(coef>1e-10 && (xmx-xmn)*coef < 1.e7 );
             
           //  map construction ----
           for (int i1=0;i1<n1;i1++)
             {
              int ie=pke1[i1];
              const BoundaryEdge & e =Th.bedges[pke1[i1]];
              if (e.lab==label1)
                 for (int ne=0;ne<2;ne++)
                  {
                   int2 i2;
                   i2[0]=ie;
                   i2[1]=-1;
                   mp->set(e[ne].x,e[ne].y);
                   double xx=GetAny<double>((*periodic[k+kk1])(stack));
                   int i0= (int) ((xx-x0)*coef);
                   map<int,int2>::iterator im=closeto(m,i0);
                   if (im==m.end())
                    {
                     if (verbosity >50)
                      cout << xx << " " << i0 << " " << ie << endl;
                     im=m.insert(make_pair<int,int2>(i0,i2)).first;
                    }
                   else {
                     if (verbosity >50)
                     cout << xx << " " << i0 << " " << ie << " :  " << im->second[0] << " " << im->second[1] << endl;
                    assert( (im->second[1] < 0) && (im->second[0] >=0) );
                    im->second[1]=ie;}
                   
               }                                
             }
            
            for (int i2=0;i2<n2;i2++)
             {
              int ie2=pke2[i2];
              const BoundaryEdge & e =Th.bedges[ie2];
              if (e.lab==label2)
                {
		if (verbosity >50)
                    cout << i2 << " : " <<Th(e[0]) << " " << Th(e[1]) << ":: ";
                 mp->set(e[0].x,e[0].y);
                 double xx0=GetAny<double>((*periodic[k+kk2])(stack));
                 mp->set(e[1].x,e[1].y);
                 double xx1=GetAny<double>((*periodic[k+kk2])(stack));
		 if(verbosity>5 && !(verbosity >50))
		      cout << "lab2:  e[" << pke2[i2] << "]  v0:   " <<  e[0].x << " " << e[0].y << "  s = " << xx0
		      << "\t v1 " <<  e[0].x << " " << e[0].y << "  s = " << xx1 << endl;
		    
                 int i0= int((xx0-x0)*coef);
                 int i1= int((xx1-x0)*coef);
                 map<int,int2>::iterator im0=closeto(m,i0);
                 map<int,int2>::iterator im1=closeto(m,i1);
                 if(im0 == m.end() || im1 == m.end() )
                 	{ 
			
			  cout << "Abscisse: s0 = "<< xx0 << " <==> s1 " << xx1  <<endl; 
			  ExecError("periodic: Sorry one vertex of edge is losted "); }
                  int ie1=-1;
                 if      (((ie1=im0->second[0])==im1->second[1]) && (ie1>=0)) ;
                 else if (((ie1=im0->second[0])==im1->second[1]) && (ie1>=0)) ;
                 else if (((ie1=im0->second[1])==im1->second[1]) && (ie1>=0)) ;
                 else if (((ie1=im0->second[1])==im1->second[0]) && (ie1>=0)) ;
                 else if (((ie1=im0->second[0])==im1->second[0]) && (ie1>=0)) ;
                 else
                  {
                   cout << ie2 << " ~ " << im0->second[0] << " " << im0->second[1] << ", " 
                                        << im1->second[0] << " " << im1->second[1] << endl;
                   ExecError("periodic: Sorry one egde is losted "); }
		 if(verbosity>50)
		   cout << " ( " << im0->second << " , " << im1->second << " ) .. ";
		 ffassert(ie1>=0 && ie1 < Th.neb );
                 const BoundaryEdge & ep =Th.bedges[ie1];
                 mp->set(ep[0].x,ep[0].y);
                 double yy0=GetAny<double>((*periodic[k+kk1])(stack));
                 mp->set(ep[1].x,ep[1].y);
                 double yy1=GetAny<double>((*periodic[k+kk1])(stack));
		if(verbosity>50)
			cout << " e0: s  "<< xx0 << " " << xx1 << "e1 s "<< yy0 << " " << yy1  ; 
                 
                 pke1[i2]=ie1*2+ ( ( (yy1-yy0) < 0)  == ( (xx1-xx0) < 0) ) ;
		    
                 if (verbosity >50)
                 cout << " \t  edge " << ie1 << " <=> " << ie2 << " " 
                      << ( ( (yy1-yy0) < 0)  == ( (xx1-xx0) < 0) ) << "; "
                      << xx0 << " " <<xx1<<" <=> "  << yy0 << " " <<yy1<< 
                      "  ::  " <<  Th(ep[0]) << " " << Th(ep[1]) << endl ;

                }
              }
              
             
          }
          
        *mp = smp;
        for (int i=0;i<Th.neb;i++)
            ndfe[i]=i;// circular link
        for (int i=0;i<Th.nv;i++)
          ndfv[i]=i;// circular link
        for (int i=0;i<nblink1;i++)
         {
           int ie1=link1[i]/2;
           int sens = link1[i]%2;
           int ie2=link2[i];
           assert(ie1!=ie2);
	   if(!InCircularList(ndfe,ie1,ie2))   // merge of two list 
	      Exchange(ndfe[ie1],ndfe[ie2]);
	   for (int ke2=0;ke2<2;ke2++)
             {
                int ke1=ke2;
                if(!sens) ke1=1-ke1; 
                int iv1=Th(Th.bedges[ie1][ke1]);
                int iv2=Th(Th.bedges[ie2][ke2]);
                if (!InCircularList(ndfv,iv1,iv2)) {  // merge of two list 
		   Exchange(ndfv[iv2],ndfv[iv1]);
                   if (verbosity >50)
		   { 
		     cout << "  vertex " << iv1 <<  "<==> " << iv2 << " list : " << iv1;
		     int i=iv1,k=0;
		     while ( (i=ndfv[i]) != iv1 && k++<10)
			 cout << ", "<< i ; 
		       cout << endl;
		   }}                  
             }
            
         } 
         // generation de numero de dlt
         
          nbdfv = numeroteclink(ndfv) ; 
          nbdfe = numeroteclink(ndfe) ; 
          if (verbosity>2) 
            cout << " -- nb df on vertices " << nbdfv << endl;
        delete [] link1;
        delete [] link2;
        return true; //new FESpace(**ppTh,*tef,nbdfv,ndfv,nbdfe,ndfe);
      }
   }
   return false;   
}


bool  v_fes::buildperiodic(Stack stack,int & nbdfv, KN<int> & ndfv,int & nbdfe, KN<int> & ndfe) { 
  return BuildPeriodic(nbcperiodic,periodic,**ppTh,stack,nbdfv,ndfv,nbdfe,ndfe);

}
#ifdef ZZZZZZZZ
FESpace * pfes_tef::update() { 
   typedef Smallvect<int,2> int2; 
    
    if (nbcperiodic ) {
       const Mesh &Th(**ppTh);
       KN<int> ndfv(Th.nv);
       KN<int> ndfe(Th.neb);
       int nbdfv,nbdfe;
        return  new FESpace(**ppTh,*tef,nbdfv,ndfv,nbdfe,ndfe);
      }
     else 
       return  new FESpace(**ppTh,*tef);
}
#endif

struct OpMake_pfes_np {
  static const int n_name_param =1;
  static basicAC_F0::name_and_type name_param[] ;
};

basicAC_F0::name_and_type  OpMake_pfes_np::name_param[]= {
  "periodic", &typeid(E_Array) 
};

template<class pfes,class Mesh,class TypeOfFE,class pfes_tefk>
struct OpMake_pfes: public OneOperator , public OpMake_pfes_np { 
  
  struct Op: public E_F0mps { 
  public:
    
    Expression eppTh;
    Expression eppfes;
    const E_Array & atef;
    int nb;
    int nbcperiodic;
    Expression *periodic;
    
    Op(Expression ppfes,Expression ppTh, const E_Array & aatef,int nbp,Expression * pr) 
      : eppTh(ppTh),eppfes(ppfes),atef(aatef),nbcperiodic(nbp),periodic(pr) {       
    }
    ~Op() { if(periodic) delete []periodic;}
    AnyType operator()(Stack s)  const {       
      Mesh ** ppTh = GetAny<Mesh  **>( (*eppTh)(s) );
      AnyType r = (*eppfes)(s) ;
      const TypeOfFE ** tef= new  const TypeOfFE * [ atef.size()];
      for (int i=0;i<atef.size();i++)
	tef[i]= GetAny<TypeOfFE *>(atef[i].eval(s));
      pfes * ppfes = GetAny<pfes *>(r);
      bool same = true;
      for (int i=1;i<atef.size();i++)
	same &= atef[i].LeftValue() == atef[1].LeftValue();
      *ppfes = new pfes_tefk(ppTh,tef,atef.size(),s    ,nbcperiodic,periodic);
      (**ppfes).decrement();  //07/2008 FH
	// Add2StackOfPtr2FreeRC(s,*ppfes);  //  bug????  a verifier 06/07/2008
      //  delete [] tef;
      return r;}
  } ;
  
  E_F0 * code(const basicAC_F0 & args)  const
  { 
    int nbcperiodic=0;
    Expression *periodic=0;
    Expression nargs[n_name_param];
    
    args.SetNameParam(n_name_param,name_param,nargs);
    GetPeriodic(nargs[0],nbcperiodic,periodic);
    
    aType t_tfe= atype<TypeOfFE*>();
    const E_Array * a2(dynamic_cast<const E_Array *>(args[2].LeftValue()));
    ffassert(a2);
    int N = a2->size(); ;
    if (!N) CompileError(" We wait an array of Type of Element ");
    for (int i=0;i< N; i++) 
      if ((*a2)[i].left() != t_tfe) CompileError(" We wait an array of  Type of Element ");
    //    ffassert(0);
    return  new Op(args[0],args[1],*a2,nbcperiodic,periodic);
  } 
  OpMake_pfes() : 
    OneOperator(atype<pfes*>(),atype<pfes*>(),atype<Mesh **>(),atype<E_Array>()) {}
};

inline pfes* MakePtr2(pfes * const &p,pmesh * const &  a, TypeOfFE * const & tef)
{ *p=new pfes_tef(a,tef) ;
  (**p).decrement();
  return p;}

inline pfes3* MakePtr3(pfes3 * const &p,pmesh3 * const &  a, TypeOfFE3 * const & tef)
{ *p=new pfes3_tef(a,tef) ;
  (**p).decrement();
  return p;}



class OP_MakePtr2 { public:
    class Op : public E_F0mps  { public:
	//  static int GetPeriodic(Expression  bb, Expression & b,Expression & f);
	static const int n_name_param =1;
      static basicAC_F0::name_and_type name_param[] ;
      Expression nargs[n_name_param];
      typedef pfes * R;
      typedef pfes * A;
      typedef pmesh * B;
      typedef TypeOfFE * C;
      Expression a,b,c;
      int nbcperiodic ;
      Expression *periodic;
      Op(const basicAC_F0 & args);
      
      AnyType operator()(Stack s) const  {
	A p= GetAny<A>( (*a)(s) );
	B th= GetAny<B>( (*b)(s) );
	C tef= GetAny<C>( (*c)(s) );   
	//  cout << " ----------- " << endl;  
	*p=new pfes_tef(th,tef,s,nbcperiodic,periodic) ;
	(**p).decrement();
	return  SetAny<R>(p);
      } 
    }; // end Op class 
  
  typedef Op::R Result;
  static  E_F0 * f(const basicAC_F0 & args) { return  new Op(args);}
  static ArrayOfaType  typeargs() {
    return ArrayOfaType(       
			atype<Op::A>(),
			atype<Op::B>(),
			atype<Op::C>(),false ) ;}
};