MOAB: Mesh Oriented datABase  (version 5.3.1)
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 =
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 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines