📄 tetgen.cpp
字号:
} 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 + -