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: 2; c-tab-always-indent: t; indent-tabs-mode: nil; c-basic-offset: 2 00028 // -*- 00029 // 00030 // SUMMARY: 00031 // USAGE: 00032 // 00033 // ORIG-DATE: 24-Jan-12 00034 // LAST-MOD: 26-Jan-1 by Stephen Kennon 00035 // 00036 // 00037 // DESCRIPTION: 00038 // ============ 00039 /*! \file par_hex.cpp 00040 00041 A test of MBMesquite's parallel capabilities. Reads a split vtk file, smooths in parallel using 00042 Laplace smoothing, writes out the result (which can be compared with the "gold" copy of the same 00043 name in the meshFiles VTK directory). 00044 00045 See the MBMesquite User's Guide, section "Using MBMesquite in Parallel" - this code is very similar 00046 to the example code shown therein. 00047 00048 */ 00049 // DESCRIP-END. 00050 // 00051 00052 #include "MeshImpl.hpp" 00053 #include "MeshUtil.hpp" 00054 #include "MsqTimer.hpp" 00055 #include "Mesquite.hpp" 00056 #include "MsqError.hpp" 00057 #include "Vector3D.hpp" 00058 #include "InstructionQueue.hpp" 00059 #include "LaplaceWrapper.hpp" 00060 #include "ShapeImprovementWrapper.hpp" 00061 #include "UntangleWrapper.hpp" 00062 #include "PatchData.hpp" 00063 #include "TerminationCriterion.hpp" 00064 #include "QualityAssessor.hpp" 00065 00066 /* MBMesquite includes */ 00067 #include "ParallelMeshImpl.hpp" 00068 #include "ParallelHelper.hpp" 00069 #include "MsqDebug.hpp" 00070 #include "Settings.hpp" 00071 //#include "ShapeImprovementWrapper.hpp" 00072 //#include "UntangleWrapper.hpp" 00073 00074 #include "IdealWeightInverseMeanRatio.hpp" 00075 00076 #include "UntangleBetaQualityMetric.hpp" 00077 #include "LPtoPTemplate.hpp" 00078 #include "ConjugateGradient.hpp" 00079 #include "SteepestDescent.hpp" 00080 //#include "FeasibleNewton.hpp" 00081 00082 // algorithms 00083 #include "Randomize.hpp" 00084 #include "ConditionNumberQualityMetric.hpp" 00085 #include "UntangleBetaQualityMetric.hpp" 00086 #include "LPtoPTemplate.hpp" 00087 #include "LInfTemplate.hpp" 00088 #include "SteepestDescent.hpp" 00089 #include "ConjugateGradient.hpp" 00090 #include "PlanarDomain.hpp" 00091 #include "TestUtil.hpp" 00092 00093 #include <iostream> 00094 using std::cout; 00095 using std::endl; 00096 #include <cstdlib> 00097 00098 #include <mpi.h> 00099 #include <sstream> 00100 00101 using namespace MBMesquite; 00102 00103 #define VTK_3D_DIR ( TestDir + "/3D/vtk/hexes/tangled/" ) 00104 #define VTK_2D_DIR ( TestDir + "/2D/vtk/quads/tangled/" ) 00105 00106 using namespace std; 00107 00108 const double DEF_UNT_BETA = 1e-8; 00109 const double DEF_SUC_EPS = 1e-4; 00110 00111 class ParShapeImprover 00112 { 00113 int innerIter; 00114 double gradNorm; 00115 00116 public: 00117 ParShapeImprover( int inner_iterations = 100, double grad_norm = 1.e-8 ) 00118 : innerIter( inner_iterations ), gradNorm( grad_norm ) 00119 { 00120 } 00121 00122 class ParShapeImprovementWrapper : public Wrapper 00123 { 00124 00125 public: 00126 // Constructor sets the instructions in the queue. 00127 ParShapeImprovementWrapper( int inner_iterations = 100, double cpu_time = 0.0, double grad_norm = 1.e-8, 00128 int parallel_iterations = 10 ) 00129 : innerIter( inner_iterations ), maxTime( cpu_time ), gradNorm( grad_norm ), untBeta( DEF_UNT_BETA ), 00130 successiveEps( DEF_SUC_EPS ), parallelIterations( parallel_iterations ), m_do_untangle_only( false ) 00131 { 00132 } 00133 00134 protected: 00135 void run_wrapper( MeshDomainAssoc* mesh_and_domain, ParallelMesh* pmesh, Settings* settings, 00136 QualityAssessor* qa, MsqError& err ); 00137 00138 private: 00139 int innerIter; 00140 double maxTime, gradNorm; 00141 // constants 00142 const double untBeta; 00143 const double successiveEps; 00144 int parallelIterations; 00145 00146 public: 00147 bool m_do_untangle_only; 00148 }; 00149 00150 static int count_invalid_elements( Mesh& mesh, MeshDomain* domain = 0 ); 00151 00152 void run( Mesh& mesh, MeshDomain* domain, MsqError& err, bool always_smooth = true, int debug = 0 ); 00153 }; 00154 00155 void ParShapeImprover::ParShapeImprovementWrapper::run_wrapper( MeshDomainAssoc* mesh_and_domain, ParallelMesh* pmesh, 00156 Settings* settings, QualityAssessor* qa, MsqError& err ) 00157 { 00158 int rank, nprocs; 00159 MPI_Comm_rank( MPI_COMM_WORLD, &rank ); 00160 MPI_Comm_size( MPI_COMM_WORLD, &nprocs ); 00161 00162 // Define an untangler 00163 // UntangleBetaQualityMetric untangle_metric( untBeta ); 00164 UntangleBetaQualityMetric untangle_metric( 1.e-6 ); 00165 00166 bool check_untangle = true; 00167 if( check_untangle ) 00168 { 00169 std::cout << "\nP[" << rank << "] ParShapeImprover.... running QA with untangle_metric before... " 00170 << std::endl; 00171 InstructionQueue q1; 00172 QualityAssessor qa_untangle( &untangle_metric ); 00173 q1.add_quality_assessor( &qa_untangle, err );MSQ_ERRRTN( err ); 00174 q1.run_common( mesh_and_domain, pmesh, settings, err ); 00175 std::cout << "\nP[" << rank << "] ParShapeImprover.... running QA with untangle_metric... before... done " 00176 << std::endl; 00177 } 00178 00179 LPtoPTemplate untangle_func( 2, &untangle_metric ); 00180 ConjugateGradient untangle_solver( &untangle_func ); 00181 // untangle_solver.set_debugging_level(3); 00182 00183 // SteepestDescent untangle_solver( &untangle_func ); 00184 TerminationCriterion untangle_inner( "<type:untangle_inner>" ), untangle_outer( "<type:untangle_outer>" ); 00185 untangle_solver.use_global_patch(); 00186 00187 // untangle_inner.add_absolute_gradient_L2_norm( gradNorm ); 00188 // untangle_inner.add_absolute_successive_improvement( successiveEps ); 00189 // untangle_inner.add_relative_successive_improvement( 1.e-6 ); 00190 // untangle_inner.add_untangled_mesh(); 00191 00192 untangle_inner.write_iterations( "untangle.gpt", err ); 00193 00194 // For parallel runs, we generally need to have the inner and outer TerminationCriterion 00195 // have the same criteria else we can get an infinite loop (see VertexMover::loop_over_mesh) 00196 untangle_inner.add_absolute_quality_improvement( 0.0 ); 00197 untangle_inner.add_iteration_limit( 20 ); 00198 00199 untangle_outer.add_absolute_quality_improvement( 0.0 ); 00200 untangle_outer.add_iteration_limit( pmesh ? parallelIterations : 1 ); 00201 00202 untangle_solver.set_inner_termination_criterion( &untangle_inner ); 00203 untangle_solver.set_outer_termination_criterion( &untangle_outer ); 00204 00205 // define shape improver 00206 IdealWeightInverseMeanRatio inverse_mean_ratio; 00207 inverse_mean_ratio.set_averaging_method( QualityMetric::LINEAR ); 00208 LPtoPTemplate obj_func( 2, &inverse_mean_ratio ); 00209 00210 ConjugateGradient shape_solver( &obj_func ); 00211 TerminationCriterion term_inner( "<type:shape_inner>" ), term_outer( "<type:shape_outer>" ); 00212 term_inner.write_iterations( "shape.gpt", err ); 00213 00214 shape_solver.use_global_patch(); 00215 qa->add_quality_assessment( &inverse_mean_ratio ); 00216 00217 // For parallel runs, we generally need to have the inner and outer TerminationCriterion 00218 // have the same criteria else we can get an infinite loop (see VertexMover::loop_over_mesh) 00219 term_inner.add_absolute_gradient_L2_norm( gradNorm ); 00220 term_inner.add_absolute_vertex_movement( 0.0 ); 00221 term_inner.add_iteration_limit( innerIter ); 00222 00223 term_outer.add_absolute_gradient_L2_norm( gradNorm ); 00224 term_outer.add_absolute_vertex_movement( 0.0 ); 00225 term_outer.add_iteration_limit( pmesh ? parallelIterations : 1 ); 00226 00227 // term_outer.add_absolute_quality_improvement( 1.e-6 ); 00228 //! term_outer.add_relative_successive_improvement( successiveEps ); 00229 00230 shape_solver.set_inner_termination_criterion( &term_inner ); 00231 shape_solver.set_outer_termination_criterion( &term_outer ); 00232 00233 // Apply CPU time limit to untangler 00234 if( maxTime > 0.0 ) untangle_inner.add_cpu_time( maxTime ); 00235 00236 Timer totalTimer; 00237 00238 // Run untangler 00239 std::cout << "\nP[" << rank << "] " 00240 << " ParShapeImprovementWrapper: running untangler...\n " << std::endl; 00241 bool use_untangle_wrapper = false; 00242 if( use_untangle_wrapper ) 00243 { 00244 UntangleWrapper uw; 00245 // uw.set_untangle_metric(UntangleWrapper::BETA); 00246 uw.run_instructions( mesh_and_domain, err ); 00247 } 00248 else 00249 { 00250 InstructionQueue q1; 00251 QualityAssessor qa_untangle( &untangle_metric ); 00252 q1.add_quality_assessor( &qa_untangle, err );MSQ_ERRRTN( err ); 00253 q1.set_master_quality_improver( &untangle_solver, err );MSQ_ERRRTN( err ); 00254 q1.add_quality_assessor( &qa_untangle, err );MSQ_ERRRTN( err ); 00255 q1.run_common( mesh_and_domain, pmesh, settings, err ); 00256 } 00257 std::cout << "\nP[" << rank << "] " 00258 << " ParShapeImprovementWrapper: running untangler... done\n " << std::endl; 00259 std::cout << "\nP[" << rank << "] " 00260 << " ParShapeImprovementWrapper: MsqError after untangler: " << err << std::endl; 00261 00262 bool check_quality_after_untangler = true; 00263 if( check_quality_after_untangler ) 00264 { 00265 Mesh* mesh = mesh_and_domain->get_mesh(); 00266 MeshDomain* domain = mesh_and_domain->get_domain(); 00267 int num_invalid = count_invalid_elements( *mesh, domain ); 00268 std::cout << "\nP[" << rank << "] " 00269 << " ParShapeImprover num_invalid after untangler= " << num_invalid << " " 00270 << ( num_invalid ? " ERROR still have invalid elements after MBMesquite untangle" 00271 : " SUCCESS: untangled invalid elements " ) 00272 << std::endl; 00273 00274 if( check_untangle ) 00275 { 00276 std::cout << "\nP[" << rank << "] ParShapeImprover.... running QA with untangle_metric " << std::endl; 00277 InstructionQueue q1; 00278 QualityAssessor qa_untangle( &untangle_metric ); 00279 q1.add_quality_assessor( &qa_untangle, err );MSQ_ERRRTN( err ); 00280 q1.run_common( mesh_and_domain, pmesh, settings, err ); 00281 std::cout << "\nP[" << rank << "] ParShapeImprover.... running QA with untangle_metric... done " 00282 << std::endl; 00283 } 00284 00285 if( num_invalid ) return; 00286 } 00287 if( m_do_untangle_only ) return;MSQ_ERRRTN( err ); 00288 00289 // If limited by CPU time, limit next step to remaning time 00290 if( maxTime > 0.0 ) 00291 { 00292 double remaining = maxTime - totalTimer.since_birth(); 00293 if( remaining <= 0.0 ) 00294 { 00295 MSQ_DBGOUT( 2 ) << "Optimization is terminating without perfoming shape improvement." << std::endl; 00296 remaining = 0.0; 00297 } 00298 term_inner.add_cpu_time( remaining ); 00299 } 00300 00301 // Run shape improver 00302 InstructionQueue q2; 00303 std::cout << "\nP[" << rank << "] " 00304 << " ParShapeImprovementWrapper: running shape improver... \n" 00305 << std::endl; 00306 q2.add_quality_assessor( qa, err );MSQ_ERRRTN( err ); 00307 q2.set_master_quality_improver( &shape_solver, err );MSQ_ERRRTN( err ); 00308 q2.add_quality_assessor( qa, err );MSQ_ERRRTN( err ); 00309 q2.run_common( mesh_and_domain, pmesh, settings, err ); 00310 std::cout << "\nP[" << rank << "] " 00311 << " ParShapeImprovementWrapper: running shape improver... done \n" 00312 << std::endl;MSQ_ERRRTN( err ); 00313 } 00314 00315 int ParShapeImprover::count_invalid_elements( Mesh& mesh, MeshDomain* domain ) 00316 { 00317 MsqError err; 00318 InstructionQueue q; 00319 00320 IdealWeightInverseMeanRatio metric; 00321 metric.set_averaging_method( QualityMetric::LINEAR ); 00322 00323 // Check for inverted elements in the mesh 00324 QualityAssessor inv_check( &metric ); 00325 // inv_check.disable_printing_results(); 00326 q.add_quality_assessor( &inv_check, err ); 00327 MSQ_ERRZERO( err ); 00328 Settings settings; 00329 // bug? should we pass in pmesh? 00330 Mesh* mesh_ptr = &mesh; 00331 MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( mesh_ptr, domain ); 00332 q.run_common( &mesh_and_domain, 0, &settings, err ); 00333 MSQ_ERRZERO( err ); 00334 const QualityAssessor::Assessor* inv_b = inv_check.get_results( &metric ); 00335 int num_invalid = inv_b->get_invalid_element_count(); 00336 return num_invalid; 00337 } 00338 00339 void ParShapeImprover::run( Mesh& mesh, MeshDomain* domain, MsqError&, bool always_smooth, int debug ) 00340 { 00341 int rank, nprocs; 00342 MPI_Comm_rank( MPI_COMM_WORLD, &rank ); 00343 MPI_Comm_size( MPI_COMM_WORLD, &nprocs ); 00344 00345 if( debug ) 00346 { 00347 MsqDebug::enable( 1 ); 00348 if( debug > 1 ) MsqDebug::enable( 2 ); 00349 if( debug > 2 ) MsqDebug::enable( 3 ); 00350 } 00351 00352 ParallelMesh* pmesh = dynamic_cast< ParallelMesh* >( &mesh ); 00353 std::cout << "P[" << rank << "] " 00354 << " ParShapeImprover::run: pmesh= " << pmesh << std::endl; 00355 00356 MsqError mErr; 00357 int num_invalid = 0; 00358 bool check_quality = true; 00359 if( check_quality ) 00360 { 00361 num_invalid = count_invalid_elements( mesh, domain ); 00362 std::cout << "\nP[" << rank << "] " 00363 << " ParShapeImprover num_invalid before= " << num_invalid 00364 << ( num_invalid 00365 ? " WARNING: invalid elements exist before MBMesquite smoothing" 00366 : ( !always_smooth ? "WARNING: no smoothing requested since always_smooth=false" : " " ) ) 00367 << std::endl; 00368 } 00369 00370 if( num_invalid || always_smooth ) 00371 { 00372 bool use_canned_wrapper = false; 00373 if( use_canned_wrapper ) 00374 { 00375 ShapeImprovementWrapper siw( mErr ); 00376 if( pmesh ) 00377 siw.run_instructions( pmesh, domain, mErr ); 00378 else 00379 { 00380 MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &mesh, domain ); 00381 siw.run_instructions( &mesh_and_domain, mErr ); 00382 } 00383 } 00384 else 00385 { 00386 // int msq_debug = debug; // 1,2,3 for more debug info 00387 // bool always_smooth_local = false; 00388 00389 bool do_untangle_only = false; 00390 ParShapeImprover::ParShapeImprovementWrapper siw( innerIter, 0.0, gradNorm, 100 ); 00391 siw.m_do_untangle_only = do_untangle_only; 00392 if( pmesh ) 00393 siw.run_instructions( pmesh, domain, mErr ); 00394 else 00395 { 00396 MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &mesh, domain ); 00397 siw.run_instructions( &mesh_and_domain, mErr ); 00398 } 00399 } 00400 00401 std::cout << "\nP[" << rank << "] " 00402 << " ParShapeImprover: MsqError after ShapeImprovementWrapper: " << mErr << std::endl; 00403 00404 if( check_quality ) 00405 { 00406 num_invalid = count_invalid_elements( mesh, domain ); 00407 std::cout << "\nP[" << rank << "] " 00408 << " ParShapeImprover num_invalid after= " << num_invalid << " " 00409 << ( num_invalid ? " ERROR still have invalid elements after MBMesquite smoothing" 00410 : " SUCCESS: smoothed and removed invalid elements " ) 00411 << std::endl; 00412 } 00413 00414 MSQ_ERRRTN( mErr ); 00415 } 00416 } 00417 00418 static int test( std::string filename_prefix, std::string mesh_topology_name, MeshDomain* domain = 0 ) 00419 { 00420 int rank, nprocs; 00421 MPI_Comm_rank( MPI_COMM_WORLD, &rank ); 00422 MPI_Comm_size( MPI_COMM_WORLD, &nprocs ); 00423 00424 if( nprocs > 2 ) 00425 { 00426 cerr << "parallel_untangle_shape::test(" << mesh_topology_name << " can only be run with 1 or 2 processors" 00427 << std::endl; 00428 return 0; 00429 } 00430 00431 /* create processor-specific file names */ 00432 ostringstream in_name, out_name, gold_name; 00433 in_name << filename_prefix << "par_untangle_original_" << mesh_topology_name << "_mesh." << nprocs << "." << rank 00434 << ".vtk"; 00435 gold_name << filename_prefix << "par_untangle_smoothed_" << mesh_topology_name << "_mesh." << nprocs << "." << rank 00436 << ".vtk"; 00437 out_name << "par_untangle_smoothed_" << mesh_topology_name << "_mesh." << nprocs << "." << rank << ".vtk"; 00438 00439 cout << "in_name= " << in_name.str() << " gold_name= " << gold_name.str() << " out_name= " << out_name.str() 00440 << std::endl; 00441 00442 /* load different mesh files on each processor */ 00443 MsqError err; 00444 MeshImpl mesh; 00445 mesh.read_vtk( in_name.str().c_str(), err ); 00446 if( err ) 00447 { 00448 cerr << err << endl; 00449 return 1; 00450 } 00451 00452 /* create parallel mesh instance, specifying tags 00453 * containing parallel data */ 00454 ParallelMeshImpl parallel_mesh( &mesh, "GLOBAL_ID", "PROCESSOR_ID" ); 00455 ParallelHelperImpl helper; 00456 helper.set_communicator( MPI_COMM_WORLD ); 00457 helper.set_parallel_mesh( ¶llel_mesh ); 00458 parallel_mesh.set_parallel_helper( &helper ); 00459 00460 /* do Laplacian smooth */ 00461 // LaplaceWrapper optimizer; 00462 // optimizer.run_instructions(¶llel_mesh, err); 00463 00464 int msq_debug = 0; // 1,2,3 for more debug info 00465 bool always_smooth = true; 00466 int innerIter = 100; 00467 double gradNorm = 1.e-9; 00468 00469 ParShapeImprover si( innerIter, gradNorm ); 00470 // Mesh *pmesh = ¶llel_mesh; 00471 si.run( parallel_mesh, domain, err, always_smooth, msq_debug ); 00472 if( err ) 00473 { 00474 cerr << err << endl; 00475 return 1; 00476 } 00477 00478 /* write mesh */ 00479 mesh.write_vtk( out_name.str().c_str(), err ); 00480 if( err ) 00481 { 00482 cerr << err << endl; 00483 return 1; 00484 } 00485 00486 // std::cout << "P[ " << rank <<"] reading gold..." << std::endl; 00487 00488 /* compare mesh with gold copy */ 00489 MeshImpl gold; 00490 gold.read_vtk( gold_name.str().c_str(), err ); 00491 if( err ) 00492 { 00493 cerr << err << endl; 00494 return 1; 00495 } 00496 00497 // std::cout << "P[ " << rank <<"] read gold, checking mesh diff..." << std::endl; 00498 bool do_print = true; 00499 double tol = 1.e-4; 00500 bool diff = MeshUtil::meshes_are_different( mesh, gold, err, tol, do_print ); 00501 if( err ) 00502 { 00503 cerr << err << endl; 00504 return 1; 00505 } 00506 // std::cout << "P[ " << rank <<"] read gold, checking mesh diff...done" << std::endl; 00507 00508 if( diff ) 00509 { 00510 cerr << "Error, computed mesh different from gold copy" << std::endl; 00511 return 1; 00512 } 00513 00514 print_timing_diagnostics( cout ); 00515 00516 return 0; 00517 } 00518 00519 int main( int argc, char* argv[] ) 00520 { 00521 /* init MPI */ 00522 if( MPI_SUCCESS != MPI_Init( &argc, &argv ) ) 00523 { 00524 cerr << "MPI_Init failed." << endl; 00525 return 2; 00526 } 00527 00528 Vector3D pnt( 0, 0, 0 ); 00529 Vector3D s_norm( 0, 0, 1 ); 00530 PlanarDomain msq_geom( s_norm, pnt ); 00531 00532 int t1 = test( VTK_2D_DIR, "quad", &msq_geom ); 00533 if( t1 ) return t1; 00534 t1 = test( VTK_3D_DIR, "hex" ); 00535 if( t1 ) return t1; 00536 00537 MPI_Finalize(); 00538 return 0; 00539 }