|
mFES - molecular Finite Element Solver
0.4
|
00001 00012 //#include <solve.hpp> 00013 #include <iostream> 00014 #include <fstream> 00015 #include <string> 00016 #include <sstream> 00017 #include <boost/timer.hpp> 00018 00019 #include "LSMS.h" 00020 #include "Voxel.h" 00021 #include "VCG.h" 00022 00023 00024 using namespace std; 00025 using namespace vcg; 00026 00027 00028 namespace nglib { 00029 #include <nglib.h> 00030 } 00031 using namespace nglib; 00032 00033 00034 typedef boost::property_tree::ptree INI; 00035 00056 class Model { 00057 public: 00060 static bool exclusionSurface; 00061 00078 void calcIonLayer(vector<Atom> &atomList, INI &ini, string fname = "exclusion.vol"){ 00079 LSMS lmSurface; 00080 mMesh mSurface; 00081 int mask = 1; 00082 00085 if (!boost::filesystem::exists( "exclusion.stl" )){ 00086 cout << "Calculating ion exclusion layer " << fname << "..." << endl; 00087 boost::timer t; 00088 ofstream time; 00091 time.open ("times"); 00093 if (lmSurface.calcMC(mSurface, atomList, ini, "exclusion")){ 00095 clean(mSurface); 00098 smooth(mSurface, ini); 00101 tri::io::ExporterSTL<mMesh>::Save(mSurface,"exclusion.stl",false); 00102 00104 double currentTime = t.elapsed(); 00105 time << "exclusion_surface " << currentTime << " s" << endl; 00106 cout << "mFES: ion exclusion layer " << fname << "surface calculation took " << currentTime << " seconds." <<endl; 00107 time.close(); 00108 } 00109 } 00110 } 00111 00112 00128 void calcModel(vector<Atom> &atomList, INI &ini, string fname = "", bool cavity = false){ 00129 00131 LSMS lmSurface; 00133 mMesh mSurface; 00134 00136 string volume_vol = ini.get<string>("model.volume_vol"); 00138 string surface_stl = ini.get<string>("model.surface_stl"); 00142 string generator = ini.get<string>("model.generator"); 00144 int generatorResolution = atoi(ini.get<string>("model.grid_resolution").c_str()); 00145 00147 string generatorResidue = ini.get_optional<string>("model.generator_residue").get_value_or("standard"); 00150 int generatorModelResolution = atoi(ini.get_optional<string>("model.grid_residue_resolution").get_value_or("256").c_str()); 00151 00153 boost::optional<string> ionc = ini.get_optional<string>("experiment.ionc"); 00154 00157 if(ionc){ 00158 string ionc = ini.get<string>("experiment.ionc"); 00159 if (atof(ionc.c_str()) != 0){ 00160 exclusionSurface = true; 00161 } 00162 } 00163 00164 00165 int mask = 1; 00166 if (generator == "standard" && cavity && !boost::filesystem::exists( "cavity.vol" )){ 00169 cout << "Calculating cavities..." << endl; 00170 00172 boost::timer t; 00173 ofstream time; 00174 time.open ("times"); 00175 00177 if (lmSurface.calcMC(mSurface, atomList, ini, "cavity")){ 00179 clean(mSurface); 00180 smooth(mSurface, ini); 00181 00183 if (surface_stl != "" ) 00184 tri::io::ExporterSTL<mMesh>::Save(mSurface,"cavity.stl",false); 00185 00186 double currentTime = t.elapsed(); 00187 time << "cavity_surface " << currentTime << " s" << endl; 00188 cout << "mFES: cavity model surface calculation took " << currentTime << " seconds." <<endl; 00189 00190 t.restart(); 00193 convert(mSurface, ini, "cavity.vol", "cavity"); 00194 00196 currentTime = t.elapsed(); 00197 time << "cavity_volume " << currentTime << " s" << endl; 00198 cout << "mFES: cavity model VOL meshing took " << currentTime << " seconds." <<endl; 00199 } 00200 time.close(); 00201 00202 } else if ( fname == "" && !boost::filesystem::exists( volume_vol ) ){ 00203 // Protein model is calculated 00204 boost::timer t; 00205 ofstream time; 00206 time.open ("times"); 00207 00209 if (boost::filesystem::exists( surface_stl )) { 00211 cout << surface_stl << " found. " << endl; 00212 tri::io::ImporterSTL<mMesh>::Open(mSurface, surface_stl.c_str(), mask); 00213 } else { 00215 if (generator == "standard"){ 00218 LSMS tSurface; 00219 mMesh tmSurface; 00222 tSurface.calcMC(tmSurface, atomList, ini); 00224 clean(tmSurface); 00225 smooth(tmSurface, ini); 00228 if (surface_stl == "" ) 00229 surface_stl = "protein.stl"; 00230 00232 tri::io::ExporterSTL<mMesh>::Save(tmSurface,surface_stl.c_str(),false); 00236 tri::io::ImporterSTL<mMesh>::Open(mSurface, surface_stl.c_str(), mask); 00237 } else { 00239 Voxel vSurface; 00242 if (!boost::filesystem::exists( "protein.stl" )){ 00245 vSurface.calcSurface(mSurface, atomList, ini, "protein.stl", generatorResolution); 00246 00248 clean(mSurface); 00249 00250 } 00253 tri::io::ImporterSTL<mMesh>::Open(mSurface, string("protein.stl").c_str(), mask); 00254 } 00255 00256 } 00257 00259 double currentTime = t.elapsed(); 00260 time << "protein_surface " << currentTime << " s" << endl; 00261 cout << "mFES: protein model surface calculation took " << currentTime << " seconds." <<endl; 00262 00263 t.restart(); 00267 convert(mSurface, ini); 00268 currentTime = t.elapsed(); 00269 time << "protein_volume " << currentTime << " s" << endl; 00270 cout << "mFES: protein model VOL meshing took " << currentTime << " seconds." <<endl; 00271 time.close(); 00272 00273 } else if (fname != ""){ 00276 string exclusionstlFile = fname+"_exclusion.stl"; 00277 string exclusionGroups = ini.get<string>("pka.explicit_models"); 00278 00282 if (!boost::filesystem::exists( exclusionstlFile.c_str() ) ){ 00283 cout << "group model calculation " << fname << endl; 00284 00285 if (generatorResidue == "standard"){ 00288 lmSurface.calcMC(mSurface, atomList, ini, "residue_exclusion"); 00289 clean(mSurface); 00290 smooth(mSurface, ini); 00293 tri::io::ExporterSTL<mMesh>::Save(mSurface,exclusionstlFile.c_str(),false); 00294 } else { 00297 Voxel vSurface; 00298 vSurface.calcSurface(mSurface, atomList, ini, exclusionstlFile, generatorModelResolution, true); 00299 clean(mSurface); 00301 tri::io::ImporterSTL<mMesh>::Open(mSurface, exclusionstlFile.c_str(), mask); 00302 } 00303 00304 } 00305 00308 if ( !boost::filesystem::exists( fname ) ){ 00309 string stlFile = fname+".stl"; 00310 ofstream time; 00311 00312 time.open ("times", ios::app); 00313 boost::timer t; 00314 // group model calculation 00315 if ( !boost::filesystem::exists( stlFile )){ 00316 00319 if (generatorResidue == "standard"){ 00320 lmSurface.calcMC(mSurface, atomList, ini, "residue"); 00321 clean(mSurface); 00322 smooth(mSurface, ini); 00325 tri::io::ExporterSTL<mMesh>::Save(mSurface,stlFile.c_str(),false); 00326 } else { 00329 Voxel vSurface; 00330 vSurface.calcSurface(mSurface, atomList, ini, stlFile, generatorModelResolution); 00331 clean(mSurface); 00332 tri::io::ImporterSTL<mMesh>::Open(mSurface, stlFile.c_str(), mask); 00333 } 00334 00335 double currentTime = t.elapsed(); 00336 time << "model_surface " << currentTime << " s" << endl; 00337 cout << "mFES: residue model surface calculation took " << currentTime << " seconds." <<endl; 00338 00339 t.restart(); 00343 convert(mSurface, ini, fname, "residue"); 00346 currentTime = t.elapsed(); 00347 time << "model_volume " << currentTime << " s" << endl; 00348 cout << "mFES: residue model VOL meshing took " << currentTime << " seconds." <<endl; 00349 time.close(); 00350 } 00351 else { 00355 cout << "model surface found: " << stlFile << endl; 00356 tri::io::ImporterSTL<mMesh>::Open(mSurface, stlFile.c_str(), mask); 00357 t.restart(); 00361 convert(mSurface, ini, fname, "residue"); 00362 double currentTime = t.elapsed(); 00363 time << "model_volume " << currentTime << " s" << endl; 00364 cout << "mFES: residue model VOL meshing took " << currentTime << " seconds." <<endl; 00365 time.close(); 00366 } 00367 } 00368 } 00369 } 00370 00371 00372 private: 00373 00386 void printMeshingOptions(Ng_Meshing_Parameters &mp, string prefix){ 00387 cout << endl; 00388 cout << prefix << endl; 00389 cout << "========================" << endl; 00390 cout << "options.localh = " << mp.uselocalh << endl; 00391 cout << "options.meshsize = " << mp.maxh << endl; 00392 cout << "options.minmeshsize = " << mp.minh << endl; 00393 cout << "meshoptions.fineness = " << mp.fineness << endl; 00394 cout << "meshoptions.blockfill = " << mp.blockfill << endl; 00395 cout << "meshoptions.filldist = " << mp.filldist << endl; 00396 cout << "options.grading = " << mp.grading << endl; 00397 cout << "options.curvaturesafety = " << mp.elementspercurve << endl; 00398 cout << "options.segmentsperedge = " << mp.elementsperedge << endl; 00399 cout << "options.secondorder = " << mp.second_order << endl; 00400 cout << "options.quad = " << mp.quad_dominated << endl; 00401 cout << "stloptions.resthcloseedgeenable = " << mp.closeedgeenable << endl; 00402 cout << "stloptions.resthcloseedgefac = " << mp.closeedgefact << endl; 00403 cout << "options.optsteps2d = " << mp.optsteps_2d << endl; 00404 cout << "options.optsteps3d = " << mp.optsteps_3d << endl; 00405 cout << "options.inverttets = " << mp.invert_tets << endl; 00406 cout << "options.inverttrigs = " << mp.invert_trigs << endl; 00407 cout << "options.checkoverlap = " << mp.check_overlap << endl << endl; 00408 } 00409 00424 void setMeshingOptions(Ng_Meshing_Parameters &mp, string optFile){ 00426 ifstream in(optFile.c_str()); 00427 if (!in) { 00428 in.close(); 00429 cout << "Cannot open mesh options file:" << optFile << endl; 00430 exit(0); 00431 } 00432 string pqrline; 00433 string variable; 00434 string value; 00435 00436 while( !in.eof() ) { 00437 getline(in, pqrline); 00438 00439 istringstream ss(pqrline); 00441 ss >> variable >> value; 00442 00445 if (variable == "meshoptions.blockfill"){ 00446 mp.blockfill = atoi(value.c_str()); 00447 } 00448 00449 if (variable == "meshoptions.filldist"){ 00450 mp.filldist = atof(value.c_str()); 00451 } 00452 00453 00454 if (variable == "options.localh"){ 00456 mp.uselocalh = atoi(value.c_str()); 00457 continue; 00458 } 00459 00460 if (variable == "options.meshsize"){ 00462 mp.maxh = atof(value.c_str()); 00463 continue; 00464 } 00465 00466 if (variable == "options.minmeshsize"){ 00468 mp.minh = atof(value.c_str()); 00469 continue; 00470 } 00471 00472 if (variable == "meshoptions.fineness"){ 00474 mp.fineness = atof(value.c_str()); 00475 continue; 00476 } 00477 00478 if (variable == "options.grading"){ 00480 mp.grading = atof(value.c_str()); 00481 continue; 00482 } 00483 00484 00485 if (variable == "options.curvaturesafety"){ 00487 mp.elementspercurve = atoi(value.c_str()); 00488 continue; 00489 } 00490 00491 if (variable == "options.segmentsperedge"){ 00493 mp.elementsperedge = atof(value.c_str()); 00494 continue; 00495 } 00496 00497 if (variable == "options.secondorder"){ 00499 mp.second_order = atoi(value.c_str()); 00500 continue; 00501 } 00502 00503 if (variable == "options.quad"){ 00505 mp.quad_dominated = atoi(value.c_str()); 00506 continue; 00507 } 00508 00509 // if (variable == "options.meshsizefilename"){ 00510 // //!< Optional external mesh size file 00511 // mp.meshsize_filename = value.c_str(); 00512 // continue; 00513 // } 00514 00515 if (variable == "stloptions.resthcloseedgeenable"){ 00517 mp.closeedgeenable = atoi(value.c_str()); 00518 continue; 00519 } 00520 00521 if (variable == "stloptions.resthcloseedgefac"){ 00523 mp.closeedgefact = atof(value.c_str()); 00524 continue; 00525 } 00526 00527 if (variable == "options.optsteps2d"){ 00529 mp.optsteps_2d = atoi(value.c_str()); 00530 continue; 00531 } 00532 00533 if (variable == "options.optsteps3d"){ 00535 mp.optsteps_3d = atoi(value.c_str()); 00536 continue; 00537 } 00538 00539 if (variable == "options.inverttets"){ 00541 mp.invert_tets = atoi(value.c_str()); 00542 continue; 00543 } 00544 00545 if (variable == "options.inverttrigs"){ 00547 mp.invert_trigs = atoi(value.c_str()); 00548 continue; 00549 } 00550 00551 if (variable == "options.checkoverlap"){ 00553 mp.check_overlap = atoi(value.c_str()); 00555 mp.check_overlapping_boundary = atoi(value.c_str()); 00556 } 00557 00558 00559 } 00560 // optimize by default is 1, because if number of 00561 // optsteps is 0, then no optimization made 00562 00563 in.close(); 00564 } 00565 00587 int convert(mMesh &mSurface, INI &ini, string fname = "", string mode = "protein") 00588 { 00589 00591 Ng_STL_Geometry *stl_geom = Ng_STL_NewGeometry(); 00592 00594 Ng_Init(); 00595 00597 ngVolume = Ng_NewMesh(); 00598 00599 int np, ne; 00600 00602 00603 double p1[3]; 00604 double p2[3]; 00605 double p3[3]; 00606 double n[3]; 00607 00608 mMesh::FaceIterator fi; 00609 00610 Point3f p; 00611 int triangles = 0; 00612 for(fi=mSurface.face.begin(); fi!=mSurface.face.end(); ++fi) if( !(*fi).IsD() ){ 00613 // For each triangle write the normal, the three coords and a short set to zero 00614 p.Import(vcg::NormalizedNormal(*fi)); 00615 n[0] = p[0]; n[1] = p[1]; n[2] = p[2]; 00616 00617 p.Import((*fi).V(0)->P()); 00618 p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2]; 00619 p.Import((*fi).V(1)->P()); 00620 p2[0] = p[0]; p2[1] = p[1]; p2[2] = p[2]; 00621 p.Import((*fi).V(2)->P()); 00622 p3[0] = p[0]; p3[1] = p[1]; p3[2] = p[2]; 00623 triangles++; 00624 Ng_STL_AddTriangle(stl_geom, p1, p2, p3, n); 00625 } 00626 00627 cout << triangles << " triangles in surface" << endl; 00628 if(!stl_geom) { 00629 cout << "Error reading VCG STL File" << endl; 00630 exit(1); 00631 } 00632 cout << "Successfully loaded VCG STL File" << endl; 00633 00635 string debug = ini.get<string>("model.debug"); 00636 string jobname = ini.get<string>("general.jobname"); 00637 string meshMoleculeSurface, meshMoleculeVolume; 00638 if (mode == "protein" || mode == "cavity" ){ 00639 meshMoleculeSurface = ini.get<string>("meshing.molecule_surface"); 00640 meshMoleculeVolume = ini.get<string>("meshing.molecule_volume"); 00641 } else if (mode == "residue") { 00642 meshMoleculeSurface = ini.get<string>("meshing.residue_surface"); 00643 meshMoleculeVolume = ini.get<string>("meshing.residue_volume"); 00644 } 00645 cout << "using: " << meshMoleculeSurface << endl; 00646 cout << "using: " << meshMoleculeVolume << endl; 00647 00648 Ng_Meshing_Parameters mp; 00649 setMeshingOptions(mp, meshMoleculeSurface); 00650 mp.optsurfmeshenable = 1; 00651 mp.optvolmeshenable = 1; 00652 00653 printMeshingOptions(mp, "Meshing options for molecular surface"); 00654 00655 cout << "Initialise the STL Geometry structure...." << endl; 00656 ngSurface = Ng_STL_InitSTLGeometry(stl_geom); 00657 if(ngSurface != NG_OK) { 00658 cout << "Error Initialising the STL Geometry....Aborting!!" << endl; 00659 exit(1); 00660 } 00661 00662 cout << "Start Edge Meshing...." << endl; 00663 ngSurface = Ng_STL_MakeEdges(stl_geom, ngVolume, &mp); 00664 if(ngSurface != NG_OK) { 00665 cout << "Error in Edge Meshing....Aborting!!" << endl; 00666 exit(1); 00667 } 00668 00669 cout << "Start Surface Meshing...." << endl; 00670 ngSurface = Ng_STL_GenerateSurfaceMesh(stl_geom, ngVolume, &mp); 00671 if(ngSurface != NG_OK) { 00672 cout << "Error in Surface Meshing....Aborting!!" << endl; 00673 exit(1); 00674 } 00675 00676 if (ini.get_optional<string>("meshing.second_order_surface").get_value_or("no") == "yes"){ 00677 cout << "Generating second order geometry mesh" << endl; 00678 Ng_STL_Generate_SecondOrder(stl_geom, ngVolume); 00679 } 00680 00681 if (debug == "analyze" && mode == "protein"){ 00682 stringstream fileName; 00683 fileName << jobname << "_surface_so.vol"; 00684 Ng_SaveMesh(ngVolume,fileName.str().c_str()); 00685 } 00686 00687 if (mode == "protein" && boost::filesystem::exists( "cavity.vol" )){ 00688 string cavity = "cavity.vol"; 00689 Ng_Mesh* cSurface; 00690 cSurface = nglib::Ng_LoadMesh(cavity.c_str()); 00691 00692 Ng_SetProperties(ngVolume, 1, 1, 1, 0); 00693 Ng_SetProperties(cSurface, 1, 1, 1, 0); 00694 00695 00696 cout << "Merging Mesh with cavity....." << endl; 00697 ngSurface = Ng_MergeMesh( ngVolume, cSurface ); 00698 if(ngSurface != NG_OK) { 00699 cout << "Error in cavity merging....Aborting!!" << endl; 00700 exit(1); 00701 } 00702 00703 } 00704 00705 // Example: Set up a global refinement if using a local h file 00706 // float globalH = 1e6; 00707 // cout << "setting global H to " << globalH << endl; 00708 // Ng_RestrictMeshSizeGlobal(ngVolume, globalH); 00709 00710 setMeshingOptions(mp, meshMoleculeVolume); 00711 mp.optsurfmeshenable = 1; 00712 mp.optvolmeshenable = 1; 00713 00714 printMeshingOptions(mp, "Meshing options for molecular volume"); 00715 00716 cout << "Start Volume Meshing of molecule...." << endl; 00717 ngSurface = Ng_GenerateVolumeMesh (ngVolume, &mp); 00718 if(ngSurface != NG_OK) { 00719 cout << "Error in Volume Meshing....Aborting!!" << endl; 00720 exit(1); 00721 } 00722 00723 // Example: Saving a volume file (VOL, NETGEN format) in current state 00724 // Ng_SaveMesh(ngVolume,"meshed.vol"); 00725 00726 00727 // Example: Generate a second order surface mesh 00728 // if (ini.get<string>("meshing.second_order_surface") == "yes") 00729 // Ng_STL_Generate_SecondOrder (stl_geom, ngVolume); 00730 00732 if (exclusionSurface){ 00733 00734 string exclusionVol = "exclusion.vol"; 00735 string exclusionStl = "exclusion.stl"; 00736 00737 if (fname != "") 00738 exclusionVol = fname+"_exclusion.vol"; 00739 if (fname != "") 00740 exclusionStl = fname+"_exclusion.stl"; 00741 00742 Ng_Mesh *exclusionVolume; 00743 exclusionVolume = Ng_NewMesh(); 00744 setMeshingOptions(mp, meshMoleculeSurface); 00745 mp.optsurfmeshenable = 1; 00746 mp.optvolmeshenable = 0; 00747 00748 if (!boost::filesystem::exists( exclusionVol )){ 00749 Ng_STL_Geometry *exclusion_geom = Ng_STL_NewGeometry(); 00750 00751 cout << "Loading exclusion " << exclusionStl << endl; 00752 exclusion_geom = Ng_STL_LoadGeometry( exclusionStl.c_str()); 00753 if(!exclusion_geom) 00754 { 00755 cout << "Error reading in STL File: " << exclusionStl << endl; 00756 return 1; 00757 } 00758 cout << "Successfully loaded STL File: " << exclusionStl << endl; 00759 00760 cout << "Initialise the STL Geometry structure for ion exclusion layer...." << endl; 00761 ngSurface = Ng_STL_InitSTLGeometry(exclusion_geom); 00762 if(ngSurface != NG_OK) { 00763 cout << "Error Initialising the STL Geometry for ion exclusion layer....Aborting!!" << endl; 00764 exit(1); 00765 } 00766 00767 cout << "Start Edge Meshing...." << endl; 00768 ngSurface = Ng_STL_MakeEdges(exclusion_geom, exclusionVolume, &mp); 00769 if(ngSurface != NG_OK) { 00770 cout << "Error in Edge Meshing ion exclusion layer....Aborting!!" << endl; 00771 exit(1); 00772 } 00773 00774 cout << "Start Surface Meshing...." << endl; 00775 ngSurface = Ng_STL_GenerateSurfaceMesh(exclusion_geom, exclusionVolume, &mp); 00776 if(ngSurface != NG_OK) { 00777 cout << "Error in Surface Meshing ion exclusion layer....Aborting!!" << endl; 00778 exit(1); 00779 } 00780 Ng_SaveMesh(exclusionVolume,exclusionVol.c_str()); 00781 00782 } else { 00783 exclusionVolume = nglib::Ng_LoadMesh(exclusionVol.c_str()); 00784 } 00785 00786 Ng_SetProperties(ngVolume, 1, 1, 1, 0); 00787 Ng_SetProperties(exclusionVolume, 1, 1, 1, 0); 00788 00789 cout << "Merging Mesh with ion exclusion layer....." << endl; 00790 ngSurface = Ng_MergeMesh( exclusionVolume, ngVolume ); 00791 if(ngSurface != NG_OK) { 00792 cout << "Error in exclusion layer merging....Aborting!!" << endl; 00793 exit(1); 00794 } 00795 00796 mp.optsurfmeshenable = 1; 00797 mp.optvolmeshenable = 1; 00798 00799 cout << "Start Volume meshing of exclusion layer...." << endl; 00800 ngSurface = Ng_GenerateVolumeMesh (exclusionVolume, &mp); 00801 if(ngSurface != NG_OK) { 00802 cout << "Error in Volume Meshing of exclusion layer....Aborting!!" << endl; 00803 exit(1); 00804 } 00805 cout << "Meshing successfully completed....!!" << endl; 00806 00807 Ng_Mesh* bSurface; 00808 00809 cout << "Loading boundary settings ....." << endl; 00810 string boundary = ini.get<string>("model.boundary"); 00811 bSurface = nglib::Ng_LoadMesh(boundary.c_str()); 00812 00813 Ng_SetProperties(bSurface, 1, 1, 1, 0); 00814 00815 cout << "Merging Mesh with boundary....." << endl; 00816 ngSurface = Ng_MergeMesh( bSurface, exclusionVolume ); 00817 if(ngSurface != NG_OK) { 00818 cout << "Error in Surface merging with ion exclusion....Aborting!!" << endl; 00819 exit(1); 00820 } 00821 00822 cout << "Merging complete" << endl; 00823 00824 string meshBoundaryVolume = ini.get<string>("meshing.boundary_volume"); 00825 setMeshingOptions(mp, meshBoundaryVolume); 00826 00827 mp.optsurfmeshenable = 1; 00828 mp.optvolmeshenable = 1; 00829 printMeshingOptions(mp, "Meshing options for boundary volume"); 00830 00831 00832 cout << "Start Volume meshing of whole model...." << endl; 00833 ngSurface = Ng_GenerateVolumeMesh (bSurface, &mp); 00834 if(ngSurface != NG_OK) { 00835 cout << "Error in Volume Meshing....Aborting!!" << endl; 00836 exit(1); 00837 } 00838 00839 // volume mesh output 00840 np = Ng_GetNP(bSurface); 00841 cout << "Points: " << np << endl; 00842 00843 ne = Ng_GetNE(bSurface); 00844 cout << "Elements: " << ne << endl; 00845 00846 string volumeVol = ini.get<string>("model.volume_vol"); 00847 00848 if (fname != "") 00849 volumeVol = fname; 00850 00851 if (volumeVol != ""){ 00852 cout << "Saving Mesh in VOL Format...." << endl; 00853 Ng_SaveMesh(bSurface,volumeVol.c_str()); 00854 } 00855 } else { 00856 00858 Ng_Mesh* bSurface; 00859 00860 if (mode != "cavity"){ 00861 cout << "Loading boundary settings ....." << endl; 00862 string boundary = ini.get<string>("model.boundary"); 00863 bSurface = nglib::Ng_LoadMesh(boundary.c_str()); 00864 00865 if (ini.get_optional<string>("meshing.second_order_surface").get_value_or("no") == "yes"){ 00866 cout << "Generating second order for volume mesh" << endl; 00867 Ng_Generate_SecondOrder(bSurface); 00868 } 00869 00870 Ng_SetProperties(ngVolume, 1, 1, 1, 0); 00871 Ng_SetProperties(bSurface, 1, 1, 1, 0); 00872 00873 cout << "Merging Mesh with boundary....." << endl; 00874 ngSurface = Ng_MergeMesh( bSurface, ngVolume ); 00875 if(ngSurface != NG_OK) { 00876 cout << "Error in Surface merging....Aborting!!" << endl; 00877 exit(1); 00878 } 00879 00880 cout << "Merging complete" << endl; 00881 00882 00883 string refineFile = ini.get_optional<string>("model.refine_file").get_value_or(""); 00884 00887 if (refineFile != ""){ 00888 // global refinement 00889 // void Ng_RestrictMeshSizeGlobal (Ng_Mesh * mesh, double h); 00890 cout << "Setting local refinement ....." << endl; 00891 ifstream in(refineFile.c_str()); 00892 if (!in) { 00893 in.close(); 00894 cout << "Cannot open refinement file:" << refineFile << endl; 00895 exit(0); 00896 } 00897 00898 string currentLine; 00899 double p[3]; double h; 00900 unsigned int refPoints = 0; 00901 while( !in.eof() ) { 00902 getline(in, currentLine); 00903 if (currentLine != ""){ 00904 istringstream ss(currentLine); 00905 ss >> p[0] >> p[1] >> p[2] >> h; 00906 Ng_RestrictMeshSizePoint (bSurface, p, h); 00907 refPoints++; 00908 } 00909 } 00910 cout << refPoints << " local refinement point(s) set." << endl; 00911 } 00912 00914 string meshBoundaryVolume = ini.get<string>("meshing.boundary_volume"); 00915 setMeshingOptions(mp, meshBoundaryVolume); 00916 mp.optsurfmeshenable = 1; 00917 mp.optvolmeshenable = 1; 00918 printMeshingOptions(mp, "Meshing options for boundary volume"); 00919 00920 00922 cout << "Start Volume meshing of whole model...." << endl; 00923 ngSurface = Ng_GenerateVolumeMesh (bSurface, &mp); 00924 if(ngSurface != NG_OK) { 00925 cout << "Error in Volume Meshing....Aborting!!" << endl; 00926 exit(1); 00927 } 00928 } else { 00929 bSurface = ngVolume; 00930 Ng_SetProperties(ngVolume, 2, 2, 2, 1); 00931 00932 } 00933 00934 00935 if (ini.get_optional<string>("meshing.second_order_surface").get_value_or("no") == "yes"){ 00936 cout << "Last second order meshing" << endl; 00937 Ng_Generate_SecondOrder(bSurface); 00938 } 00939 00940 cout << "Meshing successfully completed....!!" << endl; 00941 00942 // volume mesh output 00943 np = Ng_GetNP(bSurface); 00944 cout << "Points: " << np << endl; 00945 00946 ne = Ng_GetNE(bSurface); 00947 cout << "Elements: " << ne << endl; 00948 00949 string volumeVol = ini.get<string>("model.volume_vol"); 00950 00952 if (fname != "") 00953 volumeVol = fname; 00954 00955 if (volumeVol != ""){ 00956 cout << "Saving Mesh in VOL Format...." << endl; 00957 Ng_SaveMesh(bSurface,volumeVol.c_str()); 00958 } 00959 } 00960 00961 return 1; 00962 00963 } 00964 00965 00977 void clean(mMesh &mSurface){ 00978 // some cleaning to get rid of bad file formats like stl that duplicate vertexes.. 00979 int dup = tri::Clean<mMesh>::RemoveDuplicateVertex(mSurface); 00980 int unref = tri::Clean<mMesh>::RemoveUnreferencedVertex(mSurface); 00981 printf("Removed %i duplicate and %i unreferenced vertices from mesh\n",dup,unref); 00982 tri::UpdateTopology<mMesh>::VertexFace(mSurface); 00983 } 00984 01001 void smooth(mMesh &mSurface, INI &ini){ 01002 01003 string line = ini.get<string>("model.smoothing"); 01004 istringstream ss(line); 01005 string mode = ""; 01006 unsigned int steps = 0; 01007 string progress = "*"; 01008 01009 while (ss >> mode >> steps){ 01010 if (mode == "t"){ 01011 cout << "taubin smoothing steps: " << steps << endl; 01012 for (unsigned int i = 0; i < steps; i++){ 01013 cout << "\r" << i+1 << "/" << steps << flush; 01014 tri::UpdateNormal<mMesh>::PerFace(mSurface); 01015 tri::Smooth<mMesh>::VertexCoordTaubin(mSurface,1,TAUBIN_LAMBDA,TAUBIN_MU); 01016 } 01017 } else if (mode == "lap"){ 01018 cout << "laplace smoothing steps: " << steps << endl; 01019 for (unsigned int i = 0; i < steps; i++){ 01020 cout << "\r" << i+1 << "/" << steps << flush; 01021 tri::UpdateNormal<mMesh>::PerFace(mSurface); 01022 tri::Smooth<mMesh>::VertexCoordLaplacian(mSurface,1); 01023 } 01024 } else if (mode == "hc"){ 01025 cout << "hc laplace smoothing steps: " << steps << endl; 01026 for (unsigned int i = 0; i < steps; i++){ 01027 cout << "\r" << i+1 << "/" << steps << flush; 01028 tri::UpdateNormal<mMesh>::PerFace(mSurface); 01029 tri::Smooth<mMesh>::VertexCoordLaplacianHC(mSurface,1); 01030 } 01031 } else if (mode == "aw"){ 01032 cout << "angle weighted laplace smoothing steps: " << steps << endl; 01033 for (unsigned int i = 0; i < steps; i++){ 01034 cout << "\r" << i+1 << "/" << steps << flush; 01035 tri::UpdateNormal<mMesh>::PerFace(mSurface); 01036 tri::Smooth<mMesh>::VertexCoordLaplacianAngleWeighted(mSurface,1, CF_LAMBDA); 01037 } 01038 } else if (mode == "lapsd"){ 01039 cout << "scale dependent laplace smoothing steps: " << steps << endl; 01040 for (unsigned int i = 0; i < steps; i++){ 01041 cout << "\r" << i+1 << "/" << steps << flush; 01042 tri::UpdateNormal<mMesh>::PerFace(mSurface); 01043 // tri::UpdateNormal<mMesh>::PerFaceNormalized(mSurface); 01044 tri::Smooth<mMesh>::VertexCoordScaleDependentLaplacian_Fujiwara(mSurface,1, 0.0025); 01045 } 01046 } 01047 cout << endl; 01048 01049 } 01050 01051 // Example code to test for self intersections, Abort here if true? 01052 // vector<mMesh::FaceType *> SelfIntersectList; 01053 // tri::Clean<mMesh>::SelfIntersections(mSurface, SelfIntersectList); 01054 // int SelfIntersections = SelfIntersectList.size(); 01055 // cout << SelfIntersections << " intersections found." << endl; 01056 // if (SelfIntersections > 0) 01057 // cout << "This surface will probably not mesh! Keep care." << endl; 01058 01059 } 01060 01061 Ng_Result ngSurface; 01062 Ng_Mesh *ngVolume; 01063 01064 }; 01065 01066 bool Model::exclusionSurface = false;
1.7.6.1