MOAB: Mesh Oriented datABase  (version 5.3.0)
benchmark_tests.cpp
Go to the documentation of this file.
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 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines