MOAB: Mesh Oriented datABase  (version 5.2.1)
par_hex_untangle_shape.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 
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( &parallel_mesh );
00458     parallel_mesh.set_parallel_helper( &helper );
00459 
00460     /* do Laplacian smooth */
00461     // LaplaceWrapper optimizer;
00462     // optimizer.run_instructions(&parallel_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 = &parallel_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 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines