MOAB: Mesh Oriented datABase
(version 5.2.1)
|
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 }