MeshKit
1.0
|
00001 #include "meshkit/NGTetMesher.hpp" 00002 #include "meshkit/MeshScheme.hpp" 00003 #include "meshkit/ModelEnt.hpp" 00004 #include "moab/Interface.hpp" 00005 #include "moab/ReadUtilIface.hpp" 00006 #include "moab/Range.hpp" 00007 #include "meshkit/RegisterMeshOp.hpp" 00008 namespace nglib 00009 { 00010 #include "nglib.h" 00011 } 00012 00013 #include <vector> 00014 00015 namespace MeshKit 00016 { 00017 00018 moab::EntityType NGTetMesher::meshTps[] = {moab::MBVERTEX, moab::MBTET, moab::MBMAXTYPE}; 00019 00020 NGTetMesher::NGTetMesher(MKCore *mk_core, const MEntVector &me_vec) 00021 : MeshScheme(mk_core, me_vec) 00022 { 00023 } 00024 00025 NGTetMesher::~NGTetMesher() 00026 { 00027 } 00028 00029 MeshOp *NGTetMesher::get_tri_mesher() 00030 { 00031 std::vector<MeshOpProxy *> proxies; 00032 // mk_core()->meshop_by_mesh_type(moab::MBTRI, proxies); 00033 // if (proxies.empty()) throw Error(MK_FAILURE, "Couldn't find a MeshOp capable of producing triangles."); 00034 //return mk_core()->construct_meshop(*proxies.begin()); 00035 return mk_core()->construct_meshop("CAMALTriAdvance"); 00036 } 00037 00038 void NGTetMesher::setup_this() 00039 { 00040 MeshOp *tri_mesher = NULL; 00041 MEntVector surfs; 00042 00043 for (MEntSelection::iterator sit = mentSelection.begin(); sit != mentSelection.end(); sit++) { 00044 ModelEnt *me = (*sit).first; 00045 if (me->dimension() != 3) throw Error(MK_BAD_INPUT, "Tet mesher assigned to an entity with dimension != 3."); 00046 00047 // get the bounding faces 00048 surfs.clear(); 00049 me->get_adjacencies(2, surfs); 00050 00051 // check the mesh scheme on each one; if there is one, verify it can generate tris; if there 00052 // isn't one, make one 00053 bool inserted = false; 00054 00055 for (MEntVector::iterator vit = surfs.begin(); vit != surfs.end(); vit++) { 00056 if ((*vit)->is_meshops_list_empty()) { 00057 // get a tri mesher if we haven't already 00058 if (!tri_mesher) tri_mesher = get_tri_mesher(); 00059 00060 // add this surface to it, and if first for the volume, make sure it's added upstream 00061 tri_mesher->add_modelent(*vit); 00062 if (!inserted) { 00063 mk_core()->insert_node(tri_mesher, this); 00064 inserted = true; 00065 } 00066 } // if no meshops 00067 } // over surfs 00068 } // over vols 00069 } 00070 00071 void NGTetMesher::execute_this() 00072 { 00073 nglib::Ng_Init(); 00074 00075 for (MEntSelection::iterator sit = mentSelection.begin(); sit != mentSelection.end(); sit++) { 00076 // make a me, for convenience 00077 ModelEnt *me = (*sit).first; 00078 00079 // assemble bounding mesh 00080 std::vector<moab::EntityHandle> bdy; 00081 std::vector<int> group_sizes, senses, bdy_ids; 00082 me->boundary(2, bdy, &senses, &group_sizes); 00083 00084 // get connectivity 00085 00086 // convert the handles to integers for input to CMLTetMesher 00087 moab::Range bdy_vrange; 00088 std::vector<double> coords; 00089 me->get_indexed_connect_coords(bdy, &senses, NULL, bdy_ids, coords, &bdy_vrange, 1); 00090 00091 nglib::Ng_Mesh *ngmesh = nglib::Ng_NewMesh(); 00092 unsigned int numvs = bdy_vrange.size(); 00093 00094 // add the points 00095 for (unsigned int i = 0; i < numvs; i++) 00096 nglib::Ng_AddPoint(ngmesh, &coords[3*i]); 00097 00098 // add the tris 00099 unsigned int numtris = bdy.size(); 00100 for (unsigned int i = 0; i < numtris; i++) 00101 nglib::Ng_AddSurfaceElement(ngmesh, nglib::NG_TRIG, &bdy_ids[3*i]); 00102 00103 double my_size = me->mesh_interval_size(); 00104 if (0 > my_size) my_size = 1.0; 00105 nglib::Ng_RestrictMeshSizeGlobal(ngmesh, my_size); 00106 00107 nglib::Ng_Meshing_Parameters ngp; 00108 ngp.maxh = my_size; 00109 ngp.fineness = 0.5; 00110 00111 // Use variable name second_order if using newer version of NetGen 00112 ngp.secondorder = 0; 00113 // ngp.second_order = 0; 00114 // nglib.h Rev. 663 http://netgen-mesher.svn.sourceforge.net/viewvc/netgen-mesher/netgen/nglib/nglib.h?revision=663&view=markup 00115 //uses: second_order instead of secondorder(4.9.13) 00116 00117 00118 nglib::Ng_Result result = nglib::Ng_GenerateVolumeMesh(ngmesh, &ngp); 00119 if (nglib::NG_OK != result) ECERRCHK(MK_FAILURE, "Netgen mesher returned !ok."); 00120 00121 moab::ErrorCode rval; 00122 moab::Range new_ents; 00123 int num_pts = nglib::Ng_GetNP(ngmesh); 00124 assert(num_pts >= (int)numvs); 00125 if (num_pts > (int)numvs) { 00126 coords.resize(3*(num_pts-numvs)); 00127 00128 for (unsigned int i = (unsigned int) numvs; i < (unsigned int) num_pts; i++) 00129 nglib::Ng_GetPoint(ngmesh, i, &coords[3*(i-numvs)]); 00130 00131 // create the new vertices' entities 00132 rval = mk_core()->moab_instance()->create_vertices(&coords[0], num_pts-numvs, new_ents); 00133 MBERRCHK(rval, mk_core()->moab_instance()); 00134 } 00135 00136 // for tets, pre-allocate connectivity 00137 moab::ReadUtilIface *iface; 00138 rval = mk_core()-> moab_instance() -> query_interface(iface); 00139 MBERRCHK(rval, mk_core()->moab_instance()); 00140 00141 //create the tris, get a direct ptr to connectivity 00142 int num_tets = (int) nglib::Ng_GetNE(ngmesh); 00143 moab::EntityHandle starth, *connect; 00144 rval = iface->get_element_connect(num_tets, 4, moab::MBTET, 1, starth, connect); 00145 MBERRCHK(rval, mk_core()->moab_instance()); 00146 00147 // read connectivity directly into that array, as int's 00148 int *connecti = (int*) connect; 00149 for (unsigned int i = 0; i < (unsigned int) num_tets; i++) 00150 nglib::Ng_GetVolumeElement(ngmesh, i+1, connecti+4*i); 00151 00152 // put vertex handles into an indexible array 00153 std::vector<moab::EntityHandle> bdy_vvec; 00154 std::copy(bdy_vrange.begin(), bdy_vrange.end(), std::back_inserter(bdy_vvec)); 00155 std::copy(new_ents.begin(), new_ents.end(), std::back_inserter(bdy_vvec)); 00156 00157 // now convert vertex indices into handles in-place, working from the back 00158 for (int i = 4*num_tets-1; i >= 0; i--) { 00159 assert(connecti[i] > 0 && connecti[i] <= (int)bdy_vvec.size()); 00160 connect[i] = bdy_vvec[connecti[i]-1]; 00161 } 00162 00163 // put new tris into new entity range, then commit the mesh 00164 new_ents.insert(starth, starth+num_tets-1); 00165 me->commit_mesh(new_ents, COMPLETE_MESH); 00166 00167 nglib::Ng_DeleteMesh(ngmesh); 00168 } 00169 00170 nglib::Ng_Exit(); 00171 } 00172 00173 } // namespace MeshKit