⭐ 欢迎来到虫虫下载站! | 📦 资源下载 📁 资源专辑 ℹ️ 关于我们
⭐ 虫虫下载站

📄 tetgen.cpp

📁 FreeFem++可以生成高质量的有限元网格。可以用于流体力学
💻 CPP
📖 第 1 页 / 共 5 页
字号:
      }     AnyType operator()(Stack stack)  const ;};basicAC_F0::name_and_type ReconstructionRefine_Op::name_param[]= {  {  "switch", &typeid(string*)},  {  "reftet", &typeid(KN_<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>)},  {  "sizeofvolume", &typeid(double)}};class  ReconstructionRefine : public OneOperator { public:    ReconstructionRefine() : OneOperator(atype<pmesh3>(),atype<pmesh3>()) {}    E_F0 * code(const basicAC_F0 & args) const   {     return  new ReconstructionRefine_Op( args, t[0]->CastTo(args[0]) );   }};AnyType ReconstructionRefine_Op::operator()(Stack stack)  const {  MeshPoint *mp(MeshPointStack(stack)) , mps=*mp;  Mesh3 * pTh= GetAny<Mesh3 *>((*eTh)(stack));  //double msvol= GetAny<double>((*maxvol)(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 << "refine tetgen: Vertex Triangle Border " << nbv<< "  "<< nbt << " "<< nbe << endl;    KN<long> zzempty;  //int intempty=0;  string stringempty= string("rqaAAYQC");  string* switch_tet(arg(0,stack,&stringempty));  KN<long> nrtet(arg(1,stack,zzempty));    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];      }    }  ffassert( nrtet.N() %2 ==0);  map<int,int> maptet;  for(int i=0;i<nrtet.N();i+=2)    {      if(nrtet[i] != nrtet[i+1]){	maptet[nrtet[i]]=nrtet[i+1];      }    }  KN<double> tsizevol(nbt);  MeshPoint *mp3(MeshPointStack(stack)); 	  R3 Cdg_hat = R3(1./4.,1./4.,1./4.);  for (int it=0; it<nbt; it++){    Tet & K(Th.elements[it]);    mp3->set( Th, K(Cdg_hat), Cdg_hat, K, K.lab);    if( nargs[9]){      tsizevol[it]=  GetAny< double >( (*nargs[9])(stack) );    }    else if(tabregion.N()==0){      for(int i=0; i< nbregion;i++){	if( K.lab == tabregion[3+5*i] ){	  tsizevol[it] = tabregion[4+5*i];	}      }          }    else{      tsizevol[it] = K.mesure();    }  }  cout << "Before reconstruction:" << " nbhole="   << nbhole << " nbregion=" << nbregion << endl;   Mesh3 *Th3 = ReconstructionRefine_tetgen(switch_tetgen, Th, nbhole, tabhole, nbregion, tabregion, nbfacecl,tabfacecl,tsizevol);	  cout << "finish reconstruction " << endl;  // changement de label 1   if( nrtet.N() > 0){    for(int ii=0; ii< Th3->nt; ii++){      const Tet & K(Th3->elements[ii]);      int lab;      int iv[4];			      iv[0] = Th3->operator()(K[0]);      iv[1] = Th3->operator()(K[1]);      iv[2] = Th3->operator()(K[2]);      iv[3] = Th3->operator()(K[3]);      map< int, int> :: const_iterator imap;      imap = maptet.find(K.lab);      if( imap != maptet.end() ){	lab = imap -> second;      }      else{	lab = K.lab;      }      Th3->elements[ii].set(Th3->vertices,iv,lab);    }  }  if( nrf.N() > 0){    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 << "fin du changement de label " << endl;    Th3->BuildBound();  Th3->BuildAdj();  Th3->Buildbnormalv();    Th3->BuildjElementConteningVertex();  Th3->BuildGTree();  //Th3->decrement();      Add2StackOfPtr2FreeRC(stack,Th3);  delete [] switch_tetgen;  *mp=mps;  cout << "FreeFem++: End check mesh given by tetgen" << endl;    return Th3;}// ConvexHull3D_tetg_Opclass ConvexHull3D_tetg_Op : public E_F0mps {public:  Expression numofpts;  Expression xx,yy,zz;  static const int n_name_param =3; //   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:  ConvexHull3D_tetg_Op(const basicAC_F0 &  args, Expression nop,   Expression ffxx, Expression ffyy, Expression ffzz )     : numofpts(nop), xx(ffxx), yy(ffyy), zz(ffzz)  {    if(verbosity) cout << "Convex Hull with TetGen" << endl;    args.SetNameParam(n_name_param,name_param,nargs);  }     AnyType operator()(Stack stack)  const ;};basicAC_F0::name_and_type  ConvexHull3D_tetg_Op::name_param[]= {  {  "switch", &typeid(string*)},  {  "reftet", &typeid(long)},   {  "refface", &typeid(long)}};class ConvexHull3D_tetg : public OneOperator { public:       ConvexHull3D_tetg() : OneOperator( atype<pmesh3>(), atype<long>(),      atype< KN_<double> >(), atype< KN_<double> >(), atype< KN_<double> >() ) {}    E_F0 * code(const basicAC_F0 & args) const   { 	return  new  ConvexHull3D_tetg_Op( args,t[0]->CastTo(args[0]), 	      t[1]->CastTo(args[1]), t[2]->CastTo(args[2]), t[3]->CastTo(args[3]) );   }};AnyType  ConvexHull3D_tetg_Op::operator()(Stack stack)  const {  int nbv = (int) GetAny<long>((*numofpts)(stack));  KN<double> cxx(nbv),cyy(nbv),czz(nbv);    cxx = GetAny< KN<double> > ((*xx)(stack));  cyy = GetAny< KN<double> > ((*yy)(stack));  czz = GetAny< KN<double> > ((*zz)(stack));    assert( cxx.N() == nbv );  assert( cyy.N() == nbv );  assert( czz.N() == nbv );       KN<long> zzempty;  //int intempty=0;  string stringempty= string("YqaAAQC");  string* switch_tet(arg(0,stack,&stringempty));  int label_tet(arg(1,stack,0));    int label_face(arg(2,stack,1));  //====================================  //  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);   //======================================   Mesh3 *Th3 =  Convexhull_3Dpoints ( switch_tetgen, nbv, cxx, cyy, czz, label_tet, label_face );	  Th3->BuildBound();  Th3->BuildAdj();  Th3->Buildbnormalv();    Th3->BuildjElementConteningVertex();  Th3->BuildGTree();  //Th3->decrement();     Add2StackOfPtr2FreeRC(stack,Th3);    delete [] switch_tetgen;  cout << "FreeFem++: End check mesh given by tetgen" << endl;    return Th3;}// ConvexHull3D_tetg_file_Opclass ConvexHull3D_tetg_file_Op: public E_F0mps {public:  Expression filename;  static const int n_name_param =3; //   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:  ConvexHull3D_tetg_file_Op(const basicAC_F0 &  args, Expression zfilename )     : filename(zfilename)  {    if(verbosity) cout << "Convex Hull with TetGen" << endl;    args.SetNameParam(n_name_param,name_param,nargs);  }     AnyType operator()(Stack stack)  const ;};basicAC_F0::name_and_type  ConvexHull3D_tetg_file_Op::name_param[]= {  {  "switch", &typeid(string*)},  {  "reftet", &typeid(long)},   {  "refface", &typeid(long)}};class ConvexHull3D_tetg_file : public OneOperator { public:       ConvexHull3D_tetg_file() : OneOperator( atype<pmesh3>(), atype<string *>() ) {}    E_F0 * code(const basicAC_F0 & args) const   { 	return  new  ConvexHull3D_tetg_file_Op( args,t[0]->CastTo(args[0]) );   }};AnyType  ConvexHull3D_tetg_file_Op::operator()(Stack stack)  const {  string*  pointsfile = GetAny<string*>((*filename)(stack));  // lecture du fichier contenant les points  int nbv;  //int lec;   ifstream fp( pointsfile->c_str() );  if(!fp) {    cerr << "  -- tetgconvexhull : Erreur openning " << pointsfile <<endl;     exit(1);  }  if(verbosity>1)    cout << "  -- tetgconvexhull:  Read On file \"" << pointsfile <<"\""<<  endl;  fp >> nbv;   cout << "  -- Nb of Points " << nbv << endl;  KN<double> cxx(nbv), cyy(nbv), czz(nbv);    for(int lec=0;lec<nbv;lec++)    {      fp >> cxx[lec] >> cyy[lec] >> czz[lec];    }//  if( lec !=nbv ) {//     cerr << "  -- tetgconvexhull : Erreur Reading File " << pointsfile <<endl; //     cerr << " number of points " << nbv << " " <<"number of lecture" << lec << endl;  //     exit(1);//   }//   assert(lec==nbv);  fp.close();      KN<long> zzempty;  //int intempty=0;  string stringempty= string("YQV");  string* switch_tet(arg(0,stack,&stringempty));  int label_tet(arg(1,stack,0));    int label_face(arg(2,stack,1));  //====================================  //  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);   //======================================   Mesh3 *Th3 =new Mesh3;     Th3 = Convexhull_3Dpoints( switch_tetgen, nbv, cxx, cyy, czz, label_tet, label_face );	  Th3->BuildBound();  Th3->BuildAdj();  Th3->Buildbnormalv();    Th3->BuildjElementConteningVertex();  Th3->BuildGTree();  //Th3->decrement();      Add2StackOfPtr2FreeRC(stack,Th3);  delete [] switch_tetgen;  cout << "FreeFem++: End check mesh given by tetgen" << endl;    return Th3;}class Init1 { public:  Init1();};static Init1 init1;  //  une variable globale qui serat construite  au chargement dynamique Init1::Init1(){  // le constructeur qui ajoute la fonction "splitmesh3"  a freefem++     //if (verbosity)  if(verbosity) cout << " load: tetgen  " << endl;    Global.Add("tetgconvexhull","(",new ConvexHull3D_tetg_file);  Global.Add("tetgconvexhull","(",new ConvexHull3D_tetg);  Global.Add("tetgtransfo","(",new Build2D3D);  Global.Add("tetg","(",new Remplissage);  Global.Add("tetgreconstruction","(",new ReconstructionRefine);}// because i include this file in tetgen.cpp (very bad) FH// a will correct this in next version ...#define  WITH_NO_INIT#include "msh3.cpp" 

⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -