MeshKit
1.0
|
00001 00007 #include "meshkit/MKCore.hpp" 00008 #include "meshkit/EdgeMesher.hpp" 00009 #include "meshkit/SizingFunction.hpp" 00010 #include "meshkit/SizingFunctionVar.hpp" 00011 #include "meshkit/ModelEnt.hpp" 00012 #include "meshkit/Matrix.hpp" 00013 #include "meshkit/VertexMesher.hpp" 00014 00015 using namespace MeshKit; 00016 00017 #include "TestUtil.hpp" 00018 00019 MKCore *mk = NULL; 00020 00021 void edgemesh_hole(); 00022 void edgemesh_square(); 00023 void edgemesh_brick(); 00024 void edgemesh_var(); 00025 void edgemesh_dual(); 00026 00027 #ifdef HAVE_ACIS 00028 std::string extension = ".sat"; 00029 #elif HAVE_OCC 00030 std::string extension = ".stp"; 00031 #endif 00032 00033 int main(int argc, char **argv) 00034 { 00035 00036 // start up MK and load the geometry 00037 mk = new MKCore(); 00038 00039 int num_fail = 0; 00040 num_fail += RUN_TEST(edgemesh_hole); 00041 num_fail += RUN_TEST(edgemesh_square); 00042 num_fail += RUN_TEST(edgemesh_brick); 00043 num_fail += RUN_TEST(edgemesh_var); 00044 num_fail += RUN_TEST(edgemesh_dual); 00045 #if HAVE_OCC 00046 return 0; 00047 #else 00048 return num_fail; 00049 #endif 00050 } 00051 00052 void edgemesh_hole() 00053 { 00054 std::string file_name = TestDir + "/holysurf" + extension; 00055 mk->load_geometry(file_name.c_str()); 00056 00057 // get the surface 00058 MEntVector surfs, curves, loops; 00059 mk->get_entities_by_dimension(2, surfs); 00060 CHECK_EQUAL(1, (int)surfs.size()); 00061 00062 // test getting loops 00063 std::vector<int> senses, loop_sizes; 00064 surfs[0]->boundary(0, loops, &senses, &loop_sizes); 00065 CHECK_EQUAL(4, (int)loop_sizes.size()); 00066 00067 // add up the loop sizes, should add to 8 00068 unsigned int l = 0; 00069 for (unsigned int i = 0; i < loop_sizes.size(); i++) l += loop_sizes[i]; 00070 CHECK_EQUAL(8, (int)l); 00071 00072 // make an edge mesher 00073 mk->get_entities_by_dimension(1, curves); 00074 EdgeMesher *em = (EdgeMesher*) mk->construct_meshop("EdgeMesher", curves); 00075 00076 // make a sizing function and set it on the surface 00077 SizingFunction esize(mk, -1, 0.25); 00078 surfs[0]->sizing_function_index(esize.core_index()); 00079 00080 // mesh the edges, by calling execute 00081 mk->setup_and_execute(); 00082 00083 // make sure we got the right number of edges 00084 moab::Range edges; 00085 moab::ErrorCode rval = mk->moab_instance()->get_entities_by_dimension(0, 1, edges); 00086 CHECK_EQUAL(moab::MB_SUCCESS, rval); 00087 CHECK_EQUAL(79, (int)edges.size()); 00088 00089 // clean up 00090 delete em; 00091 delete mk->vertex_mesher(); 00092 } 00093 00094 void edgemesh_square() 00095 { 00096 std::string file_name = TestDir + "/squaresurf" + extension; 00097 mk->load_geometry(file_name.c_str()); 00098 00099 // get the surface 00100 MEntVector surfs, curves, loops; 00101 mk->get_entities_by_dimension(2, surfs); 00102 ModelEnt *this_surf = (*surfs.rbegin()); 00103 00104 // test getting loops 00105 std::vector<int> senses, loop_sizes; 00106 this_surf->boundary(0, loops, &senses, &loop_sizes); 00107 CHECK_EQUAL(1, (int)loop_sizes.size()); 00108 00109 // make an edge mesher 00110 this_surf->get_adjacencies(1, curves); 00111 mk->construct_meshop("EdgeMesher", curves); 00112 00113 // make a sizing function and set it on the surface 00114 SizingFunction esize(mk, 1, 1.0); 00115 this_surf->sizing_function_index(esize.core_index()); 00116 00117 // mesh the edges, by calling execute 00118 mk->setup_and_execute(); 00119 00120 senses.clear(); 00121 loop_sizes.clear(); 00122 std::vector<moab::EntityHandle> mloops; 00123 this_surf->boundary(0, mloops, &senses, &loop_sizes); 00124 CHECK_EQUAL(4, (int)mloops.size()); 00125 00126 // print the vertex positions 00127 std::vector<double> coords(12); 00128 moab::ErrorCode rval = mk->moab_instance()->get_coords(&mloops[0], mloops.size(), &coords[0]); 00129 MBERRCHK(rval, mk->moab_instance()); 00130 // std::cout << "Vertex positions:" << std::endl; 00131 // for (unsigned int i = 0; i < 4; i++) 00132 // std::cout << coords[3*i] << ", " << coords[3*i+1] << ", " << coords[3*i+2] << std::endl; 00133 00134 // compute the cross product vector and print that 00135 double tmp[] = {0.0, 0.0, 1.0}; 00136 Vector<3> p0(&coords[0]), p1(&coords[3]), p2(&coords[6]), unitz(tmp); 00137 p0 -= p1; 00138 p2 -= p1; 00139 p1 = vector_product(p2, p0); 00140 double dot = inner_product(unitz, p1); 00141 CHECK_REAL_EQUAL(1.0, dot, 1.0e-6); 00142 00143 mk->clear_graph(); 00144 } 00145 00146 void edgemesh_brick() 00147 { 00148 std::string file_name = TestDir + "/brick" + extension; 00149 mk->load_geometry(file_name.c_str()); 00150 00151 // get the vol, surfs, curves 00152 MEntVector vols, surfs, curves, loops; 00153 mk->get_entities_by_dimension(3, vols); 00154 ModelEnt *this_vol = (*vols.rbegin()); 00155 this_vol->get_adjacencies(2, surfs); 00156 this_vol->get_adjacencies(1, curves); 00157 00158 // mesh all the edges 00159 mk->construct_meshop("EdgeMesher", curves); 00160 SizingFunction esize(mk, 1, 1.0); 00161 this_vol->sizing_function_index(esize.core_index()); 00162 mk->setup_and_execute(); 00163 00164 // test loop directionality for each surface 00165 for (MEntVector::iterator vit = surfs.begin(); vit != surfs.end(); vit++) { 00166 00167 std::vector<int> senses, loop_sizes; 00168 std::vector<moab::EntityHandle> mloops; 00169 (*vit)->boundary(0, mloops, &senses, &loop_sizes); 00170 CHECK_EQUAL(4, (int)mloops.size()); 00171 00172 // get the vertex positions 00173 std::vector<double> coords(12); 00174 moab::ErrorCode rval = mk->moab_instance()->get_coords(&mloops[0], mloops.size(), &coords[0]); 00175 MBERRCHK(rval, mk->moab_instance()); 00176 00177 // compute the cross product vector and compare it to the surface normal 00178 Vector<3> p0(&coords[0]), p1(&coords[3]), p2(&coords[6]), p3; 00179 p0 -= p1; 00180 p2 -= p1; 00181 p1 = vector_product(p2, p0); 00182 (*vit)->evaluate(coords[0], coords[1], coords[2], NULL, p3.data()); 00183 double dot = inner_product(p3, p1); 00184 CHECK_REAL_EQUAL(1.0, dot, 1.0e-6); 00185 } 00186 00187 mk->clear_graph(); 00188 } 00189 void edgemesh_var() 00190 { 00191 // delete existing mesh 00192 mk->moab_instance()->delete_mesh(); 00193 00194 std::string file_name = TestDir + "/squaresurf" + extension; 00195 mk->load_geometry(file_name.c_str()); 00196 00197 // get the surface 00198 MEntVector surfs, curves, loops; 00199 mk->get_entities_by_dimension(2, surfs); 00200 ModelEnt *this_surf = (*surfs.rbegin()); 00201 00202 // test getting loops 00203 std::vector<int> senses, loop_sizes; 00204 this_surf->boundary(0, loops, &senses, &loop_sizes); 00205 CHECK_EQUAL(1, (int)loop_sizes.size()); 00206 00207 // make an edge mesher 00208 this_surf->get_adjacencies(1, curves); 00209 EdgeMesher * emesher = (EdgeMesher *)mk->construct_meshop("EdgeMesher", curves); 00210 00211 emesher->set_edge_scheme(EdgeMesher::VARIABLE); 00212 00213 // make a sizing function and set it on the surface 00214 SizingFunctionVar * svar = new SizingFunctionVar(mk, -1, 0.1); 00215 00216 // these could be read from a file, or something 00217 double point0[3] = {0, 0, 0}; 00218 double coeffs[4] = {0.05, 0.05, 0.05, 0.1}; 00219 svar->set_linear_coeff(point0, coeffs); 00220 00221 this_surf->sizing_function_index(svar->core_index()); 00222 00223 mk->setup_and_execute(); 00224 00225 //mk->moab_instance()->write_file("var1.vtk"); 00226 mk->clear_graph(); 00227 00228 } 00229 void edgemesh_dual() 00230 { 00231 // delete existing mesh 00232 mk->moab_instance()->delete_mesh(); 00233 00234 std::string file_name = TestDir + "/squaresurf" + extension; 00235 mk->load_geometry(file_name.c_str()); 00236 00237 // get the surface 00238 MEntVector surfs, curves, loops; 00239 mk->get_entities_by_dimension(2, surfs); 00240 ModelEnt *this_surf = (*surfs.rbegin()); 00241 00242 // test getting loops 00243 std::vector<int> senses, loop_sizes; 00244 this_surf->boundary(0, loops, &senses, &loop_sizes); 00245 CHECK_EQUAL(1, (int)loop_sizes.size()); 00246 00247 // make an edge mesher 00248 this_surf->get_adjacencies(1, curves); 00249 EdgeMesher * emesher = (EdgeMesher *)mk->construct_meshop("EdgeMesher", curves); 00250 00251 emesher->set_edge_scheme(EdgeMesher::DUAL); 00252 00253 // make a sizing function and set it on the surface 00254 SizingFunction * svar = new SizingFunction(mk, 10, -1); 00255 00256 emesher->set_ratio(1.3); 00257 00258 this_surf->sizing_function_index(svar->core_index()); 00259 00260 mk->setup_and_execute(); 00261 00262 //mk->moab_instance()->write_file("dual1.1.vtk"); 00263 mk->clear_graph(); 00264 00265 }