fespace.cpp

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

  data1 = data + n*5+N; // april 2006  add 2 array ????
  
  int c=0;
  //  int ki= 0; 
// recherche des noeuds
   KN<int> w(7),nn(7);
   w=0;
   nn=0; 
   
         
   for ( j=0;j<kk;j++)
     for ( i=0;i<teb[j]->NbDoF;i++)
         nn[teb[j]->DFOnWhat[i]]++;
   nbn=0;      
   for( j=0;j<7;j++)
     if (nn[j]) nn[j]=nbn++;
     else nn[j]=-1;
   KN<int> dln(7);
   dln=0;
  // nn donne numero de noeud sur what            
   for ( j=0;j<kk;j++)
     for ( i=0;i<teb[j]->NbDoF;i++)
       data[c++] = teb[j]->DFOnWhat[i];
    
   for ( j=0;j<kk;j++)
    {
     int  cc=c;
     for ( i=0;i<teb[j]->NbDoF;i++)
       data[c++] = teb[j]->DFOfNode[i]+dln[teb[j]->DFOnWhat[i]];
     for ( i=0;i<teb[j]->NbDoF;i++)
       dln[teb[j]->DFOnWhat[i]]=Max(dln[teb[j]->DFOnWhat[i]],data[cc++]+1);      
    }
        
       
   for ( j=0;j<kk;j++)
    { 
     //  w renumerotation des noeuds 
     //  Ok si un noeud par what 
     for ( i=0;i<teb[j]->NbDoF;i++)
       data[c++] = nn[teb[j]->DFOnWhat[i]];
    }
    
   for ( j=0;j<kk;j++)
     for ( i=0;i<teb[j]->NbDoF;i++)
       data[c++] = j; //  node from of FE

    
   for ( j=0;j<kk;j++)
     for ( i=0;i<teb[j]->NbDoF;i++)
       data[c++] = i; //  node from of df in FE
  // error -- here 
  //in case of [P2,P2],P1  
   // we expect 0,0,1   and we get 0 1 2 
   // => wrong BC ???? 
   int xx=0;
   for (j=0;j<kk;j++)
     { 
      int xxx=xx;
      for (i=0;i<teb[j]->N;i++)
       { 
         data[c] = teb[j]->dim_which_sub_fem[i]+xx;
         xxx=Max(xxx,data[c]+1);
         c++;
       }
       xx=xxx;
     }
    

  //  ou dans la partie miminal element finite atomic 
 
   int ci=n;
   int cf=2*n;
   int cl=cf+N;;
   int cj=0;
   int ccc=0;
   for ( j=0;j<kk;ccc+=teb[j++]->nb_sub_fem)
     for ( i=0;i<teb[j]->NbDoF;i++)
       {
	 int il= teb[j]->fromASubDF[i];
	 int jl= teb[j]->fromASubFE[i];
	 data1[ci++]=il;
	 data1[cj++]=ccc+jl;      
       }

   for (int  j=0,ccn=0 ; j<kk ; ccn += teb[j++]->NbDoF)     
     for(int k=0;k<teb[j]->N;++k)
       {
	 data1[cf++] = ccn + teb[j]->begin_dfcomp[k];
	 data1[cl++] = ccn + teb[j]->end_dfcomp[k];
       }
   ffassert(cl==2*n+2*N);    
  
  ffassert(c== 5*n+N);      
/*  int cc=0;
   cout << " Data : " << endl;
  for ( i=0;i<5;i++)    {
    for (j=0;j<n;j++)
      cout << " " << data[cc++];
     cout << endl;}
 cout << " which " ;
 for (i=0;i<N;i++)
   cout << " " << data[cc++];
  cout << endl;*/
}

class TypeOfFE_P1Lagrange : public  TypeOfFE { public:  
  static int Data[];
  static double Pi_h_coef[];
   TypeOfFE_P1Lagrange(): TypeOfFE(1,0,0,1,Data,1,1,3,3,Pi_h_coef)
    { const R2 Pt[] = { R2(0,0), R2(1,0), R2(0,1) }; 
      for (int i=0;i<NbDoF;i++) {
       pij_alpha[i]= IPJ(i,i,0);
       P_Pi_h[i]=Pt[i]; }
     }
  // void FB(const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
   void FB(const bool * whatd,const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
   
//   void D2_FB(const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
  // void Pi_h(const baseFElement & K,RN_ & val, InterpolFunction f, R* v,int, void *) const;
virtual R operator()(const FElement & K,const  R2 & PHat,const KN_<R> & u,int componante,int op) const ;
   
} ;
///////////////////////////////////////////////////////////////////////////////
// FH pour tester des idee de schema ----  Juin 2005 ---
///////////////////////////////////////////////////////////////////////////////

//  un VF cell centre 
class TypeOfFE_P0VF : public  TypeOfFE { public:  
  static int Data[];
  static double Pi_h_coef[];
   TypeOfFE_P0VF(): TypeOfFE(1,0,0,1,Data,1,1,3,3,Pi_h_coef)
    { const R2 Pt[] = { R2(0,0), R2(1,0), R2(0,1) }; 
      for (int i=0;i<NbDoF;i++) {
       pij_alpha[i]= IPJ(i,i,0);
       P_Pi_h[i]=Pt[i]; }
     }
   void FB(const bool * whatd,const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
   virtual R operator()(const FElement & K,const  R2 & PHat,const KN_<R> & u,int componante,int op) const ;
   
} ;
int TypeOfFE_P0VF::Data[]={0,1,2,       0,0,0,       0,1,2,       0,0,0,        0,1,2,       0, 0,3};
double TypeOfFE_P0VF::Pi_h_coef[]={1.,1.,1.}; //  bofbof a verifier ...

 R TypeOfFE_P0VF::operator()(const FElement & K,const  R2 & PHat,const KN_<R> & u,int componante,int op) const 
{ 
   R u0(u(K(0))), u1(u(K(1))), u2(u(K(2)));
   R r=0;
   if (op==0)
    {
      R l0=0,l1=PHat.x,l2=PHat.y;
      l1 = l1 * 3. < 1;
      l2 = l2 * 3. < 1;
      l0 = 1 - l0 -l2;
      
      r = u0*l0+u1*l1+l2*u2;
    }
   else
    { 
      r =0; 
    }
 //  cout << r << "\t";
   return r;
}


void TypeOfFE_P0VF::FB(const bool *whatd,const Mesh & ,const Triangle & K,const R2 & P,RNMK_ & val) const
{
//  const Triangle & K(FE.T);
  if (whatd[op_id]) 
   {
  R2 A(K[0]), B(K[1]),C(K[2]);
  R l0=1-P.x-P.y,l1=P.x,l2=P.y; 
  l1 = l1 * 3. < 1;
  l2 = l2 * 3. < 1;
  l0 = 1 - l0 -l2;
  
  if (val.N() <3) 
   throwassert(val.N() >=3);
  throwassert(val.M()==1 );
//  throwassert(val.K()==3 );
  
  val=0; 
  RN_ f0(val('.',0,op_id)); 
  
    f0[0] = l0;
    f0[1] = l1;
    f0[2] = l2;}
 
}


///////////////////////////////////////////////////////////////////////////////
////////////////////////////////// NEW ////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

class TypeOfFE_P1Bubble : public  TypeOfFE { public:  
  static int Data[];
  static double Pi_h_coef[];
   TypeOfFE_P1Bubble(): TypeOfFE(1,0,1,1,Data,1,1,4,4,Pi_h_coef)
    { const R2 Pt[] = { R2(0,0), R2(1,0), R2(0,1), R2(1./3.,1./3.) }; 
      for (int i=0;i<NbDoF;i++) {
       pij_alpha[i]= IPJ(i,i,0);
       P_Pi_h[i]=Pt[i]; }
     }
  // void FB(const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
   void FB(const bool * whatd,const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
   
//   void D2_FB(const Me禨eria哪哪膕h & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
 //  void Pi_h(const baseFElement & K,RN_ & val, InterpolFunction f, R* v,int, void *) const;
  //virtual R operator()(const FElement & K,const  R2 & PHat,const KN_<R> & u,int componante,int op) const ;
   
} ;
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

class TypeOfFE_P2Lagrange : public  TypeOfFE { public:  
  static int Data[];
  static double Pi_h_coef[];
  
   TypeOfFE_P2Lagrange(): TypeOfFE(1,1,0,1,Data,3,1,6,6,Pi_h_coef)
    { const R2 Pt[] = { R2(0,0), R2(1,0), R2(0,1),R2(0.5,0.5),R2(0,0.5),R2(0.5,0) };
      for (int i=0;i<NbDoF;i++) {
       pij_alpha[i]= IPJ(i,i,0);
       P_Pi_h[i]=Pt[i]; }
     }
   
  // void FB(const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
   void FB(const bool * whatd,const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
 //  void D2_FB(const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
  // void Pi_h(const baseFElement & K,RN_ & val, InterpolFunction f, R* v,int, void *) const;
} ;


class TypeOfFE_P2bLagrange : public  TypeOfFE { public:  
  static int Data[];
  static double Pi_h_coef[];
  
   TypeOfFE_P2bLagrange(): TypeOfFE(1,1,1,1,Data,3,1,7,7,Pi_h_coef)
    { const R2 Pt[] = { R2(0,0), R2(1,0), R2(0,1),R2(0.5,0.5),R2(0,0.5),R2(0.5,0), R2(1./3.,1./3.) };
      for (int i=0;i<NbDoF;i++) {
       pij_alpha[i]= IPJ(i,i,0);
       P_Pi_h[i]=Pt[i]; }
     }
   
  // void FB(const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
   void FB(const bool * whatd,const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
 //  void D2_FB(const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
  // void Pi_h(const baseFElement & K,RN_ & val, InterpolFunction f, R* v,int, void *) const;
} ;

int TypeOfFE_P1Lagrange::Data[]={0,1,2,       0,0,0,       0,1,2,       0,0,0,        0,1,2,       0, 0,3};
int TypeOfFE_P1Bubble::Data[]={0,1,2,6,     0,0,0,0,     0,1,2,3,     0,0,0,0,        0,1,2,3,     0, 0,4};
int TypeOfFE_P2Lagrange::Data[]={0,1,2,3,4,5, 0,0,0,0,0,0, 0,1,2,3,4,5, 0,0,0,0,0,0,  0,1,2,3,4,5, 0 ,0,6};
int TypeOfFE_P2bLagrange::Data[]={0,1,2,3,4,5,6, 0,0,0,0,0,0,0, 0,1,2,3,4,5,6, 0,0,0,0,0,0,0,  0,1,2,3,4,5,6, 0,0,7};
double TypeOfFE_P1Lagrange::Pi_h_coef[]={1.,1.,1.};
double TypeOfFE_P1Bubble::Pi_h_coef[]={1.,1.,1.,1.};
double TypeOfFE_P2Lagrange::Pi_h_coef[]={1.,1.,1.,1.,1.,1.};
double TypeOfFE_P2bLagrange::Pi_h_coef[]={1.,1.,1.,1.,1.,1.,1.};

inline void dump(char *m,int n,int * p)
{
  cout << m ;
  for (int i=0;i<n;i++) cout << " " << p[i] ;
  cout << endl;
}



ConstructDataFElement::~ConstructDataFElement()
{
  if((*counter)--==0) 
   {
   // cout << " delete ConstructDataFElement " <<   NodesOfElement << " " <<  FirstNodeOfElement << " "<< FirstDfOfNode << "  " << counter << endl;
    delete [] NodesOfElement;
    delete []  FirstNodeOfElement;
    delete [] FirstDfOfNode;
    (*counter)--; // correct bug oct 2008
    delete counter;
  }
//  else (*counter)--; // correct bug oct 2008
  // else 
     // cout << " no delete ConstructDataFElement " <<   NodesOfElement << " " <<  FirstNodeOfElement << " "<< FirstDfOfNode << "  " << counter << endl;
// (*counter)--; // correction mai 2006 bug in counter incrementation 
}

 ConstructDataFElement::ConstructDataFElement(const ConstructDataFElement * t,int k)
  ://thecounter(0), 
   counter(t->counter),
   MaxNbNodePerElement(t->MaxNbNodePerElement),
   MaxNbDFPerElement(t->MaxNbDFPerElement*k),
   NodesOfElement(t->NodesOfElement),
   FirstNodeOfElement(t->FirstNodeOfElement),
   FirstDfOfNode(0),
   NbOfElements(t->NbOfElements),
   NbOfDF(t->NbOfDF*k),
   NbOfNode(t->NbOfNode),
   Nproduit(t->Nproduit*k)
 {
   throwassert(t==0 || t->FirstDfOfNode==0);
   (*counter)++;      // correction mai 2006 bug in counter incrementation  

 }

ConstructDataFElement::ConstructDataFElement (const Mesh &Th,/*int NbDfOnSommet,int NbDfOnEdge,int NbDfOnElement*/
const  KN<const TypeOfFE *> & TFEs,const TypeOfMortar *tm,
int nbdfv,const int *ndfv,int nbdfe,const int *ndfe)
: counter(NewCounter())
{ 
 Make(Th,TFEs,/*NbDfOnSommet,NbDfOnEdge,NbDfOnElement,*/ tm,nbdfv,ndfv,nbdfe,ndfe);
}

ConstructDataFElement::ConstructDataFElement(const FESpace ** l,int k,const KN<const TypeOfFE *>  & TFEs) 
: counter(NewCounter()) 
{
 int NbDfOnSommet=0;
 int NbDfOnEdge=0;
 int NbDfOnElement=0;
 const Mesh & Th(l[0]->Th);
 for  (int i=0;i<k;i++)
   {
     NbDfOnSommet += l[i]->TFE[0]->NbDfOnVertex;
     NbDfOnEdge += l[i]->TFE[0]->NbDfOnEdge;
     NbDfOnElement += l[i]->TFE[0]->NbDfOnElement;
     ffassert( &Th== &l[i]->Th); 
     ffassert( l[i]->TFE.constant());
   }
   
 Make(Th,TFEs);//NbDfOnSommet,NbDfOnEdge,NbDfOnElement,0);   
}
 

void ConstructDataFElement::Make(const Mesh &Th,
const  KN<const TypeOfFE *> & TFEs,

/*int NbDfOnSommet,int NbDfOnEdge,int NbDfOnElement,*/const TypeOfMortar *tm,
int nb_dfv,const int *ndfv,int nb_dfe,const int *ndfe) 
/* pour le condition de periodicit