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