tetgen.cpp

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

    iv[0] = out.trifacelist[3*ibe]-1;
    iv[1] = out.trifacelist[3*ibe+1]-1;
    iv[2] = out.trifacelist[3*ibe+2]-1; 
    for(int jj=0; jj<3; jj++){
      if(iv[jj]>= Th3.nv || iv[jj]< 0 ) cout << "iv[jj]=" << iv[jj] << " triangle" << ibe << endl;
      assert( iv[jj] >=0 && iv[jj] <Th3.nv );
    }
    Th3.be(ibe).set( Th3.vertices, iv, out.trifacemarkerlist[ibe]);
  }

  /*
  if( out.numberoftetrahedronattributes != 1 ){
    cout << "out.numberoftetrahedronattributes" << out.numberoftetrahedronattributes  << endl;
    exit(1);
   }
  */
}


void mesh3_tetgenio_out(const tetgenio &out, const int & label_tet, Mesh3 & Th3)
{ 
  int i;

// All indices start from 1.
  if(out.firstnumber != 1){
    cout << " probleme ???" << endl;
    exit(1);
  }   
  
  if(out.numberoffacets !=0){
    cout << "tetgen: faces non triangulaire" << endl;
    exit(1);
  }
  
  if(out.numberofcorners !=4){
    cout << "tetgen: element subparametric of order 2" <<endl;
    exit(1);
  }
  
  cout << "Th3 :: Vertex Element Border :: " << out.numberofpoints << " " <<out.numberoftetrahedra  << " " << out.numberoftrifaces << endl;
  Th3.set(out.numberofpoints, out.numberoftetrahedra, out.numberoftrifaces);
   
  i=0;
  for(int nnv=0; nnv < Th3.nv; nnv++){
    Th3.vertices[nnv].x=out.pointlist[i];
    Th3.vertices[nnv].y=out.pointlist[i+1];
    Th3.vertices[nnv].z=out.pointlist[i+2];       
    Th3.vertices[nnv].lab=out.pointmarkerlist[nnv];
    i=i+3;    
  }
    
  i=0;
  for(int nnt=0; nnt < Th3.nt; nnt++){
    int iv[4],lab;
    iv[0] = out.tetrahedronlist[i]-1;
    iv[1] = out.tetrahedronlist[i+1]-1;
    iv[2] = out.tetrahedronlist[i+2]-1;
    iv[3] = out.tetrahedronlist[i+3]-1;
    lab   = label_tet;
    //lab = out.tetrahedronmarkerlist[nnt];
    Th3.elements[nnt].set( Th3.vertices, iv, lab);
    i=i+4;
  }
 
  for(int ibe=0; ibe < Th3.nbe; ibe++){
    int iv[3];
    iv[0] = out.trifacelist[3*ibe]-1;
    iv[1] = out.trifacelist[3*ibe+1]-1;
    iv[2] = out.trifacelist[3*ibe+2]-1; 
    Th3.be(ibe).set( Th3.vertices, iv, out.trifacemarkerlist[ibe]);
  }

}

void mesh3_tetgenio_out(const tetgenio &out, const int & label_tet, const int & label_face, Mesh3 & Th3)
{ 
  int i;

// All indices start from 1.
  if(out.firstnumber != 1){
    cout << " probleme ???" << endl;
    exit(1);
  }   
  
  if(out.numberoffacets !=0){
    cout << "tetgen: faces non triangulaire" << endl;
    exit(1);
  }
  
  if(out.numberofcorners !=4){
    cout << "tetgen: element subparametric of order 2" <<endl;
    exit(1);
  }
  
  cout << "Th3 :: Vertex Element Border :: " << out.numberofpoints << " " <<out.numberoftetrahedra  << " " << out.numberoftrifaces << endl;
  Th3.set(out.numberofpoints, out.numberoftetrahedra, out.numberoftrifaces);
   
  i=0;
  for(int nnv=0; nnv < Th3.nv; nnv++){
    Th3.vertices[nnv].x=out.pointlist[i];
    Th3.vertices[nnv].y=out.pointlist[i+1];
    Th3.vertices[nnv].z=out.pointlist[i+2];       
    Th3.vertices[nnv].lab=out.pointmarkerlist[nnv];
    i=i+3;    
  }
    
  i=0;
  for(int nnt=0; nnt < Th3.nt; nnt++){
    int iv[4],lab;
    iv[0] = out.tetrahedronlist[i]-1;
    iv[1] = out.tetrahedronlist[i+1]-1;
    iv[2] = out.tetrahedronlist[i+2]-1;
    iv[3] = out.tetrahedronlist[i+3]-1;
    lab   = label_tet;
    //lab = out.tetrahedronmarkerlist[nnt];
    Th3.elements[nnt].set( Th3.vertices, iv, lab);
    i=i+4;
  }

  if(verbosity) cout << &out.trifacemarkerlist << endl;
  
  for(int ibe=0; ibe < Th3.nbe; ibe++){
    int iv[3];
    iv[0] = out.trifacelist[3*ibe]-1;
    iv[1] = out.trifacelist[3*ibe+1]-1;
    iv[2] = out.trifacelist[3*ibe+2]-1; 
    Th3.be(ibe).set( Th3.vertices, iv, label_face);
  }

}
// verison Mesh3 *

Mesh3 * mesh3_tetgenio_out(const tetgenio &out)
{ 
  int i;

// All indices start from 1.
  if(out.firstnumber != 1){
    cout << " probleme ???" << endl;
    exit(1);
  }   
  
  if(out.numberoffacets !=0){
    cout << "tetgen: faces non triangulaire" << endl;
    exit(1);
  }
  
  if(out.numberofcorners !=4){
    cout << "tetgen: element subparametric of order 2" <<endl;
    exit(1);
  }

  cout << "Th3 :: Vertex Element Border :: " << out.numberofpoints << " " <<out.numberoftetrahedra  << " " << out.numberoftrifaces << endl;
  //Th3.set(out.numberofpoints, out.numberoftetrahedra, out.numberoftrifaces);

  // new parameter
  if( out.numberoftetrahedronattributes !=  1){
    cout << "out.numberoftetrahedronattributes" << out.numberoftetrahedronattributes << endl;
  }
   
  Vertex3 *v = new Vertex3[out.numberofpoints];
  Tet *t  = new Tet[out.numberoftetrahedra];
  Tet *tt = t;
  Triangle3 *b  = new Triangle3[out.numberoftrifaces];
  Triangle3 *bb = b;

  i=0;
  for(int nnv=0; nnv<out.numberofpoints; nnv++){
    v[nnv].x=out.pointlist[i];
    v[nnv].y=out.pointlist[i+1];
    v[nnv].z=out.pointlist[i+2];       
    v[nnv].lab=out.pointmarkerlist[nnv];
    i=i+3;    
  }

  // test pour la distance minimale entre les points 
//   {
//     double dist,dist1;
//     dist = 1000000000000.;
//     for(int nnv=0; nnv<out.numberofpoints; nnv++){
//       for(int nnv1=nnv+1; nnv1< out.numberofpoints; nnv1++){
// 	dist1=(v[nnv].x-v[nnv1].x)*(v[nnv].x-v[nnv1].x)+(v[nnv].y-v[nnv1].y)*(v[nnv].y-v[nnv1].y)
// 	  +(v[nnv].z-v[nnv1].z)*(v[nnv].z-v[nnv1].z);
// 	dist=min(dist,sqrt(dist1));
// 	if( sqrt(dist1) < 1e-12){ 
// 	  cout << "point confondus" << nnv << "<--->" <<nnv1 << endl;
// 	  if( sqrt( pow(v[nnv].x,2)+pow(v[nnv].y,2)+pow(v[nnv].z,2) ) > 1e-10   )  cout << v[nnv] << " " << v[nnv1] << endl;
// 	}
//       }
//     }
//     cout << "dist entre les points du maillage tetgen" << dist << endl;
//   }
    
  i=0;
  for(int nnt=0; nnt < out.numberoftetrahedra; nnt++){
    int iv[4],lab;
    iv[0] = out.tetrahedronlist[i]-1;
    iv[1] = out.tetrahedronlist[i+1]-1;
    iv[2] = out.tetrahedronlist[i+2]-1;
    iv[3] = out.tetrahedronlist[i+3]-1;
    //lab   = label_tet;
    for(int jj=0; jj<4; jj++){
      assert( iv[jj] >=0 && iv[jj] < out.numberofpoints );
    }
    
    //cout << "nnt= " <<  nnt << " " << lab  << " " << out.tetrahedronattributelist[nnt] << endl;
    lab =  out.tetrahedronattributelist[nnt];
    //cout << "nnt= " <<  lab  << " " << out.tetrahedronattributelist[nnt] << endl;
    //Th3.elements[nnt].set( Th3.vertices, iv, lab);
    (*tt++).set( v, iv, lab);   
    i=i+4;
  }
 
  for(int ibe=0; ibe < out.numberoftrifaces; ibe++){
    int iv[3];
    iv[0] = out.trifacelist[3*ibe]-1;
    iv[1] = out.trifacelist[3*ibe+1]-1;
    iv[2] = out.trifacelist[3*ibe+2]-1; 
    for(int jj=0; jj<3; jj++){
      if(iv[jj]>= out.numberofpoints || iv[jj]< 0 ) cout << "iv[jj]=" << iv[jj] << " triangle" << ibe << endl;
      assert( iv[jj] >=0 && iv[jj] < out.numberofpoints );
    }
    //Th3.be(ibe).set( Th3.vertices, iv, out.trifacemarkerlist[ibe]);
    (*bb++).set( v, iv, out.trifacemarkerlist[ibe]);   
  }

  Mesh3 *T_TH3 = new Mesh3(out.numberofpoints, out.numberoftetrahedra, out.numberoftrifaces, v, t, b);
  cout << "FreeFem++: Check mesh given by tetgen" << endl;
  //return T_TH3;

  if( TestElementMesh3(*T_TH3) != 1){  
    return T_TH3;
  }
  else{
    //Mesh3 *T2_TH3 = TestElementMesh3_patch( *T_TH3 );
    //return T2_TH3;
    exit(1);
  }

}


Mesh3* mesh3_tetgenio_out(const tetgenio &out, const int & label_tet)
{ 
  int i;

// All indices start from 1.
  if(out.firstnumber != 1){
    cout << " probleme ???" << endl;
    exit(1);
  }   
  
  if(out.numberoffacets !=0){
    cout << "tetgen: faces non triangulaire" << endl;
    exit(1);
  }
  
  if(out.numberofcorners !=4){
    cout << "tetgen: element subparametric of order 2" <<endl;
    exit(1);
  }
  
  cout << "Th3 :: Vertex Element Border :: " << out.numberofpoints << " " <<out.numberoftetrahedra  << " " << out.numberoftrifaces << endl;
  //Th3.set(out.numberofpoints, out.numberoftetrahedra, out.numberoftrifaces);
   
  Vertex3 *v = new Vertex3[out.numberofpoints];
  Tet *t  = new Tet[out.numberoftetrahedra];
  Tet *tt = t;
  Triangle3 *b  = new Triangle3[out.numberoftrifaces];
  Triangle3 *bb = b;
    
   
  i=0;
  for(int nnv=0; nnv < out.numberofpoints ; nnv++){
    v[nnv].x=out.pointlist[i];
    v[nnv].y=out.pointlist[i+1];
    v[nnv].z=out.pointlist[i+2];       
    v[nnv].lab=out.pointmarkerlist[nnv];
    i=i+3;    
  }
    
  i=0;
  for(int nnt=0; nnt < out.numberoftetrahedra; nnt++){
    int iv[4],lab;
    iv[0] = out.tetrahedronlist[i]-1;
    iv[1] = out.tetrahedronlist[i+1]-1;
    iv[2] = out.tetrahedronlist[i+2]-1;
    iv[3] = out.tetrahedronlist[i+3]-1;
    lab   = label_tet;
    //lab = out.tetrahedronmarkerlist[nnt];
    //Th3.elements[nnt].set( Th3.vertices, iv, lab);
    (*tt++).set( v, iv, lab);   
    i=i+4;
  }
 
  for(int ibe=0; ibe < out.numberoftrifaces ; ibe++){
    int iv[3];
    iv[0] = out.trifacelist[3*ibe]-1;
    iv[1] = out.trifacelist[3*ibe+1]-1;
    iv[2] = out.trifacelist[3*ibe+2]-1; 
    //Th3.be(ibe).set( Th3.vertices, iv, out.trifacemarkerlist[ibe]);
    (*bb++).set( v, iv, out.trifacemarkerlist[ibe]);  
  }
  Mesh3 *T_TH3 = new Mesh3(out.numberofpoints, out.numberoftetrahedra, out.numberoftrifaces, v, t, b);
  //return T_TH3;
  cout << "FreeFem++: Check mesh given by tetgen" << endl;
  if( TestElementMesh3(*T_TH3) != 1){  
    return T_TH3;
  }
  else{ 
    //cout << "patch pour tetgen " << endl;
    //Mesh3 *T2_TH3 = TestElementMesh3_patch( *T_TH3 );
    //return T2_TH3;
    exit(1);
  }
}

Mesh3* mesh3_tetgenio_out(const tetgenio &out, const int & label_tet, const int & label_face)
{ 
  int i;

// All indices start from 1.
  if(out.firstnumber != 1){
    cout << " probleme ???" << endl;
    exit(1);
  }   
  
  if(out.numberoffacets !=0){
    cout << "tetgen: faces non triangulaire" << endl;
    exit(1);
  }
  
  if(out.numberofcorners !=4){
    cout << "tetgen: element subparametric of order 2" <<endl;
    exit(1);
  }
  
  cout << "Th3 :: Vertex Element Border :: " << out.numberofpoints << " " <<out.numberoftetrahedra  << " " << out.numberoftrifaces << endl;
  //Th3.set(out.numberofpoints, out.numberoftetrahedra, out.numberoftrifaces);

  Vertex3 *v = new Vertex3[out.numberofpoints];
  Tet *t  = new Tet[out.numberoftetrahedra];
  Tet *tt = t;
  Triangle3 *b  = new Triangle3[out.numberoftrifaces];
  Triangle3 *bb = b;

  i=0;
  for(int nnv=0; nnv < out.numberofpoints; nnv++){
    v[nnv].x=out.pointlist[i];
    v[nnv].y=out.pointlist[i+1];
    v[nnv].z=out.pointlist[i+2];       
    v[nnv].lab=out.pointmarkerlist[nnv];
    i=i+3;    
  }
    
  i=0;
  for(int nnt=0; nnt < out.numberoftetrahedra; nnt++){
    int iv[4],lab;
    iv[0] = out.tetrahedronlist[i]-1;
    iv[1] = out.tetrahedronlist[i+1]-1;
    iv[2] = out.tetrahedronlist[i+2]-1;
    iv[3] = out.tetrahedronlist[i+3]-1;
    lab   = label_tet;
    //lab = out.tetrahedronmarkerlist[nnt];
    //Th3.elements[nnt].set( Th3.vertices, iv, lab);
    (*tt++).set( v, iv, lab); 
    i=i+4;
  }

  if(verbosity) cout << &out.trifacemarkerlist << endl;
  
  for(int ibe=0; ibe < out.numberoftrifaces ; ibe++){
    int iv[3];
    iv[0] = out.trifacelist[3*ibe]-1;
    iv[1] = out.trifacelist[3*ibe+1]-1;
    iv[2] = out.trifacelist[3*ibe+2]-1; 
    //Th3.be(ibe).set( Th3.vertices, iv, label_face);
    (*bb++).set( v, iv, label_face);  
  }

  Mesh3 *T_TH3 = new Mesh3(out.numberofpoints, out.numberoftetrahedra, out.numberoftrifaces, v, t, b);

  if( TestElementMesh3( *T_TH3 ) != 1){  
    return T_TH3;
  }
  else{
    exit(1);
    //Mesh3 *T2_TH3 = TestElementMesh3_patch( *T_TH3 );
    //return T2_TH3;
  }
}



Mesh3 * Convexhull_3Dpoints(char *switch_tetgen, const int &nv_t, const double *Xcoord, const double *Ycoord, 
const double *Zcoord, const int &label_tet, const int &label_face){
	
  //Mesh3 *T_Th3= new Mesh3;
		
  tetgenio in,out;
  //tetgenio::facet *f;
  //tetgenio::polygon *p;

  cout << " tetgenio: vertex " << endl;
  int itet,jtet;

  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]   = Xcoord[nnv];
      in.pointlist[itet+1] = Ycoord[nnv];
      in.pointlist[itet+2] = Zcoord[nnv];       
      in.pointmarkerlist[nnv] =  0;
      
      itet=itet+3;