fespace-v0.cpp

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

    delete [] FirstDfOfNode;
  }
 else counter--;
}

 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++;
 }

ConstructDataFElement::ConstructDataFElement (const Mesh &Th,int NbDfOnSommet,int NbDfOnEdge,int NbDfOnElement,const TypeOfMortar *tm,
int nbdfv,const int *ndfv,int nbdfe,const int *ndfe)
: counter(&thecounter),thecounter(0) 
{ 
 Make(Th,NbDfOnSommet,NbDfOnEdge,NbDfOnElement, tm,nbdfv,ndfv,nbdfe,ndfe);
}

ConstructDataFElement::ConstructDataFElement(const FESpace ** l,int k) 
: thecounter(0),counter(&thecounter) 
{
 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;
     throwassert( &Th== &l[i]->Th); 
     throwassert( l[i]->TFE.constant());
   }
   
 Make(Th,NbDfOnSommet,NbDfOnEdge,NbDfOnElement,0);   
}
 

void ConstructDataFElement::Make(const Mesh &Th,int NbDfOnSommet,int NbDfOnEdge,int NbDfOnElement,const TypeOfMortar *tm,
int nb_dfv,const int *ndfv,int nb_dfe,const int *ndfe)   
{
  *counter=0;
  Nproduit =1;
  int ndf=0,samendf=1;
  
  int nbdfe=3*NbDfOnSommet+3*NbDfOnEdge+NbDfOnElement;
  int nbne = 0;
  int nn=0;  
  int firstmul=0;
  throwassert( tm || Th.NbMortars==0);
  NbOfElements = Th.nt; //  by default
  // if mortars 
  //
  int NbOfNodeL=0;
  NbOfElements += Th.NbMortars;  
  FirstDfOfNode =0;
  FirstNodeOfElement=0;
  MaxNbDFPerElement=3*NbDfOnSommet+3*NbDfOnEdge+NbDfOnElement;

  int ks=0,ke=0,kt=0;
  if(NbDfOnSommet) { nbne+=3;
     ks=1;
     ndf=NbDfOnSommet;}

  if(NbDfOnEdge) {  nbne+=3;
     ke=1;
     samendf &= !ndf || ndf == NbDfOnEdge;
     ndf=NbDfOnEdge;}
     
  if(NbDfOnElement) {  nbne+=1;
     kt=1;
     samendf &= !ndf || ndf == NbDfOnElement;
     ndf=NbDfOnElement;}

  int NbDFonNode[7],NodeIsOn[7];
   {
     int j=0,k=0;
     if(ks) { NbDFonNode[j++]=NbDfOnSommet; NbDFonNode[j++]=NbDfOnSommet; NbDFonNode[j++]=NbDfOnSommet;}
     if(ke) { NbDFonNode[j++]=NbDfOnEdge; NbDFonNode[j++]=NbDfOnEdge; NbDFonNode[j++]=NbDfOnEdge;}
     if(kt) { NbDFonNode[j++]=NbDfOnElement;}

     if (ks) {NodeIsOn[k++]=0;NodeIsOn[k++]=1;NodeIsOn[k++]=2;}
     if (ke) {NodeIsOn[k++]=3;NodeIsOn[k++]=4;NodeIsOn[k++]=5;}
     if (kt) {NodeIsOn[k++]=6;}
     
     throwassert(j == nbne);
  }
     
    MaxNbNodePerElement=nbne;

//  
   if ( ks && (!ke && ! kt) && (ndfv==0 && ndfe==0))
     {nn=Th.nv;
      NodesOfElement=0;
      }
   else {
//    constuction du tableau  NodesOfElement  bofbof 

//     computation of the length lne of array NodesOfElement        
       int lne=  Th.nt*nbne;
       if (Th.NbMortars)
        {
          samendf= false;
          NbOfNodeL=Th.NbMortars;
          throwassert(tm);
          FirstNodeOfElement = new int[NbOfElements+1];
          int k=0,kk=0;          
          for (k=0;k<Th.nt;k++,kk+=nbne)
            FirstNodeOfElement[k] = kk;
           
          for (int im=0;im<Th.NbMortars;im++) // (www) code 
            {
             FirstNodeOfElement[k++]=lne;           
             lne += tm->NbOfNodes(Th,Th.mortars[im]);
            }
          FirstNodeOfElement[k++]=lne;           
        }       
        
       NodesOfElement = new int[lne];
       
       for (int i=0;i<lne;i++) 
         NodesOfElement[i]=-1;
       int i=0;
       int oe=0; 
       if(ks) { oe=3;nn=ndfv ? nb_dfv : Th.nv;}
       
       if (ke && ndfe) { 
        for (int be=0;be<Th.neb;be++)
         {
          int j,k=Th.BoundaryElement(be,j);
          int jj=j;
          int kk=Th.ElementAdj(k,jj);
          if ( kk >=0 && jj>=0)   
            NodesOfElement[kk*nbne+oe+jj] = nn + ndfe[be]   ; // adj
          NodesOfElement[k*nbne+oe+j]   = nn + ndfe[be]   ; // new         
         }
         nn += nb_dfe;
        }
       for (int k=0;k<Th.nt;k++)
        {
          int iold=i;
          if(ks) {
            for (int j=0;j<3;j++)
             NodesOfElement[i++]= ndfv ?  ndfv[Th(k,j)] : Th(k,j) ;
          }
          if(ke) {
            for (int j=0;j<3;j++)
             if (NodesOfElement[i]<0) {
               int jj=j;
               int kk=Th.ElementAdj(k,jj);
               
               NodesOfElement[kk*nbne+oe+jj] = nn   ; // adj   
               NodesOfElement[i++]           = nn++ ; // new
             }  else i++;
          }
          if(kt) 
             NodesOfElement[i++]= nn++;
         // cout << k ;
         // dump(" ",i-iold, NodesOfElement+iold);
        }
       // cout << i << " " << Th.nt*nbne << endl;
        firstmul=nn;
        if (Th.NbMortars)
          {  
            //  construction of the mortars element 
           int * color= new int [firstmul]; 
           //  
           int thecolor=0;
           for (int j=0;j<firstmul;j++) 
             color[j]=thecolor;
             
           for (int im=0;im<Th.NbMortars;im++)
            {   
             int iold=i;         
             thecolor++; // get a new color
             // tm->ConstructionOfNode(Th,im,NodesOfElement,FirstNodeOfElement,nn);  
             Mortar & M(Th.mortars[im]);
              NodesOfElement[i++] = nn++; // the first node is the lag. mul.
             int n=M.NbT();
             for (int kk=0;kk<n;kk++)
             {
              int K,e;
               K=M.T_e(kk,e); 
               int kb=FirstNodeOfElement[K];
               int ke=FirstNodeOfElement[K+1];             
               for (int j=kb,jj=0;j<ke;j++,jj++)
                if (onWhatIsEdge[e][NodeIsOn[jj]])
                 { // the node jj is on edge e
                    int node=NodesOfElement[j];
                   // cout << "." << jj << " K=" << K <<" "<<e << " n=" << node ;
                    throwassert(node<firstmul);
                    if (color[node] != thecolor) //  new node => adding
                     { 
                    //   cout << "a, ";
                       color[node] = thecolor;
                       NodesOfElement[i++] = node;
                     }
                   // else cout << ", ";
                 }
               }
               //cout << endl;
               
               //cout << im ;
               //dump(": ",i-iold, NodesOfElement+iold);
               throwassert(i==FirstNodeOfElement[im+Th.nt+1]);
              }
             delete [] color;
          } 
        else
          throwassert(i==Th.nt*nbne);
        NbOfNode=nn;
                
        
   }
  NbOfNode=nn;  
  int NbOfDFL=0;  
   if (! samendf) 
       { throwassert(NodesOfElement);
         FirstDfOfNode= new int [nn+1];
         for (int i=0;i<=nn;i++) FirstDfOfNode[i]=-1; 
         int i=0;
         //  the classical part (FEM)
         for (int k=0;k<Th.nt;k++)
           for (int j=0;j<nbne;j++) // thanks to student
             FirstDfOfNode[ NodesOfElement[i++]+1]=NbDFonNode[j];
         //  the mortars parts juste the mulplicator 
         
         for (int km=0,k=Th.nt;km<Th.NbMortars;km++,k++) 
            {  //  the lag. mult. is the first node of the mortar --
              throwassert(FirstNodeOfElement);
              //  hack 
              int fk=FirstNodeOfElement[k];
              int lk=FirstNodeOfElement[k+1];
              int ndlmul = tm->NbLagrangeMult(Th,Th.mortars[km]);  //  On node par 
              int nodemul = NodesOfElement[fk]; // the first node is the lagr. mul.
              throwassert(FirstDfOfNode[nodemul+1]==-1);      
              FirstDfOfNode[nodemul+1]= ndlmul;
              NbOfDFL += ndlmul;
              int nbdle=0;
              int nbnm=lk-fk;
              for (int j=fk;j<lk;j++)
               nbdle+=FirstDfOfNode[NodesOfElement[j]+1];
              MaxNbDFPerElement = Max(MaxNbDFPerElement,nbdle);
            }
             
         FirstDfOfNode[0]=0;
         for (int i=0;i<=nn;i++) throwassert(FirstDfOfNode[i]!=-1); 

         for (int i=0;i<nn;i++) 
              FirstDfOfNode[i+1] += FirstDfOfNode[i] ;
           NbOfDF=  FirstDfOfNode[nn];       
        }
    else
       {
         NbOfDF = nn*ndf; 
         Nproduit = ndf;        
       }
   MaxNbNodePerElement=nbne;
   
   cout << " Nb Of Nodes = " << nn << endl;   
   if(NbOfNodeL)        
     cout << " Nb of Lagrange Mul Node = " << NbOfNodeL << endl;        
   cout << " Nb of DF = " << NbOfDF << endl;   
   if(NbOfDFL) {  
     cout << " Nb of Lagrange Mul DF = "   << NbOfDFL << endl;  
     cout << " MaxNbDFPerElement     =   " << MaxNbDFPerElement << endl;
   };

}

FESpace::FESpace(const FESpace & Vh,int k )
 :
     ptrTFE(new TypeOfFEProduit(k,*Vh.TFE[0])),
     TFE(1,0,ptrTFE),
     cmesh(Vh.Th),
     cdef(Vh.cdef?new ConstructDataFElement(Vh.cdef,k):0),
     N(Vh.N*k),
     Nproduit(Vh.Nproduit*k),
     nb_sub_fem(TFE[0]->nb_sub_fem),
     dim_which_sub_fem(TFE[0]->dim_which_sub_fem),
     
     Th(Vh.Th),
     NbOfDF(Vh.NbOfDF*k),
     NbOfElements(Vh.NbOfElements),
     NbOfNodes(Vh.NbOfNodes),
     MaxNbNodePerElement(Vh.MaxNbNodePerElement),
     MaxNbDFPerElement(Vh.MaxNbDFPerElement*k),
     NodesOfElement(Vh.NodesOfElement),
     FirstDfOfNodeData(cdef?cdef->FirstDfOfNode:0),
     FirstNodeOfElement(Vh.FirstNodeOfElement),
     tom(0) {
     if(cdef) renum();}
 
FESpace::FESpace(const FESpace ** Vh,int k )
 :
     ptrTFE(new TypeOfFESum(Vh,k)),
     TFE(1,0,ptrTFE),
     cmesh((**Vh).Th),
     cdef(new ConstructDataFElement(Vh,k)),
     N(sum(Vh,&FESpace::N,k)),
     Nproduit(cdef->Nproduit),     
     nb_sub_fem(TFE[0]->nb_sub_fem),
     dim_which_sub_fem(TFE[0]->dim_which_sub_fem),

     Th((**Vh).Th),
     NbOfDF(cdef->NbOfDF),
     NbOfElements(cdef->NbOfElements),
     NbOfNodes(cdef->NbOfNode),
     MaxNbNodePerElement(cdef->MaxNbNodePerElement),
     MaxNbDFPerElement(cdef->MaxNbDFPerElement),
     NodesOfElement(cdef->NodesOfElement),
     FirstDfOfNodeData(cdef->FirstDfOfNode),
     FirstNodeOfElement(cdef->FirstNodeOfElement),
     tom(0) {
     if(cdef) renum(); }
     
FESpace::FESpace(const Mesh & TTh,const TypeOfFE ** tef,int k,int nbdfv,const int *ndfv,int nbdfe,const int *ndfe )
 :
     ptrTFE(new TypeOfFESum(tef,k)),
     TFE(1,0,ptrTFE),
     cmesh(TTh),
     cdef(new ConstructDataFElement(TTh,sum(tef,&TypeOfFE::NbDfOnVertex,k),
                                        sum(tef,&TypeOfFE::NbDfOnEdge,k),
                                        sum(tef,&TypeOfFE::NbDfOnElement,k),
                                        0,nbdfv,ndfv,nbdfe,ndfe)),
     N(sum(tef,&TypeOfFE::N,k)),
     Nproduit(cdef->Nproduit),     
     nb_sub_fem(TFE[0]->nb_sub_fem),
     dim_which_sub_fem(TFE[0]->dim_which_sub_fem),

     Th(TTh),
     NbOfDF(cdef->NbOfDF),
     NbOfElements(cdef->NbOfElements),
     NbOfNodes(cdef->NbOfNode),
     MaxNbNodePerElement(cdef->MaxNbNodePerElement),
     MaxNbDFPerElement(cdef->MaxNbDFPerElement),
     NodesOfElement(cdef->NodesOfElement),
     FirstDfOfNodeData(cdef->FirstDfOfNode),
     FirstNodeOfElement(cdef->FirstNodeOfElement),
     tom(0) {
     if(cdef) renum(); }


 FESpace::FESpace(const Mesh & TTh,const TypeOfFE & tef,int nbdfv,const int *ndfv,int nbdfe,const int *ndfe)
   :  
     ptrTFE(0),
     TFE(1,0,&tef),
     cmesh(TTh),
     cdef(new ConstructDataFElement(TTh,tef.NbDfOnVertex,tef.NbDfOnEdge,tef.NbDfOnElement,0,nbdfv,ndfv,nbdfe,ndfe)),
     N(tef.N),
     Nproduit(cdef->Nproduit),
     nb_sub_fem(TFE[0]->nb_sub_fem),
     dim_which_sub_fem(TFE[0]->dim_which_sub_fem),
     Th(TTh),
     NbOfDF(cdef->NbOfDF),
     NbOfElements(cdef->NbOfElements),
     NbOfNodes(cdef->NbOfNode),
     MaxNbNodePerElement(cdef->MaxNbNodePerElement),
     MaxNbDFPerElement(cdef->MaxNbDFPerElement),
     NodesOfElement(cdef->NodesOfElement),
     FirstDfOfNodeData(cdef->FirstDfOfNode),
     FirstNodeOfElement(cdef->FirstNodeOfElement),
     tom(0) {
     if(tef.NbDfOnVertex || tef.NbDfOnEdge) renum();
     }
     
 FESpace::~FESpace()
   {
     SHOWVERB(cout << " FESpace::~FESpace() " << endl);
      delete  cdef;
      if(ptrTFE) 
        delete  ptrTFE;
   }

 FESpace::FESpace(const Mesh & TTh,const TypeOfFE & tef,const TypeOfMortar & tm)
   :  
     ptrTFE(0),
     TFE(1,0,&tef),
     cmesh(TTh),
     cdef(new ConstructDataFElement(TTh,tef.NbDfOnVertex,tef.NbDfOnEdge,tef.NbDfOnElement,&tm)),
     N(tef.N),
     Nproduit(1),
     nb_sub_fem(TFE[0]->nb_sub_fem),
     dim_which_sub_fem(TFE[0]->dim_which_sub_fem),     
     Th(TTh),
     NbOfDF(cdef->NbOfDF),
     NbOfElements(cdef->NbOfElements),
     NbOfNodes(cdef->NbOfNode),
     MaxNbNodePerElement(cdef->MaxNbNodePerElement),
     MaxNbDFPerElement(cdef->MaxNbDFPerElement),
     NodesOfElement(cdef->NodesOfElement),
     FirstDfOfNodeData(cdef->FirstDfOfNode),
     FirstNodeOfElement(cdef->FirstNodeOfElement),
     tom(&tm) { 
     // cout << "avant renum ="<< *this <<endl;
       renum();
     // cout << "apres renum ="<< *this <<endl;
     }
     
void ConstructDataFElement::renum(const long *r,int l)   
 { 
/*   cout << "renu=" << l << ":" << endl;
   for (int i=0;i<NbOfNode;i++)
      if (i%10) cout << r[i] << "\t";
      else cout << "\n " << i << ":\t" << r[i] << "\t";
      cout << endl; 
*/
   throwassert(this);
   if (NodesOfElement) 
     for (int i=0;i< l ; i++)
       NodesOfElement[i]=r[NodesOfElement[i]];
   if(FirstDfOfNode)  
    { int k,i,*n=new int[NbOfNode];
      for ( i=0;i<NbOfNode;i++)        
         n[r[i]]=FirstDfOfNode[i+1]-FirstDfOfNode[i]; 
      FirstDfOfNode[0]=k=0;
      for(i=0;i< NbOfNode;)
        {k+=n[i];
         FirstDfOfNode[++i]=k;}
       delete [] n;