|
MOAB: Mesh Oriented datABase
(version 5.4.0)
|
|
MOAB: Mesh Oriented datABase
(version 5.4.0)
|
00001 /* ***************************************************************** 00002 MESQUITE -- The Mesh Quality Improvement Toolkit 00003 00004 Copyright 2004 Sandia Corporation and Argonne National 00005 Laboratory. Under the terms of Contract DE-AC04-94AL85000 00006 with Sandia Corporation, the U.S. Government retains certain 00007 rights in this software. 00008 00009 This library is free software; you can redistribute it and/or 00010 modify it under the terms of the GNU Lesser General Public 00011 License as published by the Free Software Foundation; either 00012 version 2.1 of the License, or (at your option) any later version. 00013 00014 This library is distributed in the hope that it will be useful, 00015 but WITHOUT ANY WARRANTY; without even the implied warranty of 00016 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00017 Lesser General Public License for more details. 00018 00019 You should have received a copy of the GNU Lesser General Public License 00020 (lgpl.txt) along with this library; if not, write to the Free Software 00021 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00022 00023 pknupp@sandia.gov, tleurent@mcs.anl.gov, tmunson@mcs.anl.gov 00024 00025 ***************************************************************** */ 00026 // -*- Mode : c++; tab-width: 2; c-tab-always-indent: t; indent-tabs-mode: nil; c-basic-offset: 2 00027 // -*- 00028 // 00029 // SUMMARY: 00030 // USAGE: 00031 // 00032 // ORIG-DATE: 13-08-12 by Boyd Tidwell 00033 // 00034 // 00035 // DESCRIPTION: 00036 // ============ 00037 /*! \file main.cpp 00038 00039 Benchmark timing tests for various MBMesquite wrappers. 00040 00041 */ 00042 // DESCRIP-END. 00043 // 00044 00045 #include <iostream> 00046 using std::cout; 00047 using std::endl; 00048 #include <cstdlib> 00049 00050 #include "Mesquite.hpp" 00051 #include "MsqError.hpp" 00052 #include "MeshImpl.hpp" 00053 #include "Vector3D.hpp" 00054 #include "InstructionQueue.hpp" 00055 #include "PatchData.hpp" 00056 #include "TerminationCriterion.hpp" 00057 #include "QualityAssessor.hpp" 00058 #include "PlanarDomain.hpp" 00059 #include "MsqTimer.hpp" 00060 00061 // algorythms 00062 #include "LInfTemplate.hpp" 00063 #include "EdgeLengthQualityMetric.hpp" 00064 #include "LaplaceWrapper.hpp" 00065 #include "UntangleWrapper.hpp" 00066 00067 #include "ShapeImprover.hpp" 00068 #include "MsqTimer.hpp" 00069 #include "SizeAdaptShapeWrapper.hpp" 00070 #include "SphericalDomain.hpp" 00071 #include "PaverMinEdgeLengthWrapper.hpp" 00072 #include "DeformingDomainWrapper.hpp" 00073 #include "MeshDomain1D.hpp" 00074 #include "TestUtil.hpp" 00075 00076 using namespace MBMesquite; 00077 00078 std::string shape_improv_file_name_1 = 00079 TestDir + "unittest/mesquite/3D/vtk/hexes/untangled/1000hex-block-internal-bias.vtk"; 00080 std::string shape_improv_file_name_2 = TestDir + "unittest/mesquite/3D/vtk/tets/untangled/tire.vtk"; 00081 std::string laplacian_file_name_1 = TestDir + "unittest/mesquite/2D/vtk/quads/untangled/square_quad_10_rand.vtk"; 00082 std::string laplacian_file_name_2 = TestDir + "unittest/mesquite/2D/vtk//quads/untangled/shashkov_quad.vtk"; 00083 std::string untangle_file_name_1 = TestDir + "unittest/mesquite/2D/vtk/quads/tangled/tangled_horse1.vtk"; 00084 std::string untangle_file_name_2 = TestDir + "unittest/mesquite/2D/vtk/quads/untangled//shest_grid32.vtk"; 00085 std::string size_adapt_shape_file_name_1 = TestDir + "unittest/mesquite/2D/vtk/quads/untangled/bias-sphere-quads.vtk"; 00086 std::string min_edge_length_file_name_1 = TestDir + "unittest/mesquite/2D/vtk/quads/untangled/quads_4by2_bad.vtk"; 00087 std::string min_edge_length_file_name_2 = TestDir + "unittest/mesquite/2D/vtk/quads/untangled/shashkov_quad.vtk"; 00088 std::string deforming_domain_file_name_1 = TestDir + "unittest/mesquite/3D/vtk/hexes/untangled/sph-10-zsquare.vtk"; 00089 00090 void help( const char* ) 00091 { 00092 std::cerr << "Parameters are not supported for this test" << std::endl; 00093 exit( 1 ); 00094 } 00095 typedef std::vector< Mesh::ElementHandle > elem_vec_t; 00096 00097 // Assume mesh is on a sphere of radius 10 with an axis equal to Z. 00098 // Return elements near poles and elements near equator. 00099 void find_z10_extreme_elements( Mesh& mesh, elem_vec_t& polar_elems, elem_vec_t& equatorial_elems, MsqError& err ); 00100 00101 // Get areas of quads and tris 00102 void elem_areas( Mesh& mesh, const elem_vec_t& elems, double& min, double& mean, double& max, MsqError& err ); 00103 00104 void classify_boundary( Mesh* mesh, 00105 Mesh::VertexHandle corners_out[4], 00106 std::vector< Mesh::VertexHandle > curves_out[4], 00107 MsqError& err ); 00108 const double Z = 7.0; 00109 // size of new domain 00110 const double HD = 4.5; 00111 00112 int main( int, char*[] ) 00113 { 00114 std::cout << std::endl << "********* Wrappers Timing Tests **********" << std::endl << std::endl; 00115 00116 MBMesquite::MsqPrintError err( cout ); 00117 MBMesquite::MeshImpl mesh; 00118 00119 // #################### Begin ShapeImprover tests ################### 00120 00121 ShapeImprover si_wrapper; 00122 mesh.read_vtk( shape_improv_file_name_1.c_str(), err ); 00123 00124 Timer t; 00125 si_wrapper.run_instructions( &mesh, err ); 00126 if( err ) return 1; 00127 double si_s_secs = t.since_birth(); 00128 std::cout << std::endl 00129 << "ShapeImprover small file optimization completed in " << si_s_secs << " seconds" << std::endl; 00130 00131 mesh.clear(); 00132 mesh.read_vtk( shape_improv_file_name_2.c_str(), err ); 00133 00134 t.reset(); 00135 si_wrapper.run_instructions( &mesh, err ); 00136 if( err ) return 1; 00137 double si_l_secs = t.since_birth(); 00138 std::cout << std::endl 00139 << "ShapeImprover large file optimization completed in " << si_l_secs << " seconds" << std::endl; 00140 00141 // #################### Begin LaplacianWrapper tests ################### 00142 00143 Vector3D pnt1( 0, 0, 5 ); 00144 Vector3D s_norm( 0, 0, 1 ); 00145 MBMesquite::PlanarDomain msq_geom( s_norm, pnt1 ); 00146 00147 LaplaceWrapper lp_wrapper; 00148 00149 mesh.clear(); 00150 mesh.read_vtk( laplacian_file_name_1.c_str(), err ); 00151 if( err ) return 1; 00152 00153 MeshDomainAssoc mesh_and_domain4 = MeshDomainAssoc( &mesh, &msq_geom ); 00154 t.reset(); 00155 lp_wrapper.run_instructions( &mesh_and_domain4, err ); 00156 if( err ) return 1; 00157 double lp_s_secs = t.since_birth(); 00158 std::cout << std::endl 00159 << "LaplacianWrapper small file optimization completed in " << lp_s_secs << " seconds" << std::endl; 00160 00161 Vector3D pnt2( 0, 0, 0 ); 00162 MBMesquite::PlanarDomain msq_geom2( s_norm, pnt2 ); 00163 00164 mesh.clear(); 00165 mesh.read_vtk( laplacian_file_name_2.c_str(), err ); 00166 if( err ) return 1; 00167 00168 MeshDomainAssoc mesh_and_domain5 = MeshDomainAssoc( &mesh, &msq_geom2 ); 00169 t.reset(); 00170 lp_wrapper.run_instructions( &mesh_and_domain5, err ); 00171 if( err ) return 1; 00172 double lp_l1_secs = t.since_birth(); 00173 std::cout << std::endl 00174 << "LaplacianWrapper large file (term crit=0.001) completed in " << lp_l1_secs << " seconds" << std::endl; 00175 00176 mesh.clear(); 00177 mesh.read_vtk( laplacian_file_name_2.c_str(), err ); 00178 if( err ) return 1; 00179 00180 lp_wrapper.set_vertex_movement_limit_factor( 0.1 ); 00181 t.reset(); 00182 lp_wrapper.run_instructions( &mesh_and_domain5, err ); 00183 if( err ) return 1; 00184 double lp_l2_secs = t.since_birth(); 00185 std::cout << std::endl 00186 << "LaplacianWrapper large file (term crit=0.1) completed in " << lp_l2_secs << " seconds" << std::endl; 00187 00188 // #################### Begin UntangleWrapper::BETA tests ################### 00189 00190 mesh.clear(); 00191 mesh.read_vtk( untangle_file_name_1.c_str(), err ); 00192 if( err ) return 1; 00193 00194 std::vector< Mesh::VertexHandle > verts; 00195 mesh.get_all_vertices( verts, err ); 00196 if( err || verts.empty() ) return 1; 00197 MsqVertex coords; 00198 mesh.vertices_get_coordinates( arrptr( verts ), &coords, 1, err ); 00199 if( err ) return 1; 00200 Vector3D norm( 0, 0, 1 ); 00201 PlanarDomain u_domain( norm, coords ); 00202 00203 UntangleWrapper::UntangleMetric metric = UntangleWrapper::BETA; 00204 UntangleWrapper un_wrapper( metric ); 00205 un_wrapper.set_vertex_movement_limit_factor( 0.005 ); 00206 00207 MeshDomainAssoc mesh_and_domain3 = MeshDomainAssoc( &mesh, &u_domain ); 00208 t.reset(); 00209 un_wrapper.run_instructions( &mesh_and_domain3, err ); 00210 if( err ) return 1; 00211 00212 double unb_s_secs = t.since_birth(); 00213 std::cout << std::endl 00214 << "UntangleWrapper::BETA small file optimization completed in " << unb_s_secs << " seconds" << std::endl; 00215 00216 mesh.clear(); 00217 mesh.read_vtk( untangle_file_name_2.c_str(), err ); 00218 if( err ) return 1; 00219 00220 // get domain 00221 verts.clear(); 00222 mesh.get_all_vertices( verts, err ); 00223 if( err || verts.empty() ) return 1; 00224 MsqVertex coords2; 00225 mesh.vertices_get_coordinates( arrptr( verts ), &coords2, 1, err ); 00226 if( err ) return 1; 00227 00228 PlanarDomain un_domain2( norm, coords2 ); 00229 00230 MeshDomainAssoc mesh_and_domain6 = MeshDomainAssoc( &mesh, &un_domain2 ); 00231 t.reset(); 00232 un_wrapper.run_instructions( &mesh_and_domain6, err ); 00233 if( err ) return 1; 00234 00235 double unb_l1_secs = t.since_birth(); 00236 std::cout << std::endl 00237 << "UntangleWrapper::BETA large file (term crit=0.005) completed in " << unb_l1_secs << " seconds" 00238 << std::endl; 00239 00240 mesh.clear(); 00241 mesh.read_vtk( untangle_file_name_2.c_str(), err ); 00242 if( err ) return 1; 00243 00244 // get domain 00245 verts.clear(); 00246 mesh.get_all_vertices( verts, err ); 00247 if( err || verts.empty() ) return 1; 00248 MsqVertex coords3; 00249 mesh.vertices_get_coordinates( arrptr( verts ), &coords3, 1, err ); 00250 if( err ) return 1; 00251 00252 PlanarDomain un_domain3( norm, coords3 ); 00253 00254 un_wrapper.set_vertex_movement_limit_factor( 0.1 ); 00255 MeshDomainAssoc mesh_and_domain7 = MeshDomainAssoc( &mesh, &un_domain3 ); 00256 t.reset(); 00257 un_wrapper.run_instructions( &mesh_and_domain7, err ); 00258 if( err ) return 1; 00259 00260 double unb_l2_secs = t.since_birth(); 00261 std::cout << std::endl 00262 << "UntangleWrapper::BETA large file (term crit=0.1) completed in " << unb_l2_secs << " seconds" 00263 << std::endl; 00264 00265 // #################### Begin UntangleWrapper::SIZE tests ################### 00266 00267 mesh.clear(); 00268 mesh.read_vtk( untangle_file_name_1.c_str(), err ); 00269 if( err ) return 1; 00270 00271 verts.clear(); 00272 mesh.get_all_vertices( verts, err ); 00273 if( err || verts.empty() ) return 1; 00274 MsqVertex coords2a; 00275 mesh.vertices_get_coordinates( arrptr( verts ), &coords2a, 1, err ); 00276 if( err ) return 1; 00277 PlanarDomain u_domain3( norm, coords2a ); 00278 00279 UntangleWrapper::UntangleMetric metric2 = UntangleWrapper::SIZE; 00280 UntangleWrapper un_wrapper2s( metric2 ); 00281 UntangleWrapper un_wrapper2l( metric2 ); 00282 00283 MeshDomainAssoc mesh_and_domain8 = MeshDomainAssoc( &mesh, &u_domain3 ); 00284 t.reset(); 00285 un_wrapper2s.run_instructions( &mesh_and_domain8, err ); 00286 if( err ) return 1; 00287 double uns_s_secs = t.since_birth(); 00288 std::cout << std::endl 00289 << "UntangleWrapper::SIZE small file optimization completed in " << uns_s_secs << " seconds" << std::endl; 00290 00291 mesh.clear(); 00292 mesh.read_vtk( untangle_file_name_2.c_str(), err ); 00293 if( err ) return 1; 00294 00295 // get domain 00296 verts.clear(); 00297 mesh.get_all_vertices( verts, err ); 00298 if( err || verts.empty() ) return 1; 00299 MsqVertex coords4; 00300 mesh.vertices_get_coordinates( arrptr( verts ), &coords4, 1, err ); 00301 if( err ) return 1; 00302 00303 PlanarDomain un_domain4( norm, coords4 ); 00304 un_wrapper2s.set_vertex_movement_limit_factor( 0.005 ); 00305 00306 MeshDomainAssoc mesh_and_domain9 = MeshDomainAssoc( &mesh, &un_domain4 ); 00307 t.reset(); 00308 un_wrapper2s.run_instructions( &mesh_and_domain9, err ); 00309 if( err ) return 1; 00310 00311 double uns_l1_secs = t.since_birth(); 00312 std::cout << std::endl 00313 << "UntangleWrappe::SIZE large file (term crit=0.005) completed in " << uns_l1_secs << " seconds" 00314 << std::endl; 00315 00316 mesh.clear(); 00317 mesh.read_vtk( untangle_file_name_2.c_str(), err ); 00318 if( err ) return 1; 00319 00320 mesh.get_all_vertices( verts, err ); 00321 if( err || verts.empty() ) return 1; 00322 00323 mesh.vertices_get_coordinates( arrptr( verts ), &coords4, 1, err ); 00324 if( err ) return 1; 00325 00326 un_wrapper2l.set_vertex_movement_limit_factor( 0.1 ); 00327 t.reset(); 00328 un_wrapper2l.run_instructions( &mesh_and_domain9, err ); 00329 if( err ) return 1; 00330 00331 double uns_l2_secs = t.since_birth(); 00332 std::cout << std::endl 00333 << "UntangleWrappe::SIZE large file (term crit=0.1) completed in " << uns_l2_secs << " seconds" 00334 << std::endl; 00335 00336 // #################### Begin UntangleWrapper::SHAPESIZE tests ################### 00337 00338 mesh.clear(); 00339 mesh.read_vtk( untangle_file_name_1.c_str(), err ); 00340 if( err ) return 1; 00341 00342 verts.clear(); 00343 mesh.get_all_vertices( verts, err ); 00344 if( err || verts.empty() ) return 1; 00345 MsqVertex coords5; 00346 mesh.vertices_get_coordinates( arrptr( verts ), &coords5, 1, err ); 00347 if( err ) return 1; 00348 PlanarDomain u_domain5( norm, coords3 ); 00349 00350 UntangleWrapper::UntangleMetric metric3 = UntangleWrapper::SHAPESIZE; 00351 UntangleWrapper un_wrapper3( metric3 ); 00352 00353 MeshDomainAssoc mesh_and_domain10 = MeshDomainAssoc( &mesh, &u_domain5 ); 00354 t.reset(); 00355 un_wrapper3.run_instructions( &mesh_and_domain10, err ); 00356 if( err ) return 1; 00357 00358 double unss_s_secs = t.since_birth(); 00359 std::cout << std::endl 00360 << "UntangleWrapper::SHAPESIZE small file optimization completed in " << unss_s_secs << " seconds" 00361 << std::endl; 00362 00363 mesh.clear(); 00364 mesh.read_vtk( untangle_file_name_2.c_str(), err ); 00365 if( err ) return 1; 00366 00367 // get domain 00368 verts.clear(); 00369 mesh.get_all_vertices( verts, err ); 00370 if( err || verts.empty() ) return 1; 00371 MsqVertex coords6; 00372 mesh.vertices_get_coordinates( arrptr( verts ), &coords6, 1, err ); 00373 if( err ) return 1; 00374 00375 PlanarDomain un_domain6( norm, coords6 ); 00376 00377 MeshDomainAssoc mesh_and_domain11 = MeshDomainAssoc( &mesh, &un_domain6 ); 00378 t.reset(); 00379 un_wrapper3.run_instructions( &mesh_and_domain11, err ); 00380 if( err ) return 1; 00381 00382 double unss_l_secs = t.since_birth(); 00383 std::cout << std::endl 00384 << "UntangleWrapper::SHAPESIZE large file optimization completed in " << unss_l_secs << " seconds" 00385 << std::endl; 00386 00387 // #################### Begin SizeAdaptShapeWrapper tests ################### 00388 00389 mesh.clear(); 00390 mesh.read_vtk( size_adapt_shape_file_name_1.c_str(), err ); 00391 if( err ) return 1; 00392 00393 elem_vec_t polar, equatorial; 00394 find_z10_extreme_elements( mesh, polar, equatorial, err ); 00395 if( err ) return 1; 00396 00397 double eq_min, eq_max, eq_mean, pol_min, pol_max, pol_mean; 00398 elem_areas( mesh, polar, pol_min, pol_mean, pol_max, err ); 00399 if( err ) return 1; 00400 elem_areas( mesh, equatorial, eq_min, eq_mean, eq_max, err ); 00401 if( err ) return 1; 00402 00403 SphericalDomain geom( Vector3D( 0, 0, 0 ), 10.0 ); 00404 SizeAdaptShapeWrapper sas_wrapper1( 1e-2, 50 ); 00405 SizeAdaptShapeWrapper sas_wrapper2( 1e-1, 50 ); 00406 00407 MeshDomainAssoc mesh_and_domain12 = MeshDomainAssoc( &mesh, &geom ); 00408 t.reset(); 00409 sas_wrapper1.run_instructions( &mesh_and_domain12, err ); 00410 if( err ) return 1; 00411 double sas1_secs = t.since_birth(); 00412 std::cout << std::endl 00413 << "SizeAdaptShapeWrapper (term crit=0.01) completed in " << sas1_secs << " seconds" << std::endl; 00414 00415 mesh.clear(); 00416 mesh.read_vtk( size_adapt_shape_file_name_1.c_str(), err ); 00417 if( err ) return 1; 00418 00419 t.reset(); 00420 sas_wrapper2.run_instructions( &mesh_and_domain12, err ); 00421 if( err ) return 1; 00422 double sas2_secs = t.since_birth(); 00423 std::cout << std::endl 00424 << "SizeAdaptShapeWrapper (term crit=0.1) completed in " << sas2_secs << " seconds" << std::endl; 00425 00426 // #################### Begin PaverMinEdgeLengthWrapper tests ################### 00427 00428 PaverMinEdgeLengthWrapper mel_wrapper1( .005, 50 ); 00429 PaverMinEdgeLengthWrapper mel_wrapper2( .1, 50 ); 00430 00431 mesh.clear(); 00432 mesh.read_vtk( min_edge_length_file_name_1.c_str(), err ); 00433 00434 t.reset(); 00435 mel_wrapper1.run_instructions( &mesh, err ); 00436 if( err ) return 1; 00437 double mel_s_secs = t.since_birth(); 00438 std::cout << std::endl 00439 << "PaverMinEdgeLengthWrapper small file optimization completed in " << mel_s_secs << " seconds" 00440 << std::endl; 00441 00442 mesh.clear(); 00443 mesh.read_vtk( min_edge_length_file_name_2.c_str(), err ); 00444 00445 t.reset(); 00446 mel_wrapper1.run_instructions( &mesh, err ); 00447 if( err ) return 1; 00448 double mel1_l_secs = t.since_birth(); 00449 std::cout << std::endl 00450 << "PaverMinEdgeLengthWrapper large file (term crit=0.005) completed in " << mel1_l_secs << " seconds" 00451 << std::endl; 00452 00453 mesh.clear(); 00454 mesh.read_vtk( min_edge_length_file_name_2.c_str(), err ); 00455 t.reset(); 00456 mel_wrapper2.run_instructions( &mesh, err ); 00457 if( err ) return 1; 00458 double mel2_l_secs = t.since_birth(); 00459 std::cout << std::endl 00460 << "PaverMinEdgeLengthWrapper large file (term crit=0.1) completed in " << mel2_l_secs << " seconds" 00461 << std::endl; 00462 00463 // #################### Begin DeformingDomainWrapper tests ################### 00464 00465 // load mesh 00466 mesh.clear(); 00467 mesh.read_vtk( deforming_domain_file_name_1.c_str(), err ); 00468 if( MSQ_CHKERR( err ) ) return 1; 00469 00470 std::vector< Mesh::VertexHandle > curves[4]; 00471 Mesh::VertexHandle corners[4]; 00472 classify_boundary( &mesh, corners, curves, err ); 00473 if( MSQ_CHKERR( err ) ) return 1; 00474 00475 // new, "deformed" domain will be an 2HDx2HD planar square 00476 const double corner_coords[][3] = { { -HD, -HD, Z }, { HD, -HD, Z }, { HD, HD, Z }, { -HD, HD, Z } }; 00477 LineDomain lines[4] = { LineDomain( Vector3D( corner_coords[0] ), Vector3D( 1, 0, 0 ) ), 00478 LineDomain( Vector3D( corner_coords[1] ), Vector3D( 0, 1, 0 ) ), 00479 LineDomain( Vector3D( corner_coords[2] ), Vector3D( -1, 0, 0 ) ), 00480 LineDomain( Vector3D( corner_coords[3] ), Vector3D( 0, -1, 0 ) ) }; 00481 PlanarDomain surface( PlanarDomain::XY, Z ); 00482 00483 // save initial mesh state 00484 DeformingCurveSmoother curve_tool; 00485 for( int i = 0; i < 4; ++i ) 00486 { 00487 curve_tool.store_initial_mesh( &mesh, &curves[i][0], curves[i].size(), &lines[i], err ); 00488 if( MSQ_CHKERR( err ) ) return 1; 00489 } 00490 DeformingDomainWrapper dd_wrapper; 00491 dd_wrapper.store_initial_mesh( &mesh, err ); 00492 if( MSQ_CHKERR( err ) ) return 1; 00493 00494 // move corner vertices to new location 00495 for( int i = 0; i < 4; ++i ) 00496 { 00497 Vector3D vect( corner_coords[i] ); 00498 mesh.vertex_set_coordinates( corners[i], vect, err ); 00499 if( MSQ_CHKERR( err ) ) return 1; 00500 } 00501 std::vector< bool > fixed( 4, true ); 00502 mesh.vertices_set_fixed_flag( corners, fixed, 4, err ); 00503 if( MSQ_CHKERR( err ) ) return 1; 00504 00505 // smooth curves 00506 for( int i = 0; i < 4; ++i ) 00507 { 00508 curve_tool.smooth_curve( &mesh, &curves[i][0], curves[i].size(), &lines[i], 00509 DeformingCurveSmoother::PROPORTIONAL, err ); 00510 if( MSQ_CHKERR( err ) ) return 1; 00511 fixed.resize( curves[i].size(), true ); 00512 mesh.vertices_set_fixed_flag( &curves[i][0], fixed, curves[i].size(), err ); 00513 if( MSQ_CHKERR( err ) ) return 1; 00514 } 00515 00516 MeshDomainAssoc mesh_and_domain1 = MeshDomainAssoc( &mesh, &surface ); 00517 t.reset(); 00518 dd_wrapper.run_instructions( &mesh_and_domain1, err ); 00519 if( MSQ_CHKERR( err ) ) return 1; 00520 double dd_secs = t.since_birth(); 00521 std::cout << std::endl 00522 << "DeformingDomainWrapper file (term crit=0.01) completed in " << dd_secs << " seconds" << std::endl; 00523 00524 // Do it all again for the next test 00525 mesh.clear(); 00526 mesh.read_vtk( deforming_domain_file_name_1.c_str(), err ); 00527 if( MSQ_CHKERR( err ) ) return 1; 00528 00529 std::vector< Mesh::VertexHandle > curves2[4]; 00530 Mesh::VertexHandle corners2[4]; 00531 classify_boundary( &mesh, corners2, curves2, err ); 00532 if( MSQ_CHKERR( err ) ) return 1; 00533 00534 // new, "deformed" domain will be an 2HDx2HD planar square 00535 const double corner_coords2[][3] = { { -HD, -HD, Z }, { HD, -HD, Z }, { HD, HD, Z }, { -HD, HD, Z } }; 00536 LineDomain lines2[4] = { LineDomain( Vector3D( corner_coords2[0] ), Vector3D( 1, 0, 0 ) ), 00537 LineDomain( Vector3D( corner_coords2[1] ), Vector3D( 0, 1, 0 ) ), 00538 LineDomain( Vector3D( corner_coords2[2] ), Vector3D( -1, 0, 0 ) ), 00539 LineDomain( Vector3D( corner_coords2[3] ), Vector3D( 0, -1, 0 ) ) }; 00540 PlanarDomain surface2( PlanarDomain::XY, Z ); 00541 00542 // save initial mesh state 00543 DeformingCurveSmoother curve_tool2; 00544 for( int i = 0; i < 4; ++i ) 00545 { 00546 curve_tool2.store_initial_mesh( &mesh, &curves2[i][0], curves2[i].size(), &lines2[i], err ); 00547 if( MSQ_CHKERR( err ) ) return 1; 00548 } 00549 DeformingDomainWrapper dd_wrapper2; 00550 dd_wrapper2.store_initial_mesh( &mesh, err ); 00551 if( MSQ_CHKERR( err ) ) return 1; 00552 00553 // move corner vertices to new location 00554 for( int i = 0; i < 4; ++i ) 00555 { 00556 Vector3D vect( corner_coords2[i] ); 00557 mesh.vertex_set_coordinates( corners2[i], vect, err ); 00558 if( MSQ_CHKERR( err ) ) return 1; 00559 } 00560 std::vector< bool > fixed2( 4, true ); 00561 mesh.vertices_set_fixed_flag( corners2, fixed2, 4, err ); 00562 if( MSQ_CHKERR( err ) ) return 1; 00563 00564 // smooth curves 00565 for( int i = 0; i < 4; ++i ) 00566 { 00567 curve_tool2.smooth_curve( &mesh, &curves2[i][0], curves2[i].size(), &lines2[i], 00568 DeformingCurveSmoother::PROPORTIONAL, err ); 00569 if( MSQ_CHKERR( err ) ) return 1; 00570 fixed2.resize( curves2[i].size(), true ); 00571 mesh.vertices_set_fixed_flag( &curves2[i][0], fixed2, curves2[i].size(), err ); 00572 if( MSQ_CHKERR( err ) ) return 1; 00573 } 00574 00575 dd_wrapper2.set_vertex_movement_limit_factor( 0.1 ); 00576 MeshDomainAssoc mesh_and_domain2 = MeshDomainAssoc( &mesh, &surface2 ); 00577 t.reset(); 00578 dd_wrapper2.run_instructions( &mesh_and_domain2, err ); 00579 if( MSQ_CHKERR( err ) ) return 1; 00580 double dd_secs2 = t.since_birth(); 00581 std::cout << std::endl 00582 << "DeformingDomainWrapper file (term crit=0.1) completed in " << dd_secs2 << " seconds" << std::endl; 00583 00584 // Timing Summary 00585 std::cout << std::endl << "********* Wrappers Timing Summary **********" << std::endl << std::endl; 00586 std::cout << "ShapeImprover small file optimization completed in " << si_s_secs << " seconds" << std::endl; 00587 std::cout << "ShapeImprover large file optimization completed in " << si_l_secs << " seconds" << std::endl; 00588 std::cout << "LaplacianWrapper small file optimization completed in " << lp_s_secs << " seconds" << std::endl; 00589 std::cout << "LaplacianWrapper large file optimization (term crit=0.001) in " << lp_l1_secs << " seconds" 00590 << std::endl; 00591 std::cout << "LaplacianWrapper large file optimization (term crit=0.1) in " << lp_l2_secs << " seconds" 00592 << std::endl; 00593 std::cout << "UntangleWrapper::BETA small file optimization completed in " << unb_s_secs << " seconds" << std::endl; 00594 std::cout << "UntangleWrapper::BETA large file (term crit=0.005) completed in " << unb_l1_secs << " seconds" 00595 << std::endl; 00596 std::cout << "UntangleWrapper::BETA large file (term crit=0.1) completed in " << unb_l2_secs << " seconds" 00597 << std::endl; 00598 std::cout << "UntangleWrapper::SIZE small file optimization completed in " << uns_s_secs << " seconds" << std::endl; 00599 std::cout << "UntangleWrapper::SIZE large file (term crit=0.005) completed in " << uns_l1_secs << " seconds" 00600 << std::endl; 00601 std::cout << "UntangleWrapper::SIZE large file (term crit=0.1) completed in " << uns_l2_secs << " seconds" 00602 << std::endl; 00603 std::cout << "UntangleWrapper::SHAPESIZE small file optimization completed in " << unss_s_secs << " seconds" 00604 << std::endl; 00605 std::cout << "UntangleWrapper::SHAPESIZE large file optimization completed in " << unss_l_secs << " seconds" 00606 << std::endl; 00607 std::cout << "SizeAdaptShapeWrapper (term crit=0.01) completed in " << sas1_secs << " seconds" << std::endl; 00608 std::cout << "SizeAdaptShapeWrapper (term crit=0.1) completed in " << sas2_secs << " seconds" << std::endl; 00609 std::cout << "PaverMinEdgeLengthWrapper small file optimization completed in " << mel_s_secs << " seconds" 00610 << std::endl; 00611 std::cout << "PaverMinEdgeLengthWrapper large file (term crit=0.005) completed in " << mel1_l_secs << " seconds" 00612 << std::endl; 00613 std::cout << "PaverMinEdgeLengthWrapper large file (term crit=0.1) completed in " << mel2_l_secs << " seconds" 00614 << std::endl; 00615 std::cout << "DeformingDomainWrapper file (term crit=0.01) completed in " << dd_secs << " seconds" << std::endl; 00616 std::cout << "DeformingDomainWrapper file (term crit=0.1) completed in " << dd_secs2 << " seconds" << std::endl; 00617 00618 return 0; 00619 } 00620 00621 void find_z10_extreme_elements( Mesh& mesh, elem_vec_t& polar_elems, elem_vec_t& equatorial_elems, MsqError& err ) 00622 { 00623 elem_vec_t elems; 00624 mesh.get_all_elements( elems, err );MSQ_ERRRTN( err ); 00625 00626 std::vector< Mesh::VertexHandle > verts; 00627 std::vector< MsqVertex > coords; 00628 std::vector< size_t > junk; 00629 for( elem_vec_t::iterator i = elems.begin(); i != elems.end(); ++i ) 00630 { 00631 verts.clear(); 00632 junk.clear(); 00633 mesh.elements_get_attached_vertices( &*i, 1, verts, junk, err );MSQ_ERRRTN( err ); 00634 coords.resize( verts.size() ); 00635 mesh.vertices_get_coordinates( arrptr( verts ), arrptr( coords ), verts.size(), err );MSQ_ERRRTN( err ); 00636 00637 for( std::vector< MsqVertex >::iterator j = coords.begin(); j != coords.end(); ++j ) 00638 { 00639 double z = ( *j )[2]; 00640 if( fabs( z ) < 1e-6 ) 00641 { 00642 equatorial_elems.push_back( *i ); 00643 break; 00644 } 00645 else if( fabs( z ) - 10 < 1e-6 ) 00646 { 00647 polar_elems.push_back( *i ); 00648 break; 00649 } 00650 } 00651 } 00652 } 00653 00654 void elem_areas( Mesh& mesh, const elem_vec_t& elems, double& min, double& mean, double& max, MsqError& err ) 00655 { 00656 min = HUGE_VAL; 00657 max = -1; 00658 mean = 0.0; 00659 00660 std::vector< EntityTopology > types( elems.size() ); 00661 mesh.elements_get_topologies( arrptr( elems ), arrptr( types ), elems.size(), err );MSQ_ERRRTN( err ); 00662 00663 std::vector< Mesh::VertexHandle > verts; 00664 std::vector< MsqVertex > coords; 00665 std::vector< size_t > junk; 00666 for( size_t i = 0; i < elems.size(); ++i ) 00667 { 00668 verts.clear(); 00669 junk.clear(); 00670 mesh.elements_get_attached_vertices( &elems[i], 1, verts, junk, err );MSQ_ERRRTN( err ); 00671 coords.resize( verts.size() ); 00672 mesh.vertices_get_coordinates( arrptr( verts ), arrptr( coords ), verts.size(), err );MSQ_ERRRTN( err ); 00673 00674 Vector3D v1, v2; 00675 double area; 00676 if( types[i] == TRIANGLE ) 00677 { 00678 assert( coords.size() == 3 ); 00679 v1 = coords[1] - coords[0]; 00680 v2 = coords[2] - coords[0]; 00681 area = 0.5 * ( v1 * v2 ).length(); 00682 } 00683 else if( types[i] == QUADRILATERAL ) 00684 { 00685 assert( coords.size() == 4 ); 00686 v1 = coords[0] + coords[1] - coords[2] - coords[3]; 00687 v2 = coords[0] + coords[3] - coords[1] - coords[2]; 00688 area = 0.25 * ( v1 * v2 ).length(); 00689 } 00690 else 00691 { 00692 MSQ_SETERR( err ) 00693 ( "Input file contains volume elements", MsqError::UNSUPPORTED_ELEMENT ); 00694 return; 00695 } 00696 00697 if( min > area ) min = area; 00698 if( max < area ) max = area; 00699 mean += area; 00700 } 00701 00702 mean /= elems.size(); 00703 } 00704 00705 #define CHKMESH( COND, ERR ) \ 00706 do \ 00707 { \ 00708 if( !( COND ) ) \ 00709 { \ 00710 MSQ_SETERR( err )( "Unexpected mesh topology", MsqError::INVALID_MESH ); \ 00711 return; \ 00712 } \ 00713 } while( false ) 00714 00715 void classify_boundary( Mesh* mesh, 00716 Mesh::VertexHandle corners_out[4], 00717 std::vector< Mesh::VertexHandle > curves_out[4], 00718 MsqError& err ) 00719 { 00720 std::vector< Mesh::VertexHandle > verts; 00721 mesh->get_all_vertices( verts, err );MSQ_ERRRTN( err ); 00722 00723 // Find the corner vertex that has negative X and Y coordinates 00724 Mesh::VertexHandle start; 00725 bool have_start = false; 00726 std::vector< Mesh::ElementHandle > elems; 00727 std::vector< size_t > offsets; 00728 for( size_t i = 0; i < verts.size(); ++i ) 00729 { 00730 elems.clear(); 00731 offsets.clear(); 00732 mesh->vertices_get_attached_elements( &verts[i], 1, elems, offsets, err );MSQ_ERRRTN( err ); 00733 if( elems.size() == 1 ) 00734 { 00735 MsqVertex coords; 00736 mesh->vertices_get_coordinates( &verts[i], &coords, 1, err );MSQ_ERRRTN( err ); 00737 if( coords[0] < 0.0 && coords[1] < 0.0 ) 00738 { 00739 CHKMESH( !have_start, err ); 00740 have_start = true; 00741 start = verts[i]; 00742 } 00743 } 00744 } 00745 CHKMESH( have_start, err ); 00746 00747 // starting at a a corner vertex, find skin vertices 00748 std::vector< Mesh::VertexHandle > boundary; 00749 boundary.push_back( start ); 00750 elems.clear(); 00751 offsets.clear(); 00752 mesh->vertices_get_attached_elements( &start, 1, elems, offsets, err );MSQ_ERRRTN( err ); 00753 Mesh::ElementHandle prev = elems.front(); 00754 corners_out[0] = start; 00755 int ncorner = 1; 00756 do 00757 { 00758 verts.clear(); 00759 offsets.clear(); 00760 mesh->elements_get_attached_vertices( &prev, 1, verts, offsets, err ); 00761 size_t idx = std::find( verts.begin(), verts.end(), boundary.back() ) - verts.begin(); 00762 CHKMESH( idx < verts.size(), err ); 00763 00764 Mesh::VertexHandle next = verts[( idx + 1 ) % verts.size()]; 00765 elems.clear(); 00766 offsets.clear(); 00767 mesh->vertices_get_attached_elements( &next, 1, elems, offsets, err );MSQ_ERRRTN( err ); 00768 CHKMESH( elems.size() == 1 || elems.size() == 2, err ); 00769 if( elems.size() == 2 ) 00770 { 00771 idx = std::find( elems.begin(), elems.end(), prev ) - elems.begin(); 00772 CHKMESH( idx < 2, err ); 00773 prev = elems[1 - idx]; 00774 } 00775 00776 if( elems.size() == 1 ) 00777 { 00778 CHKMESH( ncorner < 5, err ); 00779 if( ncorner < 4 ) // don't add start vertx twice 00780 corners_out[ncorner] = next; 00781 ++ncorner; 00782 } 00783 00784 boundary.push_back( next ); 00785 } while( boundary.front() != boundary.back() ); 00786 CHKMESH( ncorner == 5, err ); 00787 00788 // find curve vertices 00789 size_t idx = 0; 00790 for( int c = 0; c < 4; ++c ) 00791 { 00792 Mesh::VertexHandle s, e; 00793 s = corners_out[c]; 00794 e = corners_out[( c + 1 ) % 4]; 00795 CHKMESH( idx < boundary.size(), err ); 00796 CHKMESH( boundary[idx] == s, err ); 00797 for( ; boundary[idx] != e; ++idx ) 00798 curves_out[c].push_back( boundary[idx] ); 00799 curves_out[c].push_back( e ); 00800 } 00801 }
1.7.6.1 
