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