MOAB: Mesh Oriented datABase  (version 5.3.1)
pyramid_test.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     diachin2@llnl.gov, djmelan@sandia.gov, mbrewer@sandia.gov,
00024     pknupp@sandia.gov, tleurent@mcs.anl.gov, tmunson@mcs.anl.gov
00025     kraftche@cae.wisc.edu
00026 
00027   ***************************************************************** */
00028 
00029 #define TOL 1e-5
00030 
00031 #include <iostream>
00032 using std::cout;
00033 using std::endl;
00034 #include <cstdlib>
00035 
00036 #include "Mesquite.hpp"
00037 #include "MeshImpl.hpp"
00038 #include "MsqError.hpp"
00039 #include "InstructionQueue.hpp"
00040 #include "TerminationCriterion.hpp"
00041 #include "QualityAssessor.hpp"
00042 #include "LPtoPTemplate.hpp"
00043 #include "LInfTemplate.hpp"
00044 
00045 // algorithms
00046 #include "IdealWeightMeanRatio.hpp"
00047 #include "IdealWeightInverseMeanRatio.hpp"
00048 #include "UntangleBetaQualityMetric.hpp"
00049 #include "FeasibleNewton.hpp"
00050 #include "ConjugateGradient.hpp"
00051 #include "TrustRegion.hpp"
00052 #include "ConditionNumberQualityMetric.hpp"
00053 #include "TestUtil.hpp"
00054 
00055 using namespace MBMesquite;
00056 
00057 // Use CPPUNIT_ASSERT in code so it's easy to convert to a unit test later.
00058 #define CPPUNIT_ASSERT( A )                                       \
00059     do                                                            \
00060     {                                                             \
00061         if( !( A ) )                                              \
00062         {                                                         \
00063             std::cout << "Assertion Failed: " << #A << std::endl; \
00064             std::cout << "  File: " << __FILE__ << std::endl;     \
00065             std::cout << "  Line: " << __LINE__ << std::endl;     \
00066             return 1;                                             \
00067         }                                                         \
00068     } while( false )
00069 
00070 // Given a mesh with a single free vertex located at the origin,
00071 // move the vertex to the specified position, smooth the mesh,
00072 // and verify that the vertex was moved back to the origin by
00073 // the smoother.
00074 bool smooth_mesh( Mesh* mesh,
00075                   Mesh* ref_mesh,
00076                   Mesh::VertexHandle free_vertex_at_origin,
00077                   Vector3D initial_free_vertex_position,
00078                   QualityMetric* metric );
00079 
00080 bool smooth_mixed_mesh( const char* filename );
00081 
00082 int main( int argc, char* argv[] )
00083 {
00084     unsigned i;
00085     std::string input_file = TestDir + "unittest/mesquite/3D/vtk/mixed/tangled/mixed-hex-pyr-tet.vtk";
00086     if( argc == 2 )
00087         input_file = argv[1];
00088     else if( argc != 1 )
00089     {
00090         std::cerr << "Invalid arguments.\n";
00091         return 2;
00092     }
00093 
00094     MBMesquite::MsqPrintError err( cout );
00095     IdealWeightMeanRatio m1;
00096     IdealWeightInverseMeanRatio m2( err );
00097     ConditionNumberQualityMetric m3;
00098     QualityMetric* metrics[] = { &m1, &m2, &m3, 0 };
00099 
00100     // Read Mesh
00101     std::string mesh_file  = TestDir + "unittest/mesquite/3D/vtk/pyramids/untangled/12-pyramid-unit-sphere.vtk";
00102     std::string imesh_file = TestDir + "unittest/mesquite/3D/vtk/pyramids/untangled/12-pyramid-unit-sphere.vtk";
00103     MBMesquite::MeshImpl mesh;
00104     mesh.read_vtk( mesh_file.c_str(), err );
00105     CPPUNIT_ASSERT( !err );
00106     MBMesquite::MeshImpl ideal_mesh;
00107     ideal_mesh.read_vtk( imesh_file.c_str(), err );
00108     CPPUNIT_ASSERT( !err );
00109 
00110     // Check that the mesh read correctly, and contains what is
00111     // expected later.
00112 
00113     // Get mesh data
00114     // Expecting file to contain 12 pyramid elements constructed
00115     // from 15 vertices.
00116     std::vector< Mesh::VertexHandle > vert_array;
00117     std::vector< Mesh::ElementHandle > elem_array;
00118     std::vector< size_t > conn_offsets;
00119     mesh.get_all_elements( elem_array, err );
00120     CPPUNIT_ASSERT( !err );
00121     CPPUNIT_ASSERT( elem_array.size() == 12 );
00122     mesh.elements_get_attached_vertices( arrptr( elem_array ), elem_array.size(), vert_array, conn_offsets, err );
00123     CPPUNIT_ASSERT( !err );
00124     CPPUNIT_ASSERT( vert_array.size() == 60 );
00125     CPPUNIT_ASSERT( conn_offsets.size() == 13 );
00126     EntityTopology type_array[12];
00127     mesh.elements_get_topologies( arrptr( elem_array ), type_array, 12, err );
00128     CPPUNIT_ASSERT( !err );
00129 
00130     // Verify element types and number of vertices
00131     for( i = 0; i < 12; ++i )
00132     {
00133         CPPUNIT_ASSERT( type_array[i] == PYRAMID );
00134         CPPUNIT_ASSERT( conn_offsets[i] == 5 * i );
00135     }
00136 
00137     // All pyramids should share a common apex, at the
00138     // center of the sphere
00139     Mesh::VertexHandle apex_handle = vert_array[4];
00140     for( i = 1; i < 12; ++i )
00141     {
00142         CPPUNIT_ASSERT( vert_array[5 * i + 4] == apex_handle );
00143     }
00144 
00145     // Verify that apex is at origin and all other vertices are
00146     // on unit sphere
00147     MsqVertex vertices[60];
00148     mesh.vertices_get_coordinates( arrptr( vert_array ), vertices, 60, err );
00149     CPPUNIT_ASSERT( !err );
00150     for( i = 0; i < 60; ++i )
00151     {
00152         if( vert_array[i] == apex_handle )
00153             CPPUNIT_ASSERT( vertices[i].within_tolerance_box( Vector3D( 0, 0, 0 ), 1e-6 ) );
00154         else
00155             CPPUNIT_ASSERT( fabs( 1.0 - vertices[i].length() ) < 1e-6 );
00156     }
00157 
00158     // Try smoothing w/out moving the free vertex and verify that
00159     // the smoother didn't move the vertex
00160     Vector3D position( 0, 0, 0 );
00161     for( i = 0; metrics[i] != NULL; ++i )
00162         CPPUNIT_ASSERT( !smooth_mesh( &mesh, &ideal_mesh, apex_handle, position, metrics[i] ) );
00163 
00164     // Now try moving the vertex and see if the smoother moves it back
00165     // to the origin
00166     position.set( 0.1, 0.1, 0.1 );
00167     for( i = 0; metrics[i] != NULL; ++i )
00168         CPPUNIT_ASSERT( !smooth_mesh( &mesh, &ideal_mesh, apex_handle, position, metrics[i] ) );
00169 
00170     // Now try moving the vertex further and see if the smoother moves it back
00171     // to the origin
00172     position.set( 0.3, 0.3, 0.3 );
00173     for( i = 0; metrics[i] != NULL; ++i )
00174         CPPUNIT_ASSERT( !smooth_mesh( &mesh, &ideal_mesh, apex_handle, position, metrics[i] ) );
00175 
00176     // Now try smoothing a real mixed mesh
00177     CPPUNIT_ASSERT( !smooth_mixed_mesh( input_file.c_str() ) );
00178 
00179     return 0;
00180 }
00181 
00182 bool smooth_mesh( Mesh* mesh,
00183                   Mesh*,
00184                   Mesh::VertexHandle free_vertex_at_origin,
00185                   Vector3D initial_free_vertex_position,
00186                   QualityMetric* metric )
00187 {
00188     MBMesquite::MsqPrintError err( cout );
00189     const Vector3D origin( 0, 0, 0 );
00190 
00191     // print a little output so we know when we died
00192     std::cout << "**************************************************************************" << std::endl
00193               << "* Smoothing..." << std::endl
00194               << "* Metric: " << metric->get_name() << std::endl
00195               << "* Apex position: " << initial_free_vertex_position
00196               << std::endl  //<<
00197               //"**************************************************************************"
00198               << std::endl;
00199 
00200     // Set free vertex to specified position
00201     mesh->vertex_set_coordinates( free_vertex_at_origin, initial_free_vertex_position, err );
00202     CPPUNIT_ASSERT( !err );
00203 
00204     // Create an InstructionQueue
00205     InstructionQueue Q;
00206 
00207     // Set up objective function
00208     LPtoPTemplate obj_func( metric, 1, err );
00209     CPPUNIT_ASSERT( !err );
00210 
00211     // Create solver
00212     FeasibleNewton solver( &obj_func );
00213     CPPUNIT_ASSERT( !err );
00214     solver.use_global_patch();
00215     CPPUNIT_ASSERT( !err );
00216 
00217     // Set stoping criteria for solver
00218     TerminationCriterion tc_inner;
00219     tc_inner.add_absolute_vertex_movement( 1e-6 );
00220     solver.set_inner_termination_criterion( &tc_inner );
00221 
00222     TerminationCriterion tc_outer;
00223     tc_outer.add_iteration_limit( 1 );
00224     solver.set_outer_termination_criterion( &tc_outer );
00225 
00226     // Add solver to queue
00227     Q.set_master_quality_improver( &solver, err );
00228     CPPUNIT_ASSERT( !err );
00229 
00230     // And smooth...
00231     Q.run_instructions( mesh, err );
00232     CPPUNIT_ASSERT( !err );
00233 
00234     // Verify that vertex was moved back to origin
00235     MsqVertex vtx;
00236     mesh->vertices_get_coordinates( &free_vertex_at_origin, &vtx, 1, err );
00237     CPPUNIT_ASSERT( !err );
00238     Vector3D position = vtx;
00239 
00240     // print a little output so we know when we died
00241     std::cout  //<<
00242                //"**************************************************************************"
00243         << std::endl
00244         << "* Done Smoothing:" << std::endl
00245         << "* Metric: " << metric->get_name() << std::endl
00246         << "* Apex position: " << position << std::endl
00247         << "**************************************************************************" << std::endl;
00248 
00249     CPPUNIT_ASSERT( position.within_tolerance_box( Vector3D( 0, 0, 0 ), TOL ) );
00250     return false;
00251 }
00252 
00253 bool smooth_mixed_mesh( const char* filename )
00254 {
00255     MBMesquite::MsqPrintError err( cout );
00256 
00257     // print a little output so we know when we died
00258     std::cout << "**************************************************************************" << std::endl
00259               << "* Smoothing: " << filename << std::endl
00260               << "**************************************************************************" << std::endl;
00261 
00262     // The instruction queue to set up
00263     InstructionQueue Q;
00264 
00265     // Use numeric approx of derivitives until analytic solutions
00266     // are working for pyramids
00267     IdealWeightInverseMeanRatio mr_metric( err );
00268     // sRI_DFT dft_metric;
00269     UntangleBetaQualityMetric un_metric( 0 );
00270     CPPUNIT_ASSERT( !err );
00271 
00272     // Create Mesh object
00273     MBMesquite::MeshImpl mesh;
00274     mesh.read_vtk( filename, err );
00275     CPPUNIT_ASSERT( !err );
00276 
00277     // Set up a preconditioner
00278     LInfTemplate pre_obj_func( &un_metric );
00279     ConjugateGradient precond( &pre_obj_func, err );
00280     CPPUNIT_ASSERT( !err );
00281     precond.use_element_on_vertex_patch();
00282     TerminationCriterion pre_term, pre_outer;
00283     // pre_term.add_relative_quality_improvement( 0.1 );
00284     pre_term.add_iteration_limit( 3 );
00285     pre_outer.add_iteration_limit( 1 );
00286     CPPUNIT_ASSERT( !err );
00287     precond.set_inner_termination_criterion( &pre_term );
00288     precond.set_outer_termination_criterion( &pre_outer );
00289     // precond.use_element_on_vertex_patch();
00290 
00291     // Set up objective function
00292     LPtoPTemplate obj_func( &mr_metric, 1, err );
00293     CPPUNIT_ASSERT( !err );
00294 
00295     // Create solver
00296     FeasibleNewton solver( &obj_func );
00297     CPPUNIT_ASSERT( !err );
00298     solver.use_global_patch();
00299     CPPUNIT_ASSERT( !err );
00300 
00301     // Set stoping criteria for solver
00302     TerminationCriterion tc_inner;
00303     tc_inner.add_relative_quality_improvement( 0.25 );
00304     solver.set_inner_termination_criterion( &tc_inner );
00305 
00306     TerminationCriterion tc_outer;
00307     tc_outer.add_iteration_limit( 1 );
00308     CPPUNIT_ASSERT( !err );
00309     solver.set_outer_termination_criterion( &tc_outer );
00310 
00311     // Create a QualityAssessor
00312     MBMesquite::QualityAssessor qa;
00313     qa.add_quality_assessment( &mr_metric );
00314     qa.add_quality_assessment( &un_metric );
00315     Q.add_quality_assessor( &qa, err );
00316     CPPUNIT_ASSERT( !err );
00317 
00318     // Add untangler to queue
00319     Q.add_preconditioner( &precond, err );
00320     CPPUNIT_ASSERT( !err );
00321     Q.add_quality_assessor( &qa, err );
00322     CPPUNIT_ASSERT( !err );
00323 
00324     // Add solver to queue
00325     Q.set_master_quality_improver( &solver, err );
00326     CPPUNIT_ASSERT( !err );
00327     Q.add_quality_assessor( &qa, err );
00328     CPPUNIT_ASSERT( !err );
00329 
00330     // And smooth...
00331     Q.run_instructions( &mesh, err );
00332     CPPUNIT_ASSERT( !err );
00333 
00334     return false;
00335 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines