fespace.cpp

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

// -*- Mode : c++ -*-
//
// SUMMARY  :      
// USAGE    :        
// ORG      : 
// AUTHOR   : Frederic Hecht
// E-MAIL   : [e-mail]
//

/*
 
 This file is part of Freefem++
 
 Freefem++ is free software; you can redistribute it and/or modify
 it under the terms of the GNU Lesser General Public License as published by
 the Free Software Foundation; either version 2.1 of the License, or
 (at your option) any later version.
 
 Freefem++  is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU Lesser General Public License for more details.
 
 You should have received a copy of the GNU Lesser General Public License
 along with Freefem++; if not, write to the Free Software
 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 */

#include <cmath>
#include <cstdlib>
#include "error.hpp"
#include <iostream>
#include <fstream>
#include <map>
#include "rgraph.hpp"
using namespace std;  

#include "RNM.hpp"
#include "fem.hpp"
#include "FESpacen.hpp"
#include "FESpace.hpp"
extern long verbosity ;
namespace  Fem2D {

 int  Make(const TypeOfFE ** t,int k,KN<R2> & P,KN<int> & I)
 {
   typedef  TypeOfFE::IPJ IPJ;

   int n=0,nn=0;

   for (int i=0;i<k;i++)
    {
     const KN<R2> p(t[i]->P_Pi_h);
     for (int j=0;j<p.N();j++,nn++)
     {
       P[n]=p[j]; // par defaut un nouveau => ajout
       I[nn]=n++;       
       for (int l=0;l<n-1;l++)
        {
          R2 QP(p[j],P[l]);
          if ( (QP,QP) < 1.0e-12 ) { 
            I[nn]=l;
            n--; // on a retrouver =>  detruit
            break;}
        }
       
     }
    }
   return n; // nombre de point commun
 }
 
 KN<TypeOfFE::IPJ > Makepij_alpha(const TypeOfFE ** t,int k)
 {
 // Attention les df est numerote de facon croissant 
 // en faisant une boucle sur les TypeOfFE
 // comme dans la class TypeOfFESum
  typedef TypeOfFE::IPJ IPJ;
  int n=0,m=0;
  for (int i=0;i< k;i++) {
    n += t[i]->pij_alpha.N();
    m += t[i]->P_Pi_h.N();
    }
  KN<TypeOfFE::IPJ > ij(n);
  KN<int> I(m);
  KN<R2> P(m);
  Make(t,k,P,I);
  int p0=0,i0=0,N0=0,nn=0;
  for (int i=0;i< k;i++) {
    const KN<IPJ > p(t[i]->pij_alpha);
     for (int j=0;j<p.N();j++,nn++)  
      {
        ij[nn].i= p[j].i+i0; //  comme dans TypeOfFESum
        ij[nn].j= N0+ p[j].j;
        ij[nn].p= I[p[j].p+p0];
       
      //   cout << nn << "Makepij_alpha: " << ij[nn].i << " " << ij[nn].p << " " << ij[nn].j << endl;
        
      }
    i0+=t[i]->NbDoF;
    p0+=t[i]->P_Pi_h.N();
    N0+=t[i]->N;}
  return ij; 
 }
 KN<R2 > MakeP_Pi_h(const TypeOfFE **t,int k)
 {
  int np=0;
  for (int i=0;i< k;i++) 
    np += t[i]->P_Pi_h.N();
    
  KN< R2 >  yy(np);
  KN<int> zz(np);
  int kk=Make(t,k,yy,zz);
 // cout << " MakeP_Pi_h: " << kk << " from " << np << endl; 
  return yy(SubArray(kk));
 
 }

ListOfTFE * ListOfTFE::all ; // list of all object of this type 

void init_static_FE(); //   to correct so probleme with static Library FH aout 2004
//  the list of other FE file to force the link 

ListOfTFE::ListOfTFE (const char * n,TypeOfFE *t) : name(n),tfe(t) 
{
  if(!t)
  assert(t);
  static int count=0;
  if (count++==0) 
    all=0; // init of all in dependant of the ordre of the objet file   
  next=all;
  all=this;
 //  to correct so probleme with static Library FH aout 2004 
 init_static_FE();
}

const TypeOfFE ** Make(const FESpace **l,int k) {
  const TypeOfFE** p=new const TypeOfFE*[k];
  for (int i=0;i<k;i++)
    p[i]=l[i]->TFE[0];
  return p;
}
const TypeOfFE ** Make(const TypeOfFE **l,int k) {
  const TypeOfFE** p=new const TypeOfFE*[k];
  for (int i=0;i<k;i++)
    p[i]=l[i];
  return p;
}

  
bool Same(const FESpace **l,int k)
{
   for (int i=1;i<k;i++)
    if (l[0] != l[i] ) return false;
   return true;
}   

class FESumConstruct { protected:
   const int k;
   const TypeOfFE ** teb;
   int nbn; // nb of node 
   int  *  data;
   int  * data1;
   int  * const NN; //  NN[ i:i+1[ dimension de l'element i
   int  * const DF; // DF[i:i+1[  df associe a l'element i
   int  * const comp; //  
   FESumConstruct(int kk,const TypeOfFE **t);   
   virtual ~FESumConstruct(){
     delete [] DF;
     delete [] NN;
     delete [] comp;
     delete [] data;}   
};

  void FElement::Pi_h(RN_ val,InterpolFunction f,R *v, void * arg=0) const {
   // routine: a  tester FH.
    FElement::aIPJ ipj(Pi_h_ipj()); 
    FElement::aR2  PtHat(Pi_h_R2()); 
    KN<R>   Aipj(ipj.N());
    KNM<R>  Vp(N,PtHat.N());
        
     Pi_h(Aipj);
     for (int p=0;p<PtHat.N();p++)
          { 
            f(v,T(PtHat[p]),*this,T.lab,PtHat[p],arg);
            KN_<double> Vpp(Vp('.',p));
            for (int j=0;j<N;j++)          
               Vpp[j]=v[j];              
           }
           
         for (int i=0;i<Aipj.N();i++)
          { 
           const FElement::IPJ &ipj_i(ipj[i]);
           val[ipj_i.i] += Aipj[i]*Vp(ipj_i.j,ipj_i.p);           
          }
 }

class TypeOfFESum: public FESumConstruct, public  TypeOfFE { public:
   TypeOfFESum(const FESpace **t,int kk): 
     FESumConstruct(kk,Make(t,kk)),TypeOfFE(teb,kk,data,data1) {}
       TypeOfFESum(const TypeOfFE **t,int kk): 
     FESumConstruct(kk,Make(t,kk)),TypeOfFE(teb,kk,data,data1) {}

  // 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 * arg ) const; 
   virtual void Pi_h_alpha(const baseFElement & K,KN_<double> & v) const
    {
      int k0=0;
      for (int i=0;i<k;i++) {
         //  int n=teb[i]->NbDoF; // ici BUG 28/11/2006 FH 
          int n=teb[i]->pij_alpha.N(); // ici BUG 
          KN_<R> sv(v(SubArray(n,k0)));
          teb[i]->Pi_h_alpha(K,sv);
          k0+= n;}
      assert(pij_alpha.N()==k0);
    }
   ~TypeOfFESum(){  delete []  teb;}
} ;

class FEProduitConstruct { protected:
   int k;
   const TypeOfFE & teb;
   int * data;
   int * data1;
   FEProduitConstruct(int kk,const TypeOfFE &t)  ;   
   ~FEProduitConstruct(){delete [] data;}   
};

class TypeOfFEProduit: protected FEProduitConstruct, public  TypeOfFE { public:  
  TypeOfFEProduit(int kk,const TypeOfFE &t): 
    FEProduitConstruct(kk,t),TypeOfFE(t,kk,data,data1)  {}
  
  
  void FB(const bool * whatd,const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
  virtual void Pi_h_alpha(const baseFElement & K,KN_<double> & v) const
  { int nbof=teb.NbDoF;
  for (int i=0,k0=0;i<k;i++,k0+=nbof)
    {          
      KN_<R> sv(v(SubArray(nbof,k0)));
      teb.Pi_h_alpha(K,sv);
    }
  }
  
  ~TypeOfFEProduit(){}
} ;

FEProduitConstruct::FEProduitConstruct(int kk,const TypeOfFE &t)
 :k(kk),teb(t) 
{  
  int m= teb.NbDoF;
  KN<int> nn(teb.NbNode);
  nn=0; // nb de dl par noeud 
  for (int i=0;i<m;i++)
    nn[teb.NodeOfDF[i]]++;   
    
  int n= m*kk;
  int No=teb.N;
  int N= No*kk;
  data = new int [n*(5+2)+3*N];
  data1 = data + n*(5)+N; // april 2006  add 2 array ????
  int c=0;
  
  for (int i=0;i<m;i++)
   for (int j=0;j<kk;j++)
    data[c++] = teb.DFOnWhat[i];
  
  for (int i=0;i<m;i++)
   for (int j=0;j<kk;j++) // num of df on node  for the df = j 
      data[c++] = teb.DFOfNode[i]+j*nn[teb.NodeOfDF[i]];
    
  for (int i=0;i<m;i++)
   for (int j=0;j<kk;j++) 
    data[c++] = teb.NodeOfDF[i]; //  node of df

  for (int i=0;i<m;i++)
   for (int j=0;j<kk;j++) 
    data[c++] = j; //  node from of FE
    
  for (int i=0;i<m;i++)
   for (int j=0;j<kk;j++) 
    data[c++] = i; //  node from of df in FE
    
   for (int j=0;j<kk;j++) 
    for (int i=0;i<teb.N;i++) 
      data[c++]= teb.dim_which_sub_fem[i] + teb.nb_sub_fem*j ;
 
   int ci=n;
   int cj=0;

  //  ou dans la partie miminal element finite atomic 
  for (int i=0;i<m;i++)
   for (int j=0;j<kk;j++) 
    {
      int il= teb.fromASubDF[i];
      int jl= teb.fromASubFE[i];
      data1[ci++]=il;
      data1[cj++]=j*teb.nb_sub_fem+jl;      
    }
  //  warning the numbering of  
  for(int j=0;j<kk;++j)
    for(int i=0;i<No;++i)
      data1[ci++]=0;//j*m+teb.begin_dfcomp[i];
  for(int j=0;j<kk;++j)
    for(int i=0;i<No;++i)      
      data1[ci++]=m*kk;//j*m+teb.end_dfcomp[i];
  cout << " kk "<< kk << " " << m << " : ";
  //  for(int i=0;i< N*2;++i)
  //  cout << data1[2*n+i] << " " ;
  // cout << endl;    
}

FESumConstruct::FESumConstruct(int kk,const TypeOfFE **t)
 :k(kk),teb(t),NN(new int[kk+1]),DF(new int[kk+1]) , comp(new int[kk])
{  
   map<const TypeOfFE *,int> m;
   int i=k,j;    
   while(i--) // on va a l'envert pour avoir comp[i] <=i 
      m[teb[i]]=i;
    // l'ordre comp est important comp est croissant  mais pas de pb. 
   i=k;    
   while(i--) 
     comp[i]=m[teb[i]]; //  comp[i] <=i
        
  // reservatition des intervalles en espaces
  int n=0,N=0;
   for ( j=0;j<kk;j++)
     {NN[j]=N;N+=teb[j]->N;}
   NN[kk] = N;
 //  reservation des interval en df   
   n=0;
   for ( j=0;j<kk;j++)
    { DF[j]=n;n+=teb[j]->NbDoF;}
   DF[kk] = n;
//  n = nb de DF total   
//  N the fem is in R^N 
   
  data = new int [n*(5+2) + 3*N];