MOAB: Mesh Oriented datABase  (version 5.2.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 <stdlib.h>
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 true;                                          \
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 true;                                                          \
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 true;                                                          \
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 true;                                                          \
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, Mesh* ref_mesh, Mesh::VertexHandle vertex_8, Mesh::VertexHandle vertex_9,
00113                   Vector3D delta, QualityMetric* metric );
00114 
00115 const unsigned NUM_ELEM      = 6;
00116 const unsigned NUM_VERT      = 14;
00117 const unsigned VERT_PER_ELEM = 6;
00118 
00119 int main( int argc, char*[] )
00120 {
00121     if( argc != 1 )
00122     {
00123         std::cerr << "Invalid arguments.\n";
00124         return 2;
00125     }
00126 
00127     std::string meshfile = TestDir + "/3D/vtk/prisms/untangled/6-wedge-prism.vtk";
00128     unsigned i;
00129 
00130     MBMesquite::MsqPrintError err( cout );
00131     IdealWeightMeanRatio m1;
00132     IdealWeightInverseMeanRatio m2( err );
00133     ConditionNumberQualityMetric m3;
00134     QualityMetric* metrics[] = { &m1, &m2, &m3, 0 };
00135 
00136     // Read Mesh
00137     MBMesquite::MeshImpl mesh;
00138     mesh.read_vtk( meshfile.c_str(), err );
00139     CPPUNIT_ASSERT( !err );
00140     MBMesquite::MeshImpl ideal_mesh;
00141     ideal_mesh.read_vtk( meshfile.c_str(), err );
00142     CPPUNIT_ASSERT( !err );
00143 
00144     // Check that the mesh read correctly, and contains what is
00145     // expected later.
00146 
00147     // Get mesh data
00148     // Expecting file to contain 6 wedge elements constructed
00149     // from 14 vertices.
00150     std::vector< Mesh::VertexHandle > vert_array;
00151     std::vector< Mesh::ElementHandle > elem_array;
00152     std::vector< size_t > conn_offsets;
00153     mesh.get_all_elements( elem_array, err );
00154     CPPUNIT_ASSERT( !err );
00155     CPPUNIT_ASSERT_EQUAL( elem_array.size(), NUM_ELEM );
00156     mesh.elements_get_attached_vertices( arrptr( elem_array ), elem_array.size(), vert_array, conn_offsets, err );
00157     CPPUNIT_ASSERT( !err );
00158     CPPUNIT_ASSERT_EQUAL( vert_array.size(), VERT_PER_ELEM * NUM_ELEM );
00159     CPPUNIT_ASSERT_EQUAL( conn_offsets.size(), NUM_ELEM + 1 );
00160     EntityTopology type_array[NUM_ELEM];
00161     mesh.elements_get_topologies( arrptr( elem_array ), type_array, NUM_ELEM, err );
00162     CPPUNIT_ASSERT( !err );
00163 
00164     // Verify element types and number of vertices
00165     for( i = 0; i < NUM_ELEM; ++i )
00166     {
00167         CPPUNIT_ASSERT_EQUAL( type_array[i], PRISM );
00168         CPPUNIT_ASSERT_EQUAL( conn_offsets[i], VERT_PER_ELEM * i );
00169     }
00170 
00171     // All wedges should share the 9th and 10th vertices
00172     const unsigned INDEX_8 = 1, INDEX_9 = 4;
00173     Mesh::VertexHandle handle8 = vert_array[INDEX_8];
00174     Mesh::VertexHandle handle9 = vert_array[INDEX_9];
00175     for( i = 1; i < NUM_ELEM; ++i )
00176     {
00177         CPPUNIT_ASSERT_EQUAL( vert_array[VERT_PER_ELEM * i + INDEX_8], handle8 );
00178         CPPUNIT_ASSERT_EQUAL( vert_array[VERT_PER_ELEM * i + INDEX_9], handle9 );
00179     }
00180 
00181     // The input file should be a hexagonal prism decomposed into
00182     // 6 wedges such that all wedges have one quad face on the
00183     // boundary and all wedges share a single edge which is the axis
00184     // of the prism.
00185     MsqVertex vertices[NUM_ELEM * VERT_PER_ELEM];
00186     mesh.vertices_get_coordinates( arrptr( vert_array ), vertices, NUM_ELEM * VERT_PER_ELEM, err );
00187     CPPUNIT_ASSERT( !err );
00188     for( i = 0; i < NUM_ELEM * VERT_PER_ELEM; ++i )
00189     {
00190         if( vert_array[i] == handle8 || vert_array[i] == handle9 )
00191         {
00192             CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, vertices[i][1], 1e-6 );
00193             CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, vertices[i][2], 1e-6 );
00194         }
00195         else
00196         {
00197             Vector3D xproj( vertices[i] );
00198             xproj[0] = 0;
00199             CPPUNIT_ASSERT_DOUBLES_EQUAL( 2.0, xproj.length(), 1e-6 );
00200         }
00201     }
00202 
00203     // Try smoothing w/out moving the free vertices and verify that
00204     // the smoother didn't move the vertex
00205     Vector3D delta( 0, 0, 0 );
00206     for( i = 0; metrics[i] != NULL; ++i )
00207         CPPUNIT_ASSERT( !smooth_mesh( &mesh, &ideal_mesh, handle8, handle9, delta, metrics[i] ) );
00208 
00209     // Now try moving the vertex and see if the smoother moves it back
00210     // to the origin
00211     delta.set( 0.1, 0.0, 0.1 );
00212     for( i = 0; metrics[i] != NULL; ++i )
00213         CPPUNIT_ASSERT( !smooth_mesh( &mesh, &ideal_mesh, handle8, handle9, delta, metrics[i] ) );
00214 
00215     // Now try moving the vertex further and see if the smoother moves it back
00216     // to the origin
00217     delta.set( 1.0, 0.0, 1.0 );
00218     for( i = 0; metrics[i] != NULL; ++i )
00219         CPPUNIT_ASSERT( !smooth_mesh( &mesh, &ideal_mesh, handle8, handle9, delta, metrics[i] ) );
00220 
00221     return 0;
00222 }
00223 
00224 bool smooth_mesh( MeshImpl* mesh, Mesh*, Mesh::VertexHandle vert1, Mesh::VertexHandle vert2, Vector3D delta,
00225                   QualityMetric* metric )
00226 {
00227     MBMesquite::MsqPrintError err( cout );
00228     const Vector3D origin( 0, 0, 0 );
00229 
00230     // print a little output so we know when we died
00231     std::cout << "**************************************************************************" << std::endl
00232               << "* Smoothing..." << std::endl
00233               << "* Metric: " << metric->get_name() << std::endl
00234               << "* Offset: " << delta
00235               << std::endl  //<<
00236               //"**************************************************************************"
00237               << std::endl;
00238 
00239     // Set free vertices to specified position
00240     Mesh::VertexHandle handles[] = { vert1, vert2 };
00241     MsqVertex coordinates[2];
00242     mesh->vertices_get_coordinates( handles, coordinates, 2, err );
00243     CPPUNIT_ASSERT( !err );
00244     Vector3D coord1 = coordinates[0] + delta;
00245     Vector3D coord2 = coordinates[1] + delta;
00246     mesh->vertex_set_coordinates( vert1, coord1, err );
00247     CPPUNIT_ASSERT( !err );
00248     mesh->vertex_set_coordinates( vert2, coord2, err );
00249     CPPUNIT_ASSERT( !err );
00250 
00251     // Create an InstructionQueue
00252     InstructionQueue Q;
00253 
00254     // Set up objective function
00255     LPtoPTemplate obj_func( metric, 1, err );
00256     CPPUNIT_ASSERT( !err );
00257 
00258     // Create solver
00259     FeasibleNewton solver( &obj_func );
00260     CPPUNIT_ASSERT( !err );
00261     solver.use_global_patch();
00262 
00263     // Set stoping criteria for solver
00264     TerminationCriterion tc_inner;
00265     tc_inner.add_absolute_vertex_movement( 1e-7 );
00266     solver.set_inner_termination_criterion( &tc_inner );
00267 
00268     TerminationCriterion tc_outer;
00269     tc_outer.add_iteration_limit( 1 );
00270     solver.set_outer_termination_criterion( &tc_outer );
00271 
00272 #ifdef DO_QUALITY_ASSESSOR
00273     QualityAssessor qa( metric, 10 );
00274     Q.add_quality_assessor( &qa, err );
00275     CPPUNIT_ASSERT( !err );
00276 #endif
00277 
00278     // Add solver to queue
00279     Q.set_master_quality_improver( &solver, err );
00280     CPPUNIT_ASSERT( !err );
00281 
00282 #ifdef DO_QUALITY_ASSESSOR
00283     Q.add_quality_assessor( &qa, err );
00284     CPPUNIT_ASSERT( !err );
00285 #endif
00286 
00287     // And smooth...
00288     Q.run_instructions( mesh, err );
00289     CPPUNIT_ASSERT( !err );
00290 
00291     // Verify that vertices were moved back to origin
00292     MsqVertex new_coords[2];
00293     mesh->vertices_get_coordinates( handles, new_coords, 2, err );
00294     CPPUNIT_ASSERT( !err );
00295 
00296     // print a little output so we know when we died
00297     std::cout  //<<
00298                //"**************************************************************************"
00299         << std::endl
00300         << "* Done Smoothing:" << std::endl
00301         << "* Metric: " << metric->get_name() << std::endl
00302         << "* Position1: " << new_coords[0][0] << " " << new_coords[0][1] << " " << new_coords[0][2] << std::endl
00303         << "* Position2: " << new_coords[1][0] << " " << new_coords[1][1] << " " << new_coords[1][2] << std::endl
00304         << "**************************************************************************" << std::endl;
00305 
00306     CPPUNIT_ASSERT_VECTORS_EQUAL( coordinates[0], new_coords[0], TOL );
00307     CPPUNIT_ASSERT_VECTORS_EQUAL( coordinates[1], new_coords[1], TOL );
00308     return false;
00309 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines