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 00026 ***************************************************************** */ 00027 // -*- Mode : c++; tab-width: 3; c-tab-always-indent: t; indent-tabs-mode: nil; c-basic-offset: 3 00028 // -*- 00029 // 00030 // SUMMARY: 00031 // USAGE: 00032 // 00033 // ORIG-DATE: 19-Feb-02 at 10:57:52 00034 // LAST-MOD: 23-Jul-03 at 18:09:51 by Thomas Leurent 00035 // 00036 // 00037 // DESCRIPTION: 00038 // ============ 00039 /*! \file main.cpp 00040 00041 describe main.cpp here 00042 00043 */ 00044 // DESCRIP-END. 00045 // 00046 00047 #include <iostream> 00048 using std::cout; 00049 using std::endl; 00050 #include <cstdlib> 00051 00052 #include "TestUtil.hpp" 00053 #include "Mesquite.hpp" 00054 #include "MeshImpl.hpp" 00055 #include "MsqError.hpp" 00056 #include "Vector3D.hpp" 00057 #include "InstructionQueue.hpp" 00058 #include "PatchData.hpp" 00059 #include "TerminationCriterion.hpp" 00060 #include "QualityAssessor.hpp" 00061 //#include "QuadLagrangeShape.hpp" 00062 00063 // algorithms 00064 #include "IdealShapeTarget.hpp" 00065 #include "TQualityMetric.hpp" 00066 #include "TInverseMeanRatio.hpp" 00067 #include "ConditionNumberQualityMetric.hpp" 00068 #include "LPtoPTemplate.hpp" 00069 #include "LInfTemplate.hpp" 00070 #include "SteepestDescent.hpp" 00071 #include "ConjugateGradient.hpp" 00072 00073 #include "PlanarDomain.hpp" 00074 using namespace MBMesquite; 00075 00076 /* This is the input mesh topology 00077 (0)------(16)-----(1)------(17)-----(2)------(18)-----(3) 00078 | | | | 00079 | | | | 00080 | | | | 00081 | | | | 00082 (19) 0 (20) 1 (21) 2 (22) 00083 | | | | 00084 | | | | 00085 | | | | 00086 | | | | 00087 (4)------(23)-----(5)------(24)-----(6)------(25)-----(7) 00088 | | | | 00089 | | | | 00090 | | | | 00091 | | | | 00092 (26) 3 (27) 4 (28) 5 (29) 00093 | | | | 00094 | | | | 00095 | | | | 00096 | | | | 00097 (8)------(30)-----(9)------(31)-----(10)-----(32)-----(11) 00098 | | | | 00099 | | | | 00100 | | | | 00101 | | | | 00102 (33) 6 (34) 7 (35) 8 (36) 00103 | | | | 00104 | | | | 00105 | | | | 00106 | | | | 00107 (12)-----(37)-----(13)-----(38)-----(14)-----(39)-----(15) 00108 */ 00109 00110 std::string LINEAR_INPUT_FILE_NAME = std::string( STRINGIFY( SRCDIR ) ) + "/linear_input.vtk"; 00111 std::string QUADRATIC_INPUT_FILE_NAME = std::string( STRINGIFY( SRCDIR ) ) + "/quadratic_input.vtk"; 00112 std::string EXPECTED_LINAR_FILE_NAME = std::string( STRINGIFY( SRCDIR ) ) + "/expected_linear_output.vtk"; 00113 std::string EXPECTED_QUADRATIC_FILE_NAME = std::string( STRINGIFY( SRCDIR ) ) + "/expected_quadratic_output.vtk"; 00114 std::string HOUR_INPUT_FILE_NAME = std::string( STRINGIFY( SRCDIR ) ) + "/hour-quad8.vtk"; 00115 std::string OUTPUT_FILE_NAME = "smoothed_qudratic_mesh.vtk"; 00116 const unsigned NUM_CORNER_VERTICES = 16; 00117 const unsigned NUM_MID_NODES = 24; 00118 const double SPATIAL_COMPARE_TOLERANCE = 4e-6; 00119 00120 void compare_nodes( size_t start_index, size_t end_index, Mesh* mesh1, Mesh* mesh2, MsqError& err ) 00121 { 00122 size_t i, num_verts = end_index - start_index; 00123 std::vector< MsqVertex > verts1( num_verts ), verts2( num_verts ); 00124 std::vector< Mesh::VertexHandle > handles1( num_verts ), handles2( num_verts ); 00125 00126 /* VertexIterator skips higher-order nodes. 00127 For now, just assume index == handle 00128 00129 std::vector<Mesh::VertexHandle>::iterator handle_iter1, handle_iter2; 00130 00131 // Skip start_index vertices 00132 VertexIterator* iter1 = mesh1->vertex_iterator( err ); MSQ_ERRRTN(err); 00133 VertexIterator* iter2 = mesh2->vertex_iterator( err ); MSQ_ERRRTN(err); 00134 for (i = 0; i < start_index; ++i) 00135 { 00136 if (iter1->is_at_end()) 00137 { 00138 MSQ_SETERR(err)("start index out of range for first mesh set", MsqError::INVALID_ARG); 00139 return; 00140 } 00141 if (iter2->is_at_end()) 00142 { 00143 MSQ_SETERR(err)("start index out of range for second mesh set", MsqError::INVALID_ARG); 00144 return; 00145 } 00146 iter1->operator++(); 00147 iter2->operator++(); 00148 } 00149 00150 // Get handles for vertices 00151 handle_iter1 = handles1.begin(); 00152 handle_iter2 = handles2.begin(); 00153 for (i = start_index; i < end_index; ++i) 00154 { 00155 if (iter1->is_at_end()) 00156 { 00157 MSQ_SETERR(err)("end index out of range for first mesh set", MsqError::INVALID_ARG); 00158 return; 00159 } 00160 *handle_iter1 = iter1->operator*(); 00161 iter1->operator++(); 00162 ++handle_iter1; 00163 00164 if (iter2->is_at_end()) 00165 { 00166 MSQ_SETERR(err)("end index out of range for second mesh set", MsqError::INVALID_ARG); 00167 return; 00168 } 00169 *handle_iter2 = iter2->operator*(); 00170 iter2->operator++(); 00171 ++handle_iter2; 00172 } 00173 */ 00174 for( i = start_index; i < end_index; ++i ) 00175 handles1[i - start_index] = handles2[i - start_index] = (void*)i; 00176 00177 // Get coordinates from handles 00178 mesh1->vertices_get_coordinates( arrptr( handles1 ), arrptr( verts1 ), num_verts, err );MSQ_ERRRTN( err ); 00179 mesh2->vertices_get_coordinates( arrptr( handles2 ), arrptr( verts2 ), num_verts, err );MSQ_ERRRTN( err ); 00180 00181 // Compare coordinates 00182 for( i = 0; i < num_verts; ++i ) 00183 { 00184 const double diff = ( verts1[i] - verts2[i] ).length(); 00185 if( diff > SPATIAL_COMPARE_TOLERANCE ) 00186 { 00187 MSQ_SETERR( err ) 00188 ( MsqError::INTERNAL_ERROR, "%u%s vertices differ. (%f,%f,%f) vs (%f,%f,%f)", (unsigned)( 1 + i ), 00189 i % 10 == 0 ? "st" : i % 10 == 1 ? "nd" : i % 10 == 2 ? "rd" : "th", verts1[i][0], verts1[i][1], 00190 verts1[i][2], verts2[i][0], verts2[i][1], verts2[i][2] ); 00191 return; 00192 } 00193 } 00194 } 00195 00196 // code copied from testSuite/algorithm_test/main.cpp 00197 InstructionQueue* create_instruction_queue( MsqError& err ) 00198 { 00199 // creates an intruction queue 00200 InstructionQueue* queue1 = new InstructionQueue; 00201 00202 // creates a mean ratio quality metric ... 00203 // IdealWeightInverseMeanRatio* mean = new IdealWeightInverseMeanRatio(err); MSQ_ERRZERO(err); 00204 TargetCalculator* tc = new IdealShapeTarget; 00205 TQualityMetric* mean = new TQualityMetric( tc, 0, new TInverseMeanRatio ); 00206 00207 LPtoPTemplate* obj_func = new LPtoPTemplate( mean, 1, err ); 00208 MSQ_ERRZERO( err ); 00209 00210 // creates the optimization procedures 00211 // ConjugateGradient* pass1 = new ConjugateGradient( obj_func, err ); 00212 SteepestDescent* pass1 = new SteepestDescent( obj_func ); 00213 00214 // perform optimization globally 00215 pass1->use_global_patch(); 00216 00217 // QualityAssessor* mean_qa = new QualityAssessor(mean); 00218 00219 //**************Set termination criterion**************** 00220 00221 // perform 1 pass of the outer loop (this line isn't essential as it is 00222 // the default behavior). 00223 TerminationCriterion* tc_outer = new TerminationCriterion; 00224 tc_outer->add_iteration_limit( 1 ); 00225 pass1->set_outer_termination_criterion( tc_outer ); 00226 00227 // perform the inner loop until a certain objective function value is 00228 // reached. The exact value needs to be determined (about 18095). 00229 // As a safety, also stop if the time exceeds 10 minutes (600 seconds). 00230 TerminationCriterion* tc_inner = new TerminationCriterion; 00231 tc_inner->add_absolute_vertex_movement( 1e-6 ); 00232 pass1->set_inner_termination_criterion( tc_inner ); 00233 00234 // adds 1 pass of pass1 to mesh_set1 00235 // queue1->add_quality_assessor(mean_qa,err); MSQ_ERRZERO(err); 00236 queue1->set_master_quality_improver( pass1, err ); 00237 MSQ_ERRZERO( err ); 00238 // queue1->add_quality_assessor(mean_qa,err); MSQ_ERRZERO(err); 00239 00240 return queue1; 00241 } 00242 00243 int do_test( bool slave ) 00244 { 00245 MsqPrintError err( cout ); 00246 // QuadLagrangeShape quad9; 00247 00248 // Create geometry 00249 Vector3D z( 0, 0, 1 ), o( 0, 0, 0 ); 00250 PlanarDomain geom( z, o ); 00251 00252 // Read in linear input mesh 00253 cout << "Reading " << LINEAR_INPUT_FILE_NAME << endl; 00254 MeshImpl* linear_in = new MeshImpl; 00255 linear_in->read_vtk( LINEAR_INPUT_FILE_NAME.c_str(), err ); 00256 if( MSQ_CHKERR( err ) ) return 1; 00257 00258 // Read in expected linear results 00259 cout << "Reading " << EXPECTED_LINAR_FILE_NAME << endl; 00260 MeshImpl* linear_ex = new MeshImpl; 00261 linear_ex->read_vtk( EXPECTED_LINAR_FILE_NAME.c_str(), err ); 00262 if( MSQ_CHKERR( err ) ) return 1; 00263 00264 // Read in second copy of quadratic input mesh 00265 cout << "Reading " << QUADRATIC_INPUT_FILE_NAME << " again" << endl; 00266 MeshImpl* quadratic_in_2 = new MeshImpl; 00267 quadratic_in_2->read_vtk( QUADRATIC_INPUT_FILE_NAME.c_str(), err ); 00268 if( MSQ_CHKERR( err ) ) return 1; 00269 00270 // Read in expected quadratic results 00271 cout << "Reading " << EXPECTED_QUADRATIC_FILE_NAME << endl; 00272 MeshImpl* quadratic_ex = new MeshImpl; 00273 quadratic_ex->read_vtk( EXPECTED_QUADRATIC_FILE_NAME.c_str(), err ); 00274 if( MSQ_CHKERR( err ) ) return 1; 00275 00276 // Smooth linear mesh and check results 00277 cout << "Smoothing linear elements" << endl; 00278 InstructionQueue* q1 = create_instruction_queue( err ); 00279 if( MSQ_CHKERR( err ) ) return 1; 00280 MeshDomainAssoc* mesh_and_domain = new MeshDomainAssoc( linear_in, &geom ); 00281 q1->run_instructions( mesh_and_domain, err ); 00282 if( MSQ_CHKERR( err ) ) return 1; 00283 delete mesh_and_domain; 00284 cout << "Checking results" << endl; 00285 compare_nodes( 0, NUM_CORNER_VERTICES, linear_in, linear_ex, err ); 00286 if( MSQ_CHKERR( err ) ) 00287 { 00288 MsqError tmperr; 00289 linear_in->write_vtk( "bad_mesh.vtk", tmperr ); 00290 return 1; 00291 } 00292 delete q1; 00293 delete linear_in; 00294 00295 // Smooth corner vertices and adjust mid-side nodes 00296 cout << "Smoothing quadratic elements" << endl; 00297 InstructionQueue* q3 = create_instruction_queue( err ); 00298 if( MSQ_CHKERR( err ) ) return 1; 00299 if( !slave ) q3->set_slaved_ho_node_mode( Settings::SLAVE_NONE ); 00300 // q3->set_mapping_function( &quad9 ); 00301 MeshDomainAssoc* mesh_and_domain2 = new MeshDomainAssoc( quadratic_in_2, &geom ); 00302 q3->run_instructions( mesh_and_domain2, err ); 00303 if( MSQ_CHKERR( err ) ) return 1; 00304 delete mesh_and_domain2; 00305 // Make sure corner vertices are the same as in the linear case 00306 cout << "Checking results" << endl; 00307 compare_nodes( 0, NUM_CORNER_VERTICES, quadratic_in_2, linear_ex, err ); 00308 if( MSQ_CHKERR( err ) ) return 1; 00309 delete linear_ex; 00310 00311 // Make sure mid-side vertices are updated correctly 00312 compare_nodes( NUM_CORNER_VERTICES, NUM_CORNER_VERTICES + NUM_MID_NODES, quadratic_in_2, quadratic_ex, err ); 00313 if( MSQ_CHKERR( err ) ) 00314 { 00315 MsqError tmperr; 00316 quadratic_in_2->write_vtk( "bad_mesh.vtk", tmperr ); 00317 return 1; 00318 } 00319 delete q3; 00320 delete quadratic_in_2; 00321 delete quadratic_ex; 00322 00323 if( MSQ_CHKERR( err ) ) return 1; 00324 return 0; 00325 } 00326 00327 int do_smooth_ho() 00328 { 00329 MsqPrintError err( cout ); 00330 // QuadLagrangeShape quad9; 00331 00332 // Create geometry 00333 PlanarDomain geom( PlanarDomain::XY ); 00334 00335 // Read in one copy of quadratic input mesh 00336 cout << "Reading " << HOUR_INPUT_FILE_NAME << endl; 00337 MeshImpl* quadratic_in = new MeshImpl; 00338 quadratic_in->read_vtk( HOUR_INPUT_FILE_NAME.c_str(), err ); 00339 if( MSQ_CHKERR( err ) ) return 1; 00340 00341 // Read in expected results 00342 // cout << "Reading " << HOUR_EXPECTED_FILE_NAME << endl; 00343 // MeshImpl* quadratic_ex = new MeshImpl; 00344 // quadratic_ex->read_vtk( results, err ); 00345 // if (MSQ_CHKERR(err)) return 1; 00346 00347 // Smooth linear mesh and check results 00348 cout << "Smoothing higher-order nodes" << endl; 00349 InstructionQueue* q1 = create_instruction_queue( err ); 00350 if( MSQ_CHKERR( err ) ) return 1; 00351 q1->set_slaved_ho_node_mode( Settings::SLAVE_NONE ); 00352 // q1->set_mapping_function( &quad9 ); 00353 MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( quadratic_in, &geom ); 00354 q1->run_instructions( &mesh_and_domain, err ); 00355 if( MSQ_CHKERR( err ) ) return 1; 00356 cout << "Checking results" << endl; 00357 // compare_nodes( 0, NUM_CORNER_VERTICES + NUM_MID_NODES, 00358 // quadratic_in, quadratic_ex, err ); 00359 if( MSQ_CHKERR( err ) ) 00360 { 00361 MsqError tmperr; 00362 quadratic_in->write_vtk( "bad_mesh.vtk", tmperr ); 00363 return 1; 00364 } 00365 delete q1; 00366 quadratic_in->write_vtk( "smooth_ho.vtk", err ); 00367 delete quadratic_in; 00368 00369 if( MSQ_CHKERR( err ) ) return 1; 00370 return 0; 00371 } 00372 00373 int main() 00374 { 00375 cout << "Running test with all higher-order nodes slaved." << endl; 00376 int result1 = do_test( true ); 00377 cout << "Running test with no higher-order nodes slaved." << endl; 00378 int result2 = do_test( false ); 00379 cout << "Running test with only ho-nodes free." << endl; 00380 int result3 = do_smooth_ho(); 00381 return result1 + result2 + result3; 00382 }