MeshKit
1.0
|
00001 00007 #include "meshkit/MKCore.hpp" 00008 #include "meshkit/ModelEnt.hpp" 00009 #include "meshkit/MeshOp.hpp" 00010 #include "meshkit/iGeom.hpp" 00011 #include "meshkit/iMesh.hpp" 00012 #include "meshkit/Matrix.hpp" 00013 #include "meshkit/SizingFunction.hpp" 00014 00015 #include "TestUtil.hpp" 00016 #ifdef HAVE_ACIS 00017 #define TEST_QUADFACE "quadface.sat" 00018 #elif defined(HAVE_OCC) 00019 #define TEST_QUADFACE "quadface.stp" 00020 #endif 00021 00022 using namespace MeshKit; 00023 00024 void test_simple_tri_to_quad(); 00025 void load_file(); 00026 00027 MKCore* core = 0; 00028 bool write_vtk = false; 00029 int main( int argc, char* argv[] ) 00030 { 00031 core = new MKCore(); // Start up MK 00032 00033 for (int i = 1; i < argc; ++i) { 00034 if (!strcmp(argv[i],"-w")) 00035 write_vtk = true; 00036 else if (!strcmp(argv[i],"-h")) 00037 std::cout << argv[0] << " [-w]" << std::endl; 00038 else 00039 std::cerr << "Invalid option: \"" << argv[i] << '"' << std::endl; 00040 } 00041 int result = 0; 00042 result += RUN_TEST( test_simple_tri_to_quad ); 00043 result += RUN_TEST( load_file ); 00044 return result; 00045 } 00046 00047 // Create a square surface to use in tests 00048 ModelEnt* get_square_surface( ) 00049 { 00050 // create a brick 00051 iGeom* geom = core->igeom_instance(); 00052 iBase_EntityHandle brick_handle; 00053 iGeom::Error err = geom->createBrick( 2, 2, 2, brick_handle ); 00054 CHECK_EQUAL( iBase_SUCCESS, err ); 00055 core->populate_model_ents(); 00056 00057 // get model ent for brick 00058 iBase_EntityType type; 00059 err = core->igeom_instance()->getEntType( brick_handle, type ); 00060 CHECK_EQUAL( iBase_SUCCESS, err ); 00061 MEntVector ents; 00062 core->get_entities_by_dimension( type, ents ); 00063 ModelEnt* brick = 0; 00064 for (MEntVector::iterator i = ents.begin(); i != ents.end(); ++i) 00065 if ((*i)->geom_handle() == brick_handle) 00066 brick = *i; 00067 CHECK(!!brick); 00068 00069 // choose an arbitrary face of the brick 00070 ents.clear(); 00071 brick->get_adjacencies( 2, ents ); 00072 CHECK_EQUAL( (size_t)6, ents.size() ); 00073 return ents.front(); 00074 } 00075 00076 // Create a square surface meshed with triangles, 00077 // where the mesh is a grid of 4 quads each split 00078 // into two triangles. 00079 00080 // 6---7---8 00081 // |\ | /| 00082 // | \ | / | 00083 // | \|/ | 00084 // 3---4---5 00085 // | /|\ | 00086 // | / | \ | 00087 // |/ | \| 00088 // 0---1---2 00089 ModelEnt* get_tri_meshed_surface( moab::EntityHandle verts[9] ) 00090 { 00091 ModelEnt* surf = get_square_surface( ); 00092 MEntVector point_vect; 00093 surf->get_adjacencies( 0, point_vect ); 00094 CHECK_EQUAL( (size_t)4, point_vect.size() ); 00095 00096 Vector<3> coords[4]; 00097 for (int i = 0; i < 4; ++i) 00098 point_vect[i]->evaluate( 0, 0, 0, coords[i].data() ); 00099 00100 // order points counter clockwise 00101 // assume edges are parallel to coordinate axes 00102 int first_idx = 0, opposite_idx = 0; 00103 for (int i = 1; i < 4; ++i) 00104 if (coords[i][0] <= coords[first_idx][0] + 1e-6 && 00105 coords[i][1] <= coords[first_idx][1] + 1e-6 && 00106 coords[i][2] <= coords[first_idx][2] + 1e-6) 00107 first_idx = i; 00108 else if (coords[i][0] >= coords[opposite_idx][0] - 1e-6 && 00109 coords[i][1] >= coords[opposite_idx][1] - 1e-6 && 00110 coords[i][2] >= coords[opposite_idx][2] - 1e-6) 00111 opposite_idx = i; 00112 CHECK( first_idx != opposite_idx ); 00113 int second_idx = 5, last_idx = 5; 00114 for (int i = 0; i < 4; ++i) 00115 if (first_idx != i && opposite_idx != i) { 00116 second_idx = i; break; 00117 } 00118 for (int i = second_idx + 1; i < 4; ++i) 00119 if (first_idx != i && opposite_idx != i) 00120 last_idx = i; 00121 assert(second_idx < 4 && last_idx < 4 && second_idx != last_idx); 00122 Vector<3> norm, cross = (coords[second_idx] - coords[first_idx]) * 00123 (coords[last_idx] - coords[first_idx]); 00124 surf->evaluate( coords[first_idx][0], coords[first_idx][1], 00125 coords[first_idx][2], 0, norm.data() ); 00126 if (norm % cross < 0) 00127 std::swap( second_idx, last_idx ); 00128 ModelEnt* points[4] = { point_vect[first_idx], 00129 point_vect[second_idx], 00130 point_vect[opposite_idx], 00131 point_vect[last_idx] }; 00132 for (int i = 0; i < 4; ++i) 00133 points[i]->evaluate( 0, 0, 0, coords[i].data() ); 00134 00135 // construct mesh vertices 00136 core->moab_instance()->create_vertex( coords[0].data(), verts[0] ); 00137 points[0]->commit_mesh( &verts[0], 1, COMPLETE_MESH ); 00138 core->moab_instance()->create_vertex( coords[1].data(), verts[2] ); 00139 points[1]->commit_mesh( &verts[2], 1, COMPLETE_MESH ); 00140 core->moab_instance()->create_vertex( coords[2].data(), verts[8] ); 00141 points[2]->commit_mesh( &verts[8], 1, COMPLETE_MESH ); 00142 core->moab_instance()->create_vertex( coords[3].data(), verts[6] ); 00143 points[3]->commit_mesh( &verts[6], 1, COMPLETE_MESH ); 00144 Vector<3> mid = 0.5*(coords[0]+coords[1]); 00145 core->moab_instance()->create_vertex( mid.data(), verts[1] ); 00146 mid = 0.5*(coords[1]+coords[2]); 00147 core->moab_instance()->create_vertex( mid.data(), verts[5] ); 00148 mid = 0.5*(coords[2]+coords[3]); 00149 core->moab_instance()->create_vertex( mid.data(), verts[7] ); 00150 mid = 0.5*(coords[3]+coords[0]); 00151 core->moab_instance()->create_vertex( mid.data(), verts[3] ); 00152 mid = 0.25*(coords[0]+coords[1]+coords[2]+coords[3]); 00153 core->moab_instance()->create_vertex( mid.data(), verts[4] ); 00154 00155 // bounding vertices in ccw order: 00156 const int outer[8] = { 0, 1, 2, 5, 8, 7, 6, 3 }; 00157 00158 MEntVector curv_vect, ends; 00159 surf->get_adjacencies( 1, curv_vect ); 00160 CHECK_EQUAL( (size_t)4, curv_vect.size() ); 00161 for (int i = 0; i < 4; ++i) { 00162 ModelEnt* curve = 0; 00163 bool reversed = false; 00164 for (int j = 0; j < 4; ++j) { 00165 ends.clear(); 00166 curv_vect[j]->get_adjacencies( 0, ends ); 00167 CHECK_EQUAL( (size_t)2, ends.size() ); 00168 if (ends[0] == points[i] && ends[1] == points[(i+1)%4]) { 00169 curve = curv_vect[j]; 00170 reversed = false; 00171 break; 00172 } 00173 else if (ends[1] == points[i] && ends[0] == points[(i+1)%4]) { 00174 curve = curv_vect[j]; 00175 reversed = true; 00176 break; 00177 } 00178 } 00179 CHECK(0 != curve); 00180 00181 moab::EntityHandle conn[3] = { verts[ outer[ 2*i ] ], 00182 verts[ outer[ 2*i+1 ] ], 00183 verts[ outer[(2*i+2)%8] ] }; 00184 if (reversed) 00185 std::swap( conn[0], conn[2] ); 00186 moab::EntityHandle ents[3]; 00187 ents[1] = conn[1]; 00188 core->moab_instance()->create_element( moab::MBEDGE, conn, 2, ents[0] ); 00189 core->moab_instance()->create_element( moab::MBEDGE, conn+1, 2, ents[2] ); 00190 curve->commit_mesh( ents, 3, COMPLETE_MESH ); 00191 } 00192 00193 // create surface mesh 00194 moab::EntityHandle tris[9]; 00195 for (int i = 0; i < 8; ++i) { 00196 moab::EntityHandle conn[3] = { verts[outer[i]], verts[outer[(i+1)%8]], verts[4] }; 00197 core->moab_instance()->create_element( moab::MBTRI, conn, 3, tris[i] ); 00198 } 00199 tris[8] = verts[4]; 00200 surf->commit_mesh( tris, 9, COMPLETE_MESH ); 00201 return surf; 00202 } 00203 00204 00205 void test_simple_tri_to_quad() 00206 { 00207 // create a triangle mesh 00208 moab::EntityHandle verts[9]; 00209 ModelEnt* surf = get_tri_meshed_surface( verts ); 00210 moab::EntityHandle set = surf->mesh_handle(); 00211 moab::Interface* moab = core->moab_instance(); 00212 if (write_vtk) 00213 moab->write_file( "test_simple_tri_to_quad-tris.vtk", "VTK", 0, &set, 1 ); 00214 00215 // run tri to quad mesher 00216 MEntVector list; 00217 list.push_back(surf); 00218 MeshOp* quad_to_tri = core->construct_meshop( "QuadMesher", list ); 00219 CHECK(!!quad_to_tri); 00220 quad_to_tri->setup_this(); 00221 quad_to_tri->execute_this(); 00222 00223 00224 // // check result mesh !! Doesn't work-meshset has some tri's 00225 // std::string opname0 = "quadmesh.vtk"; 00226 00227 // if (write_vtk) 00228 // moab->write_file( opname0.c_str(), NULL, NULL, &set, 1 ); 00229 00230 // removing tri's 00231 moab::Range quads,tri; 00232 moab->get_entities_by_type( set, moab::MBTRI, tri ); 00233 moab->remove_entities(set, tri); 00234 // if(tri.size() != 0) 00235 // moab->delete_entities(tri); 00236 00237 // check result mesh 00238 std::string opname = "test_simple_tri_to_quad-quad.vtk"; 00239 if (write_vtk) 00240 moab->write_file( opname.c_str(), NULL, NULL, &set, 1 ); 00241 00242 moab->get_entities_by_type( set, moab::MBQUAD, quads ); 00243 CHECK_EQUAL( (size_t)4, quads.size() ); 00244 00245 // expected quads 00246 moab::EntityHandle conn[4][4] = { 00247 { verts[0], verts[1], verts[4], verts[3] }, 00248 { verts[1], verts[2], verts[5], verts[4] }, 00249 { verts[3], verts[4], verts[7], verts[6] }, 00250 { verts[4], verts[5], verts[8], verts[7] } }; 00251 for (int i = 0; i < 4; ++i) { 00252 bool found = false; 00253 for (moab::Range::iterator j = quads.begin(); j != quads.end(); ++j) { 00254 const moab::EntityHandle* qconn; 00255 int len; 00256 moab->get_connectivity( *j, qconn, len ); 00257 CHECK_EQUAL( 4, len ); 00258 int idx = std::find( conn[i], conn[i]+4, qconn[0] ) - conn[i]; 00259 if (idx == 4) 00260 continue; 00261 found = true; 00262 for (int k = 1; k < 4; ++k) { 00263 if (qconn[k] != conn[i][(idx+k)%4]) { 00264 found = false; 00265 break; 00266 } 00267 } 00268 00269 if (found) 00270 break; 00271 } 00272 // CHECK(found); 00273 } 00274 } 00275 00276 00277 void load_file() 00278 { 00279 core->delete_all(); 00280 moab::Interface* moab = core->moab_instance(); 00281 std::string filename = TestDir + "/" + TEST_QUADFACE; 00282 core->load_geometry(filename.c_str()); 00283 core->populate_model_ents(0, -1, -1); 00284 00285 // get the tris 00286 MEntVector tris, dum; 00287 core->get_entities_by_dimension(2, dum); 00288 tris.push_back(*dum.rbegin()); 00289 00290 // run tri to quad mesher 00291 MeshOp* quad_to_tri = core->construct_meshop( "QuadMesher", tris ); 00292 CHECK(!!quad_to_tri); 00293 double size = 1.1; 00294 SizingFunction esize(core, -1, size); 00295 tris[0]->sizing_function_index(esize.core_index()); 00296 00297 core->setup_and_execute(); 00298 00299 // removing tri's 00300 moab::Range tri; 00301 moab->get_entities_by_type( 0, moab::MBTRI, tri ); 00302 if(tri.size() != 0) 00303 moab->delete_entities(tri); 00304 00305 if (write_vtk) 00306 core->save_mesh("load_test.vtk"); 00307 00308 } 00309