tetgen.cpp

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

    p->vertexlist[0] = Numero_Som[ Th2.operator()(K[0]) ]+1;
    p->vertexlist[1] = Numero_Som[ Th2.operator()(K[1]) ]+1;
    p->vertexlist[2] = Numero_Som[ Th2.operator()(K[2]) ]+1;
		
    for( int kkk=0; kkk<3; kkk++){ 
      assert( p->vertexlist[kkk]<= in.numberofpoints && p->vertexlist[kkk]> 0);
    }
    map< int, int>:: const_iterator imap;
    imap = maptri.find(K.lab); // imap= maptri.find( label_nbe_t[ibe] );
    assert( imap != maptri.end());
    in.facetmarkerlist[ibe] = imap->second; // K.lab; // before 
		
  }  

  // mise a jour des nouvelles variables

  in.numberofholes = nbhole;
  in.holelist = new REAL[3*nbhole];
  
  for(int ii=0; ii<3*in.numberofholes; ii++){
    in.holelist[ii]   = tabhole[ii];
  }


  in.numberofregions = nbregion;
  in.regionlist = new REAL[5*nbregion];
  for(int ii=0; ii<5*in.numberofregions; ii++){
    in.regionlist[ii] = tabregion[ii];
  }

  in.numberoffacetconstraints = nbfacecl;
  in.facetconstraintlist = new REAL[2*in.numberoffacetconstraints];

  for(int ii=0; ii<2*in.numberoffacetconstraints; ii++){
    in.facetconstraintlist[ii+1] = tabfacecl[ii+1];
  }
  
  cout << "tetgen: before tetrahedralize( , &in, &out);" << endl;
  
  tetrahedralize(switch_tetgen, &in, &out);
  
  cout << "tetgen: after tetrahedralize( , &in, &out);" << endl;
  //mesh3_tetgenio_out( out, *T_Th3);
  Mesh3 *T_Th3=mesh3_tetgenio_out( out);
  cout <<" Finish Mesh3 :: Vertex, Element, Border" << T_Th3->nv << " "<< T_Th3->nt << " " << T_Th3->nbe << endl;

  delete [] Numero_Som;
  delete [] ind_nv_t;
  delete [] ind_nbe_t;
  delete [] label_nbe_t;
  
  cout << "FreeFem++: End check mesh given by tetgen" << endl;  
  return T_Th3;
}

// Fonction Refine avec tetgen

Mesh3 * ReconstructionRefine_tetgen(char *switch_tetgen,const Mesh3 & Th3, 
				     const int &nbhole, const double *tabhole, 
				     const int & nbregion, const double *tabregion, 
				     const int &nbfacecl, const double *tabfacecl, const double *tsizevol){

  // verif option refine
  int i;
  //Mesh3 *T_Th3= new Mesh3;	
  assert(Th3.nt != 0 );
  {
    size_t testr, testp;
    int lenswitch;
    const char* test_tetgen = switch_tetgen;
    
    testr = strcspn(test_tetgen,"r");
    testp = strcspn(test_tetgen,"p");
    //cout << "switch_tetgen=" << switch_tetgen << endl;
    //cout << "testr=" << testr  <<" testp="<< testp <<  endl;
    //cout << "strlen(test_tetgen)" << strlen(test_tetgen) << endl; 
    if( testr == strlen(test_tetgen) )
      { 
	cout << "The option 'r' of tetgen is not used" << endl;
	exit(1);
      }
    testp = strcspn(test_tetgen,"p");
    if( testp != strlen(test_tetgen) )
      { 
	cout << "With TetGen :: the option 'p' is not possible to use with option 'r' " << endl;
	exit(1);
      }
  }

  int nv_t = Th3.nv;
  int nt_t = Th3.nt;
  int nbe_t = Th3.nbe;
	
  if(verbosity) cout << "3D RemplissageSurf3D:: Vertex  triangle2  border "  << nv_t << " "<< nt_t << " " << nbe_t<< endl;
  // Creation des tableau de tetgen
	
  tetgenio in,out;
  //tetgenio::facet *f;
  //tetgenio::polygon *p;
	
  if(verbosity) cout << " tetgenio: vertex " << endl;
  int itet,jtet;
  // All indices start from 1.
  in.firstnumber = 1;
  in.numberofpoints = nv_t;
  in.pointlist = new REAL[in.numberofpoints*3];
  in.pointmarkerlist = new int[in.numberofpoints];
  itet=0;
  jtet=0;
  for(int nnv=0; nnv < nv_t; nnv++)
    {
      in.pointlist[itet]   = Th3.vertices[nnv].x;
      in.pointlist[itet+1] = Th3.vertices[nnv].y;
      in.pointlist[itet+2] = Th3.vertices[nnv].z;       
      in.pointmarkerlist[nnv] =  Th3.vertices[nnv].lab;   
      itet=itet+3;
    }
  assert(itet==in.numberofpoints*3);

  // Tetrahedrons
  if(verbosity) cout << "tetrahedrons" << endl;
  in.numberofcorners = 4;
  in.numberoftetrahedra = Th3.nt;
  in.tetrahedronlist = new int[in.numberofcorners*in.numberoftetrahedra];
  in.numberoftetrahedronattributes = 1;
  in.tetrahedronattributelist = new REAL[in.numberoftetrahedronattributes*in.numberoftetrahedra];
  
  in.tetrahedronvolumelist = new REAL[in.numberoftetrahedra];

  i=0;
  for(int nnt=0; nnt < Th3.nt; nnt++){
    const Tet & K(Th3.elements[nnt]);

    in.tetrahedronlist[i]   = Th3.operator()(K[0])+1;
    in.tetrahedronlist[i+1] = Th3.operator()(K[1])+1;
    in.tetrahedronlist[i+2] = Th3.operator()(K[2])+1;
    in.tetrahedronlist[i+3] = Th3.operator()(K[3])+1;

    /*
    if( Th3.vertices[Th3.operator()(K[0])].y > 3.14) {
    in.tetrahedronvolumelist[nnt] = 0.001;
    }
    else{
    in.tetrahedronvolumelist[nnt] = 0.01;
    }
    */
    
    in.tetrahedronvolumelist[nnt] = tsizevol[nnt];

    in.tetrahedronattributelist[nnt] = K.lab;
    
    i=i+4;
  }

  
  if(verbosity) cout << "lecture des facettes" << endl;
  in.numberoftrifaces = Th3.nbe;
  in.trifacelist = new int[3*in.numberoftrifaces];
  in.trifacemarkerlist = new int[in.numberoftrifaces];

  for(int ibe=0; ibe < Th3.nbe; ibe++){
    const Triangle3 & K(Th3.be(ibe));
    
    in.trifacelist[3*ibe]   = Th3.operator()(K[0])+1;
    in.trifacelist[3*ibe+1] = Th3.operator()(K[1])+1;
    in.trifacelist[3*ibe+2] = Th3.operator()(K[2])+1; 					      
    in.trifacemarkerlist[ibe] = K.lab;

  }
  
  
  /*
  cout << " tetgenio: facet " << endl;
  // Version avec des facettes

  in.numberoffacets = nbe_t;
  in.facetlist = new tetgenio::facet[in.numberoffacets];
  in.facetmarkerlist = new int[in.numberoffacets];
	  
  for(int ibe=0; ibe < nbe_t; ibe++){
    tetgenio::facet *f;
    tetgenio::polygon *p;
    f = &in.facetlist[ibe];
    f->numberofpolygons = 1;
    f->polygonlist = new tetgenio::polygon[f->numberofpolygons];
    f->numberofholes = 0;
    f->holelist = NULL;
		
    p = &f->polygonlist[0];
    p->numberofvertices = 3;
    p->vertexlist = new int[3];
		
    // creation of elements
    const Triangle3 & K(Th3.be(ibe)); // const Triangle2 & K(Th2.elements[ii]); // Version Mesh2  
    p->vertexlist[0] = Th3.operator()(K[0])+1;
    p->vertexlist[1] = Th3.operator()(K[1])+1;
    p->vertexlist[2] = Th3.operator()(K[2])+1;
		
    for( int kkk=0; kkk<3; kkk++){ 
      assert( p->vertexlist[kkk]<= in.numberofpoints && p->vertexlist[kkk]>0);
    }
		
    in.facetmarkerlist[ibe] = K.lab; 
		
  }  
  */

  // mise a jour des nouvelles variables
 
  in.numberofholes = nbhole;
  in.holelist = new REAL[3*nbhole];
  
  for(int ii=0; ii<3*in.numberofholes; ii++){
    in.holelist[ii]   = tabhole[ii];
    if(verbosity) cout << "in.holelist[ii]=" << in.holelist[ii] << endl;
  }

  in.numberofregions = nbregion;
  in.regionlist = new REAL[5*nbregion];
  for(int ii=0; ii<5*in.numberofregions; ii++){
    in.regionlist[ii] = tabregion[ii];
    if(verbosity) cout << "in.regionlist[ii]=" << in.regionlist[ii] << endl;
  }

  in.numberoffacetconstraints = nbfacecl;
  in.facetconstraintlist = new REAL[2*in.numberoffacetconstraints];

  for(int ii=0; ii<2*in.numberoffacetconstraints; ii++){
    in.facetconstraintlist[ii+1] = tabfacecl[ii+1];
  }

  cout << "tetgen: before tetrahedralize( , &in, &out);" << endl;
  cout << "numberof regions "<<  in.numberofregions  << endl;
  cout << "numberof hole "<<  in.numberofholes  << endl;
	
  tetrahedralize(switch_tetgen, &in, &out);
	 
  cout << "tetgen: after tetrahedralize( , &in, &out);" << endl;
  //mesh3_tetgenio_out( out, *T_Th3);
  Mesh3 *T_Th3=mesh3_tetgenio_out( out);
  cout <<" Finish Mesh3 tetgen :: Vertex, Element, Border" << T_Th3->nv << " "<< T_Th3->nt << " " << T_Th3->nbe << endl;
  cout << "FreeFem++: End check mesh given by tetgen" << endl;  
  return T_Th3;
}


// declaration pour FreeFem++

class Remplissage_Op : public E_F0mps 
{
public:
  Expression eTh;
  static const int n_name_param =9; // 
  static basicAC_F0::name_and_type name_param[] ;
  Expression nargs[n_name_param];
  KN_<long>  arg(int i,Stack stack,KN_<long> a ) const
  { return nargs[i] ? GetAny<KN_<long> >( (*nargs[i])(stack) ): a;}
  KN_<double>  arg(int i,Stack stack,KN_<double> a ) const
  { return nargs[i] ? GetAny<KN_<double> >( (*nargs[i])(stack) ): a;}
  double  arg(int i,Stack stack,double a ) const{ return nargs[i] ? GetAny< double >( (*nargs[i])(stack) ): a;}
  int  arg(int i,Stack stack, int a ) const{ return nargs[i] ? GetAny< int >( (*nargs[i])(stack) ): a;}
  string*  arg(int i,Stack stack, string* a ) const{ return nargs[i] ? GetAny< string* >( (*nargs[i])(stack) ): a;}
public:
  Remplissage_Op(const basicAC_F0 &  args,Expression tth) 
    : eTh(tth)
  {
    if(verbosity >1) cout << "Remplissage du bord" << endl;
    args.SetNameParam(n_name_param,name_param,nargs);
  } 
  
  AnyType operator()(Stack stack)  const ;
};


basicAC_F0::name_and_type Remplissage_Op::name_param[]= {
  {  "switch", &typeid(string*)},
  {  "reftet", &typeid(long)}, 
  {  "refface", &typeid(KN_<long> )},
  // new parmameters
  {  "nbofholes", &typeid(long)},
  {  "holelist", &typeid(KN_<double>)},
  {  "nbofregions", &typeid(long)},
  {  "regionlist", &typeid(KN_<double>)},
  {  "nboffacetcl", &typeid(long)},
  {  "facetcl", &typeid(KN_<double>)}
};

class  Remplissage : public OneOperator { public:  
    Remplissage() : OneOperator(atype<pmesh3>(),atype<pmesh3>()) {}
  
  E_F0 * code(const basicAC_F0 & args) const 
  { 
	return  new Remplissage_Op( args,t[0]->CastTo(args[0]) ); 
  }
};

AnyType Remplissage_Op::operator()(Stack stack)  const 
{
  MeshPoint *mp(MeshPointStack(stack)) , mps=*mp;
  Mesh3 * pTh= GetAny<Mesh3 *>((*eTh)(stack));
  ffassert( pTh );
  Mesh3 &Th=*pTh;
  Mesh3 *m= pTh;   // question a quoi sert *m ??
  int nbv=Th.nv; // nombre de sommet 
  int nbt=Th.nt; // nombre de triangles
  int nbe=Th.nbe; // nombre d'aretes fontiere
  cout << "Tetgen : Vertex Triangle Border " << nbv<< "  "<< nbt << " nbe "<< nbe << endl; 
 
  KN<long> zzempty;
  //int intempty=0;
  string stringempty= string("pqaAAYQC");
  string* switch_tet(arg(0,stack,&stringempty));
  int label_tet(arg(1,stack,0));  
  KN<long> nrf (arg(2,stack,zzempty));

  // new parameters
  KN<double> zdzempty;
  int nbhole (arg(3,stack,0));
  KN<double> tabhole (arg(4,stack,zdzempty));
  int nbregion (arg(5,stack,0));
  KN<double> tabregion (arg(6,stack,zdzempty));
  int nbfacecl (arg(7,stack,0));
  KN<double> tabfacecl (arg(8,stack,zdzempty));

  
  // assertion au niveau de la taille 
  ffassert( tabhole.N()   == 3*nbhole);
  ffassert( tabregion.N() == 5*nbregion);
  ffassert( tabfacecl.N() == 2*nbfacecl);

  

  //====================================
  //  How to change string* into char* 
  //====================================
  size_t size_switch_tet = switch_tet->size()+1;
  char* switch_tetgen =new char[size_switch_tet];
  strncpy(switch_tetgen, switch_tet->c_str(), size_switch_tet); 
   
  ffassert( nrf.N() %2 ==0);

  map<int,int> mapf;
  for(int i=0;i<nrf.N();i+=2)
    {
      if(nrf[i] != nrf[i+1]){
	mapf[nrf[i]]=nrf[i+1];
      }
    }
  //Mesh3 *Th3 = Th3 = RemplissageSurf3D_tetgen( switch_tetgen, Th, label_tet);
  
  cout << "tetgen:" << "nbhole="   << nbhole << "nbregion=" << nbregion << endl;

  cout << "check :: orientation des surfaces" << endl;
  Th.BuildSurfaceAdj();
  cout << "fin check :: orientation des surfaces" << endl;
 
  Mesh3 *Th3 = RemplissageSurf3D_tetgen_new( switch_tetgen, Th, label_tet, nbhole, tabhole, nbregion, tabregion, nbfacecl,tabfacecl);
	
  cout << "finish tetgen " << endl;
  // changement de label 
  if( nrf.N() > 0){
    cout << "changement de label" << endl;
    for(int ii=0; ii< Th3->nbe; ii++){
      const Triangle3 & K(Th3->be(ii));
      int lab;
      int iv[3];
			
      iv[0] = Th3->operator()(K[0]);
      iv[1] = Th3->operator()(K[1]);
      iv[2] = Th3->operator()(K[2]);
			
      map< int, int> :: const_iterator imap;
      imap = mapf.find(K.lab);
      if( imap != mapf.end() ){
	lab = imap -> second;
      }
      else{
	lab = K.lab;
      }
      Th3->be(ii).set(Th3->vertices,iv,lab);
    }
  }
  cout << "action sur le maillage" << endl;

  Th3->BuildBound();
  Th3->BuildAdj();
  Th3->Buildbnormalv();  
  Th3->BuildjElementConteningVertex();
  Th3->BuildGTree();
  //Th3->decrement();    
  Add2StackOfPtr2FreeRC(stack,Th3);
  
  *mp=mps;
  delete [] switch_tetgen; 
  cout << "FreeFem++: End check mesh given by tetgen" << endl;  
  return Th3;
}

// Refine et Recontruction 

class ReconstructionRefine_Op : public E_F0mps 
{
public:
  Expression eTh;
  static const int n_name_param =10; // 
  static basicAC_F0::name_and_type name_param[] ;
  Expression nargs[n_name_param];
  KN_<long>  arg(int i,Stack stack,KN_<long> a ) const
  { return nargs[i] ? GetAny<KN_<long> >( (*nargs[i])(stack) ): a;}
  KN_<double>  arg(int i,Stack stack,KN_<double> a ) const
  { return nargs[i] ? GetAny<KN_<double> >( (*nargs[i])(stack) ): a;}
  double  arg(int i,Stack stack,double a ) const{ return nargs[i] ? GetAny< double >( (*nargs[i])(stack) ): a;}
  int  arg(int i,Stack stack, int a ) const{ return nargs[i] ? GetAny< int >( (*nargs[i])(stack) ): a;}
  string*  arg(int i,Stack stack, string* a ) const{ return nargs[i] ? GetAny< string* >( (*nargs[i])(stack) ): a;}
public:
  ReconstructionRefine_Op(const basicAC_F0 &  args,Expression tth) 
    : eTh(tth)
  {
    if(verbosity >1) cout << "ReconstructionRefine du bord"<< endl;
    args.SetNameParam(n_name_param,name_param,nargs);