fespacen.cpp

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

    data0[9]=discon;
    
    return data0;
}

  
dataTypeOfFE::dataTypeOfFE(const int nitemdim[4],const KN< dataTypeOfFE const *>  &  tef)
: 
data(builddata_d(nitemdim,tef)),  
dataalloc(data),
ndfonVertex(data[0]),
ndfonEdge(data[1]),
ndfonFace(data[2]),
ndfonVolume(data[3]),
NbDoF(data[4]),
NbNode(data[5]),
N(data[6]),
nb_sub_fem(data[7]),
nbsubdivision(data[8]),
discontinue(data[9]),
DFOnWhat(data+10+0*NbDoF),
DFOfNode(data+10+1*NbDoF),
NodeOfDF(data+10+2*NbDoF),
fromFE(data+10+3*NbDoF),
fromDF(data+10+4*NbDoF),
fromASubFE(data+10+5*NbDoF),
fromASubDF(data+10+6*NbDoF) ,
dim_which_sub_fem(data+10+7*NbDoF)
{}

template<class Mesh>
void GTypeOfFESum<Mesh>::init(InterpolationMatrix<RdHat> & M,FElement * pK,int odf,int ocomp,int *pp) const
{
  // a faire ..... cas matrix invariante 
  assert(0);
}

template<class Mesh>
     GTypeOfFESum<Mesh>::GTypeOfFESum(const KN< GTypeOfFE<Mesh> const *> & t)
     : 
     GTypeOfFE<Mesh>(t),
     k(t.N()),
     teb(t),
     NN(k+1),
     DF(k+1) ,
     comp(k+1) {Build();}
     
template<class Mesh> 
static  KN< GTypeOfFE<Mesh> const *> kn(const GFESpace<Mesh> ** tt,int kk)
     {
	 KN< GTypeOfFE<Mesh> const *> r(kk);
	 for(int i=0;i<kk;++i)
	   { r[i]=tt[i]->TFE[0];ffassert(tt[i]->TFE.constant());}
	 return r;
     }
template<class Mesh> 
static     KN< GTypeOfFE<Mesh> const *> kn(const GFESpace<Mesh> & tt,int kk)
     {
	 return  KN< GTypeOfFE<Mesh> const *> (kk,tt.TFE[0]);
     }
     
template<class Mesh>
     GTypeOfFESum<Mesh>::GTypeOfFESum(const GFESpace<Mesh> ** tt,int kk)
     :	
     GTypeOfFE<Mesh>(kn(tt,kk)),
     k(kk),
     teb(kn(tt,kk)),
     NN(k+1),
     DF(k+1) ,
     comp(k+1) {Build();}
     
template<class Mesh>
     GTypeOfFESum<Mesh>::GTypeOfFESum(const GFESpace<Mesh> & tt,int kk)
     :	
     GTypeOfFE<Mesh>(kn(tt,kk)),
     k(kk),
     teb(kn(tt,kk)),
     NN(k+1),
     DF(k+1) ,
     comp(k+1) {Build();}
     
template<class Mesh>
void GTypeOfFESum<Mesh>::Build()
{
  bool debug=true;
  {
    const KN< GTypeOfFE<Mesh> const *> & t=teb;
    map<const GTypeOfFE<Mesh> *,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<k;j++)
      {NN[j]=N;N+=teb[j]->N;}
    NN[k] = N;
    //  reservation des interval en df   
    n=0;
    for ( j=0;j<k;j++)
      { DF[j]=n;n+=teb[j]->NbDoF;}
    DF[k] = n;
  }
  int ii=0;
  for (int i=0;i<k;++i)
    {
      for (int j=0;j<teb[i]->nb_sub_fem;++j)
	this->Sub_ToFE[ii++]=teb[i]->Sub_ToFE[j];
    }
  assert(ii==this->nb_sub_fem );
  
  int c=0,c0=0, fcom=0;
  for (int i=0;i<this->nb_sub_fem;i++) 
    { 
      int N=this->Sub_ToFE[i]->N;
      int ndofi=this->Sub_ToFE[i]->NbDoF;
      this->first_comp[i]= fcom;
      this->last_comp[i]= fcom+N;
      fcom += N;

      for(int j=0;j<N;++j)
	{
	  this->begin_dfcomp[c] = c0 + this->Sub_ToFE[i]->begin_dfcomp[j] ; 
	  this->end_dfcomp[c]   = c0 + this->Sub_ToFE[i]->end_dfcomp[j] ;
	  c++;
	}
      c0+=ndofi;
      
    }
  if(debug)
    {
      cout <<" NbDoF : " << this->NbDoF <<endl;
      for(int i=0;i<this->N;++i)
	cout << "      comp " << i << " ["<<this->begin_dfcomp[i]<<", "<< this->end_dfcomp[i]<< "[\n";
    }
  
  // construction de l'interpolation .
  
  int npi=0;
  int nci=0;
  bool var=true;
  for (int i=0;i<this->nb_sub_fem;i++)
    {
      npi +=this->Sub_ToFE[i]->NbPtforInterpolation;
      nci +=this->Sub_ToFE[i]->NbcoefforInterpolation;
      var = var && this->Sub_ToFE[i]->invariantinterpolationMatrix;
    }
  assert(this->NbcoefforInterpolation== nci);
  this->invariantinterpolationMatrix=var;
  // this->pInterpolation.init(nci);
  // this->cInterpolation.init(nci);
  // this->dofInterpolation.iniy(nci);
  {
    map<RdHat,int,lessRd> mpt;
    numPtInterpolation.init(npi);
    int npp=0,kkk=0;
    KN<RdHat> Ptt(npi);
    for (int i=0;i<this->nb_sub_fem;i++)
      {
	const GTypeOfFE<Mesh> &ti=*this->Sub_ToFE[i];
	
	for(int p=0;p<ti.NbPtforInterpolation;++p,++kkk)
	  {
	    Ptt[kkk]=ti.PtInterpolation[p];
	    if(verbosity>5)
	    cout << "    p= "<< p << " [ " << Ptt[kkk]<< "] ,  "<< kkk<< " "<< npp<<endl;;
	    if( mpt.find(Ptt[kkk]) == mpt.end())
	      mpt[Ptt[kkk]]=npp++;
	    numPtInterpolation[kkk]=mpt[Ptt[kkk]];
	  }
      }
    assert(this->NbPtforInterpolation==0);
    if(verbosity>5)
    cout << npp;
    this->NbPtforInterpolation=npp;
    this->PtInterpolation.init(npp);
    for(int i=0;i<npp;++i)
      this->PtInterpolation[numPtInterpolation[i]]=Ptt[i];
  }
  
  int oc=0,odof=0;
  for (int i=0,k=0;i<this->nb_sub_fem;i++)
    {
      const GTypeOfFE<Mesh> &ti=*this->Sub_ToFE[i];
      for(int j=0;j<ti.NbcoefforInterpolation; ++j,++k)
	{
	  this->pInterpolation[k]   = numPtInterpolation[ti.pInterpolation[j]];
	  this->cInterpolation[k]   = ti.cInterpolation[j]+oc;
	  this->dofInterpolation[k] = ti.dofInterpolation[j]+odof;
	  this->coefInterpolation[k]=ti.coefInterpolation[j];
	}
      oc += ti.N;
      odof += ti.NbDoF; 
    }
  
  assert(c==this->N);
}

template<class MMesh> 
GFESpace<MMesh>::GFESpace(const GFESpace & Vh,int kk)
     :
     GFESpacePtrTFE<MMesh>(new GTypeOfFESum<MMesh>(Vh,kk)),
     DataFENodeDF(Vh.Th.BuildDFNumbering(this->ptrTFE->ndfonVertex,this->ptrTFE->ndfonEdge,this->ptrTFE->ndfonFace,this->ptrTFE->ndfonVolume)),
     Th(Vh.Th),
     TFE(1,0,this->ptrTFE), 
     cmesh(Th),
     N(TFE[0]->N),
     Nproduit(TFE[0]->N),
     nb_sub_fem(TFE[0]->nb_sub_fem),
     dim_which_sub_fem(TFE[0]->dim_which_sub_fem),
     maxNbPtforInterpolation(TFE[0]->NbPtforInterpolation),
     maxNbcoefforInterpolation(TFE[0]->NbcoefforInterpolation)
     
     {
     }
    
template<class MMesh> 
     GFESpace<MMesh>::GFESpace(const GFESpace ** pVh,int kk)
     :
     GFESpacePtrTFE<MMesh>(new GTypeOfFESum<MMesh>(pVh,kk)),
     DataFENodeDF((**pVh).Th.BuildDFNumbering(this->ptrTFE->ndfonVertex,this->ptrTFE->ndfonEdge,this->ptrTFE->ndfonFace,this->ptrTFE->ndfonVolume)),
     Th((**pVh).Th),
     TFE(1,0,this->ptrTFE), 
     cmesh(Th),
     N(TFE[0]->N),
     Nproduit(TFE[0]->N),
     nb_sub_fem(TFE[0]->nb_sub_fem),
     dim_which_sub_fem(TFE[0]->dim_which_sub_fem),
     maxNbPtforInterpolation(TFE[0]->NbPtforInterpolation),
     maxNbcoefforInterpolation(TFE[0]->NbcoefforInterpolation)
     
     {
	 for(int i=0;i<kk;++i)
	     ffassert(&Th==&pVh[i]->Th);
     }
   
template<class MMesh>     
KN<R>  GFESpace<MMesh>::newSaveDraw(const KN_<R> & U,int componante,int & lg,KN<Rd> &Psub,KN<int> &Ksub,int op_U) const  
{
  const int d =  Rd::d;
  Rd *Ps;
  int *Ks;
  int nsb = TFE[0]->nbsubdivision;        
  int nvsub,nksub;
  SplitSimplex<Rd>(nsb, nvsub,  Ps,  nksub ,  Ks);
  ffassert( Psub.unset());
  ffassert( Ksub.unset());
  Psub.set(Ps,nvsub);
  Ksub.set(Ks,nksub*(d+1));
   lg= nvsub*Th.nt;
  KN<double> v(lg);
  for (int k=0,i=0;k<Th.nt;k++)
    { 
      FElement K=(*this)[k];
      for(int l=0;l<nvsub;l++)
	  v[i++] =   K(Psub[l], U, componante, op_U)  ;

    }                                                                                                                                                                                                                            
  return KN<double>(true,v);// to remove the copy.    
}

   /*
template<class MMesh>
KN<double>  GFESpace<MMesh>::newSaveDraw(const KN_<R> & U,int composante,int & lg,int & nsb) const 
     {
	 nsb = TFE[0]->nbsubdivision;
	 int nsbv = NbOfSubInternalVertices(nsb,d);
	 lg = nsbv*Th.nt;
	 cout << "newSaveDraw what: nt " << Th.nt << " " << nsbv << " " << lg << endl;
	 KN<double> v(lg);
	 ffassert(v);
	 for (int k=0,i=0;k<Th.nt;k++)
	   {
	       (*this)[k].SaveDraw( U,composante,&v[i]);	
	       i+=nsbv;
	   }
	 return KN<double>(true,v);// to remove the copy.
     }
   */     
// explicite instance..
template class GTypeOfFESum<Mesh2>;
template class GTypeOfFESum<Mesh3>;
template class GFESpace<Mesh1>;
template class GFESpace<Mesh2>;
template class GFESpace<Mesh3>;
    
}