tetgen.cpp

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

    }
  assert(itet==in.numberofpoints*3);
	
  in.numberoffacets = 0;

  cout << "tetgen: before tetrahedralize( , &in, &out): switch=" << switch_tetgen<< endl;
  tetrahedralize(switch_tetgen, &in, &out);
	 
  cout << "tetgen: finish tetrahedralize( , &in, &out);" << endl;
  //mesh3_tetgenio_out( out, label_tet, label_face,*T_Th3);
  Mesh3* T_Th3=mesh3_tetgenio_out( out, label_tet, label_face);
  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; 
}


Mesh3 * RemplissageSurf3D_tetgen(char *switch_tetgen,const Mesh3 & Th3, const int & label_tet){
	
  //Mesh3 *T_Th3= new Mesh3;
	
  assert(Th3.nt == 0 );
  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);
	
  if(verbosity) 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; 
		
  }  
  cout << "tetgen: before tetrahedralize( , &in, &out);" << endl;
	
  tetrahedralize(switch_tetgen, &in, &out);
	 
  cout << "tetgen: after tetrahedralize( , &in, &out);" << endl;
  //mesh3_tetgenio_out( out, label_tet, *T_Th3);
  Mesh3 *T_Th3 = mesh3_tetgenio_out( out, label_tet);
  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;
}

Mesh3 * RemplissageSurf3D_tetgen_new(char *switch_tetgen,const Mesh3 & Th3, const int & label_tet,
				     const int &nbhole, const double *tabhole, 
				     const int & nbregion, const double *tabregion, 
				     const int &nbfacecl, const double *tabfacecl){
	
  //Mesh3 *T_Th3= new Mesh3;
	
  assert(Th3.nt == 0 );
  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);
	
  if(verbosity) 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];
  }

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

Mesh3 * Transfo_Mesh2_tetgen(const double &precis_mesh, char *switch_tetgen,const Mesh & Th2, const double *tab_XX, const double *tab_YY, const double *tab_ZZ, 
	int &border_only, int &recollement_border, int &point_confondus_ok, 
	const int &label_tet, const map<int, int> &maptri ){
  
  //Mesh3 *T_Th3= new Mesh3;
        int nv_t,nt_t,nbe_t;
        int* Numero_Som;
	
	int* ind_nv_t;
	int* ind_nt_t;
	int* ind_nbe_t;
	
	int* label_nbe_t;
	
	//int i_som;
	Numero_Som = new int[Th2.nv];
	ind_nv_t   = new int[Th2.nv];
	ind_nbe_t  = new int[Th2.nt];
	     
	label_nbe_t = new int[Th2.nt];
	
        if(verbosity) cout << "2D: Mesh::Vertex  triangle2  border " << Th2.nv << " "<<Th2.nt<< " " << Th2.neb<< endl;
	
	for(int ii=0; ii<Th2.nv; ii++){ 
		Numero_Som[ii]=ii;
	}
	if(verbosity) cout <<" debut: SamePointElement " <<endl;
	
	SamePointElement_Mesh2( precis_mesh, tab_XX, tab_YY, tab_ZZ, Th2, recollement_border, point_confondus_ok, 
		Numero_Som, ind_nv_t, ind_nt_t, ind_nbe_t, label_nbe_t, nv_t, nt_t, nbe_t);
	
	if(verbosity) cout <<" fin: SamePointElement " <<endl;
	
	if(verbosity) cout << "2D transfo: Mesh::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++)
	{
		int & ii = ind_nv_t[nnv];
		//cout << "nnv ,  ii  =" << nnv << "  " << ii << endl;
		//cout << "tab_XX[ii], tab_YY[ii], tab_ZZ[ii]=" <<  tab_XX[ii] << " "<< tab_YY[ii] << " "<< tab_ZZ[ii] << endl;
		assert( Numero_Som[ii] == nnv );
		const Mesh::Vertex & K = Th2.vertices[ii];//const Mesh::Vertex & K(Th2.vertices[ii]); //Version Mesh2   
		in.pointlist[itet]   = tab_XX[ii];
		in.pointlist[itet+1] = tab_YY[ii];
		in.pointlist[itet+2] = tab_ZZ[ii];       
		in.pointmarkerlist[nnv] =  K.lab;   
		itet=itet+3;
	}
	assert(itet==in.numberofpoints*3);
	
	if(verbosity) 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];
		
	  int & ii=ind_nbe_t[ibe];
	  // creation of elements
	  const Mesh::Triangle & K(Th2.t(ii)); // const Triangle2 & K(Th2.elements[ii]); // Version Mesh2  
	  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 
		
	}  

	cout << "tetgen: before tetrahedralize( , &in, &out);" << endl;
	
	tetrahedralize(switch_tetgen, &in, &out);
	 
	cout << "tetgen: after tetrahedralize( , &in, &out);" << endl;
	//mesh3_tetgenio_out( out, label_tet, *T_Th3);
	Mesh3 *T_Th3=mesh3_tetgenio_out( out, label_tet);
	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;


}

Mesh3 * Transfo_Mesh2_tetgen_new(const double &precis_mesh, char *switch_tetgen,const Mesh & Th2, const double *tab_XX, const double *tab_YY, const double *tab_ZZ, 
				 int &border_only, int &recollement_border, int &point_confondus_ok, 
				 const int &label_tet, const map<int, int> &maptri, 
				 const int &nbhole, const double *tabhole, const int & nbregion, const double *tabregion, 
				 const int &nbfacecl, const double *tabfacecl){
  //Mesh3 *T_Th3= new Mesh3;
  int nv_t,nt_t,nbe_t;
  int* Numero_Som;
	
  int* ind_nv_t;
  int* ind_nt_t;
  int* ind_nbe_t;
	
  int* label_nbe_t;
	
  //int i_som;
  Numero_Som = new int[Th2.nv];
  ind_nv_t   = new int[Th2.nv];
  ind_nbe_t  = new int[Th2.nt];
	     
  label_nbe_t = new int[Th2.nt];
	
  if(verbosity) cout << "2D: Mesh::Vertex  triangle2  border " << Th2.nv << " "<<Th2.nt<< " " << Th2.neb<< endl;
	
  for(int ii=0; ii<Th2.nv; ii++){ 
    Numero_Som[ii]=ii;
  }
  if(verbosity) cout <<" debut: SamePointElement " <<endl;
	
  SamePointElement_Mesh2( precis_mesh, tab_XX, tab_YY, tab_ZZ, Th2, recollement_border, point_confondus_ok, 
			  Numero_Som, ind_nv_t, ind_nt_t, ind_nbe_t, label_nbe_t, nv_t, nt_t, nbe_t);
	
  if(verbosity) cout <<" fin: SamePointElement " <<endl;
	
  if(verbosity) cout << "2D transfo: Mesh::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++)
    {
      int & ii = ind_nv_t[nnv];
      //cout << "nnv ,  ii  =" << nnv << "  " << ii << endl;
      //cout << "tab_XX[ii], tab_YY[ii], tab_ZZ[ii]=" <<  tab_XX[ii] << " "<< tab_YY[ii] << " "<< tab_ZZ[ii] << endl;
      assert( Numero_Som[ii] == nnv );
      const Mesh::Vertex & K = Th2.vertices[ii];//const Mesh::Vertex & K(Th2.vertices[ii]); //Version Mesh2   
      in.pointlist[itet]   = tab_XX[ii];
      in.pointlist[itet+1] = tab_YY[ii];
      in.pointlist[itet+2] = tab_ZZ[ii];       
      in.pointmarkerlist[nnv] =  K.lab;   
      itet=itet+3;
    }
  assert(itet==in.numberofpoints*3);
	
  if(verbosity) 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];
		
    int & ii=ind_nbe_t[ibe];
    // creation of elements
    const Mesh::Triangle & K(Th2.t(ii)); // const Triangle2 & K(Th2.elements[ii]); // Version Mesh2