MOAB: Mesh Oriented datABase  (version 5.3.1)
wedge_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 //#define DO_QUALITY_ASSESSOR
00032 
00033 #include <iostream>
00034 using std::cout;
00035 using std::endl;
00036 #include <cstdlib>
00037 
00038 #include "Mesquite.hpp"
00039 #include "MeshImpl.hpp"
00040 #include "MsqError.hpp"
00041 #include "InstructionQueue.hpp"
00042 #include "TerminationCriterion.hpp"
00043 #include "QualityAssessor.hpp"
00044 #include "LPtoPTemplate.hpp"
00045 #include "LInfTemplate.hpp"
00046 #include "TestUtil.hpp"
00047 
00048 // algorithms
00049 #include "IdealWeightMeanRatio.hpp"
00050 #include "IdealWeightInverseMeanRatio.hpp"
00051 #include "UntangleBetaQualityMetric.hpp"
00052 #include "FeasibleNewton.hpp"
00053 #include "ConditionNumberQualityMetric.hpp"
00054 using namespace MBMesquite;
00055 
00056 // Use CPPUNIT_ASSERT in code so it's easy to convert to a unit test later.
00057 #define CPPUNIT_ASSERT( A )                                       \
00058     do                                                            \
00059     {                                                             \
00060         if( !( A ) )                                              \
00061         {                                                         \
00062             std::cout << "Assertion Failed: " << #A << std::endl; \
00063             std::cout << "  File: " << __FILE__ << std::endl;     \
00064             std::cout << "  Line: " << __LINE__ << std::endl;     \
00065             return 1;                                             \
00066         }                                                         \
00067     } while( false )
00068 
00069 #define CPPUNIT_ASSERT_DOUBLES_EQUAL( E, V, T )                                    \
00070     do                                                                             \
00071     {                                                                              \
00072         if( fabs( ( E ) - ( V ) ) > ( T ) )                                        \
00073         {                                                                          \
00074             std::cout << "Assertion Failed: " << #V << " == " << #E << std::endl;  \
00075             std::cout << "Expected: " << ( E ) << "  Got: " << ( V ) << std::endl; \
00076             std::cout << "  File: " << __FILE__ << std::endl;                      \
00077             std::cout << "  Line: " << __LINE__ << std::endl;                      \
00078             return 1;                                                              \
00079         }                                                                          \
00080     } while( false )
00081 
00082 #define CPPUNIT_ASSERT_EQUAL( E, V )                                               \
00083     do                                                                             \
00084     {                                                                              \
00085         if( ( E ) != ( V ) )                                                       \
00086         {                                                                          \
00087             std::cout << "Assertion Failed: " << #V << " == " << #E << std::endl;  \
00088             std::cout << "Expected: " << ( E ) << "  Got: " << ( V ) << std::endl; \
00089             std::cout << "  File: " << __FILE__ << std::endl;                      \
00090             std::cout << "  Line: " << __LINE__ << std::endl;                      \
00091             return 1;                                                              \
00092         }                                                                          \
00093     } while( false )
00094 
00095 #define CPPUNIT_ASSERT_VECTORS_EQUAL( E, V, T )                                    \
00096     do                                                                             \
00097     {                                                                              \
00098         if( !( E ).within_tolerance_box( V, T ) )                                  \
00099         {                                                                          \
00100             std::cout << "Assertion Failed: " << #V << " == " << #E << std::endl;  \
00101             std::cout << "Expected: " << ( E ) << "  Got: " << ( V ) << std::endl; \
00102             std::cout << "  File: " << __FILE__ << std::endl;                      \
00103             std::cout << "  Line: " << __LINE__ << std::endl;                      \
00104             return 1;                                                              \
00105         }                                                                          \
00106     } while( false )
00107 
00108 // Given a mesh with a single free vertex located at the origin,
00109 // move the vertex to the specified position, smooth the mesh,
00110 // and verify that the vertex was moved back to the origin by
00111 // the smoother.
00112 bool smooth_mesh( MeshImpl* mesh,
00113                   Mesh* ref_mesh,
00114                   Mesh::VertexHandle vertex_8,
00115                   Mesh::VertexHandle vertex_9,
00116                   Vector3D delta,
00117                   QualityMetric* metric );
00118 
00119 const unsigned NUM_ELEM      = 6;
00120 const unsigned NUM_VERT      = 14;
00121 const unsigned VERT_PER_ELEM = 6;
00122 
00123 int main( int argc, char*[] )
00124 {
00125     if( argc != 1 )
00126     {
00127         std::cerr << "Invalid arguments.\n";
00128         return 2;
00129     }
00130 
00131     std::string meshfile = TestDir + "unittest/mesquite/3D/vtk/prisms/untangled/6-wedge-prism.vtk";
00132     unsigned i;
00133 
00134     MBMesquite::MsqPrintError err( cout );
00135     IdealWeightMeanRatio m1;
00136     IdealWeightInverseMeanRatio m2( err );
00137     ConditionNumberQualityMetric m3;
00138     QualityMetric* metrics[] = { &m1, &m2, &m3, 0 };
00139 
00140     // Read Mesh
00141     MBMesquite::MeshImpl mesh;
00142     mesh.read_vtk( meshfile.c_str(), err );
00143     CPPUNIT_ASSERT( !err );
00144     MBMesquite::MeshImpl ideal_mesh;
00145     ideal_mesh.read_vtk( meshfile.c_str(), err );
00146     CPPUNIT_ASSERT( !err );
00147 
00148     // Check that the mesh read correctly, and contains what is
00149     // expected later.
00150 
00151     // Get mesh data
00152     // Expecting file to contain 6 wedge elements constructed
00153     // from 14 vertices.
00154     std::vector< Mesh::VertexHandle > vert_array;
00155     std::vector< Mesh::ElementHandle > elem_array;
00156     std::vector< size_t > conn_offsets;
00157     mesh.get_all_elements( elem_array, err );
00158     CPPUNIT_ASSERT( !err );
00159     CPPUNIT_ASSERT_EQUAL( elem_array.size(), NUM_ELEM );
00160     mesh.elements_get_attached_vertices( arrptr( elem_array ), elem_array.size(), vert_array, conn_offsets, err );
00161     CPPUNIT_ASSERT( !err );
00162     CPPUNIT_ASSERT_EQUAL( vert_array.size(), VERT_PER_ELEM * NUM_ELEM );
00163     CPPUNIT_ASSERT_EQUAL( conn_offsets.size(), NUM_ELEM + 1 );
00164     EntityTopology type_array[NUM_ELEM];
00165     mesh.elements_get_topologies( arrptr( elem_array ), type_array, NUM_ELEM, err );
00166     CPPUNIT_ASSERT( !err );
00167 
00168     // Verify element types and number of vertices
00169     for( i = 0; i < NUM_ELEM; ++i )
00170     {
00171         CPPUNIT_ASSERT_EQUAL( type_array[i], PRISM );
00172         CPPUNIT_ASSERT_EQUAL( conn_offsets[i], VERT_PER_ELEM * i );
00173     }
00174 
00175     // All wedges should share the 9th and 10th vertices
00176     const unsigned INDEX_8 = 1, INDEX_9 = 4;
00177     Mesh::VertexHandle handle8 = vert_array[INDEX_8];
00178     Mesh::VertexHandle handle9 = vert_array[INDEX_9];
00179     for( i = 1; i < NUM_ELEM; ++i )
00180     {
00181         CPPUNIT_ASSERT_EQUAL( vert_array[VERT_PER_ELEM * i + INDEX_8], handle8 );
00182         CPPUNIT_ASSERT_EQUAL( vert_array[VERT_PER_ELEM * i + INDEX_9], handle9 );
00183     }
00184 
00185     // The input file should be a hexagonal prism decomposed into
00186     // 6 wedges such that all wedges have one quad face on the
00187     // boundary and all wedges share a single edge which is the axis
00188     // of the prism.
00189     MsqVertex vertices[NUM_ELEM * VERT_PER_ELEM];
00190     mesh.vertices_get_coordinates( arrptr( vert_array ), vertices, NUM_ELEM * VERT_PER_ELEM, err );
00191     CPPUNIT_ASSERT( !err );
00192     for( i = 0; i < NUM_ELEM * VERT_PER_ELEM; ++i )
00193     {
00194         if( vert_array[i] == handle8 || vert_array[i] == handle9 )
00195         {
00196             CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, vertices[i][1], 1e-6 );
00197             CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, vertices[i][2], 1e-6 );
00198         }
00199         else
00200         {
00201             Vector3D xproj( vertices[i] );
00202             xproj[0] = 0;
00203             CPPUNIT_ASSERT_DOUBLES_EQUAL( 2.0, xproj.length(), 1e-6 );
00204         }
00205     }
00206 
00207     // Try smoothing w/out moving the free vertices and verify that
00208     // the smoother didn't move the vertex
00209     Vector3D delta( 0, 0, 0 );
00210     for( i = 0; metrics[i] != NULL; ++i )
00211         CPPUNIT_ASSERT( !smooth_mesh( &mesh, &ideal_mesh, handle8, handle9, delta, metrics[i] ) );
00212 
00213     // Now try moving the vertex and see if the smoother moves it back
00214     // to the origin
00215     delta.set( 0.1, 0.0, 0.1 );
00216     for( i = 0; metrics[i] != NULL; ++i )
00217         CPPUNIT_ASSERT( !smooth_mesh( &mesh, &ideal_mesh, handle8, handle9, delta, metrics[i] ) );
00218 
00219     // Now try moving the vertex further and see if the smoother moves it back
00220     // to the origin
00221     delta.set( 1.0, 0.0, 1.0 );
00222     for( i = 0; metrics[i] != NULL; ++i )
00223         CPPUNIT_ASSERT( !smooth_mesh( &mesh, &ideal_mesh, handle8, handle9, delta, metrics[i] ) );
00224 
00225     return 0;
00226 }
00227 
00228 bool smooth_mesh( MeshImpl* mesh,
00229                   Mesh*,
00230                   Mesh::VertexHandle vert1,
00231                   Mesh::VertexHandle vert2,
00232                   Vector3D delta,
00233                   QualityMetric* metric )
00234 {
00235     MBMesquite::MsqPrintError err( cout );
00236     const Vector3D origin( 0, 0, 0 );
00237 
00238     // print a little output so we know when we died
00239     std::cout << "**************************************************************************" << std::endl
00240               << "* Smoothing..." << std::endl
00241               << "* Metric: " << metric->get_name() << std::endl
00242               << "* Offset: " << delta
00243               << std::endl  //<<
00244               //"**************************************************************************"
00245               << std::endl;
00246 
00247     // Set free vertices to specified position
00248     Mesh::VertexHandle handles[] = { vert1, vert2 };
00249     MsqVertex coordinates[2];
00250     mesh->vertices_get_coordinates( handles, coordinates, 2, err );
00251     CPPUNIT_ASSERT( !err );
00252     Vector3D coord1 = coordinates[0] + delta;
00253     Vector3D coord2 = coordinates[1] + delta;
00254     mesh->vertex_set_coordinates( vert1, coord1, err );
00255     CPPUNIT_ASSERT( !err );
00256     mesh->vertex_set_coordinates( vert2, coord2, err );
00257     CPPUNIT_ASSERT( !err );
00258 
00259     // Create an InstructionQueue
00260     InstructionQueue Q;
00261 
00262     // Set up objective function
00263     LPtoPTemplate obj_func( metric, 1, err );
00264     CPPUNIT_ASSERT( !err );
00265 
00266     // Create solver
00267     FeasibleNewton solver( &obj_func );
00268     CPPUNIT_ASSERT( !err );
00269     solver.use_global_patch();
00270 
00271     // Set stoping criteria for solver
00272     TerminationCriterion tc_inner;
00273     tc_inner.add_absolute_vertex_movement( 1e-7 );
00274     solver.set_inner_termination_criterion( &tc_inner );
00275 
00276     TerminationCriterion tc_outer;
00277     tc_outer.add_iteration_limit( 1 );
00278     solver.set_outer_termination_criterion( &tc_outer );
00279 
00280 #ifdef DO_QUALITY_ASSESSOR
00281     QualityAssessor qa( metric, 10 );
00282     Q.add_quality_assessor( &qa, err );
00283     CPPUNIT_ASSERT( !err );
00284 #endif
00285 
00286     // Add solver to queue
00287     Q.set_master_quality_improver( &solver, err );
00288     CPPUNIT_ASSERT( !err );
00289 
00290 #ifdef DO_QUALITY_ASSESSOR
00291     Q.add_quality_assessor( &qa, err );
00292     CPPUNIT_ASSERT( !err );
00293 #endif
00294 
00295     // And smooth...
00296     Q.run_instructions( mesh, err );
00297     CPPUNIT_ASSERT( !err );
00298 
00299     // Verify that vertices were moved back to origin
00300     MsqVertex new_coords[2];
00301     mesh->vertices_get_coordinates( handles, new_coords, 2, err );
00302     CPPUNIT_ASSERT( !err );
00303 
00304     // print a little output so we know when we died
00305     std::cout  //<<
00306                //"**************************************************************************"
00307         << std::endl
00308         << "* Done Smoothing:" << std::endl
00309         << "* Metric: " << metric->get_name() << std::endl
00310         << "* Position1: " << new_coords[0][0] << " " << new_coords[0][1] << " " << new_coords[0][2] << std::endl
00311         << "* Position2: " << new_coords[1][0] << " " << new_coords[1][1] << " " << new_coords[1][2] << std::endl
00312         << "**************************************************************************" << std::endl;
00313 
00314     CPPUNIT_ASSERT_VECTORS_EQUAL( coordinates[0], new_coords[0], TOL );
00315     CPPUNIT_ASSERT_VECTORS_EQUAL( coordinates[1], new_coords[1], TOL );
00316     return false;
00317 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines