MOAB: Mesh Oriented datABase  (version 5.3.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 + "unittest/mesquite/3D/vtk/hexes/tangled/" )
00104 #define VTK_2D_DIR ( TestDir + "unittest/mesquite/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,
00128                                     double cpu_time         = 0.0,
00129                                     double grad_norm        = 1.e-8,
00130                                     int parallel_iterations = 10 )
00131             : innerIter( inner_iterations ), maxTime( cpu_time ), gradNorm( grad_norm ), untBeta( DEF_UNT_BETA ),
00132               successiveEps( DEF_SUC_EPS ), parallelIterations( parallel_iterations ), m_do_untangle_only( false )
00133         {
00134         }
00135 
00136       protected:
00137         void run_wrapper( MeshDomainAssoc* mesh_and_domain,
00138                           ParallelMesh* pmesh,
00139                           Settings* settings,
00140                           QualityAssessor* qa,
00141                           MsqError& err );
00142 
00143       private:
00144         int innerIter;
00145         double maxTime, gradNorm;
00146         // constants
00147         const double untBeta;
00148         const double successiveEps;
00149         int parallelIterations;
00150 
00151       public:
00152         bool m_do_untangle_only;
00153     };
00154 
00155     static int count_invalid_elements( Mesh& mesh, MeshDomain* domain = 0 );
00156 
00157     void run( Mesh& mesh, MeshDomain* domain, MsqError& err, bool always_smooth = true, int debug = 0 );
00158 };
00159 
00160 void ParShapeImprover::ParShapeImprovementWrapper::run_wrapper( MeshDomainAssoc* mesh_and_domain,
00161                                                                 ParallelMesh* pmesh,
00162                                                                 Settings* settings,
00163                                                                 QualityAssessor* qa,
00164                                                                 MsqError& err )
00165 {
00166     int rank, nprocs;
00167     MPI_Comm_rank( MPI_COMM_WORLD, &rank );
00168     MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
00169 
00170     // Define an untangler
00171     // UntangleBetaQualityMetric untangle_metric( untBeta );
00172     UntangleBetaQualityMetric untangle_metric( 1.e-6 );
00173 
00174     bool check_untangle = true;
00175     if( check_untangle )
00176     {
00177         std::cout << "\nP[" << rank << "]  ParShapeImprover.... running QA with untangle_metric before... "
00178                   << std::endl;
00179         InstructionQueue q1;
00180         QualityAssessor qa_untangle( &untangle_metric );
00181         q1.add_quality_assessor( &qa_untangle, err );MSQ_ERRRTN( err );
00182         q1.run_common( mesh_and_domain, pmesh, settings, err );
00183         std::cout << "\nP[" << rank << "]  ParShapeImprover.... running QA with untangle_metric... before... done "
00184                   << std::endl;
00185     }
00186 
00187     LPtoPTemplate untangle_func( 2, &untangle_metric );
00188     ConjugateGradient untangle_solver( &untangle_func );
00189     // untangle_solver.set_debugging_level(3);
00190 
00191     // SteepestDescent untangle_solver( &untangle_func );
00192     TerminationCriterion untangle_inner( "<type:untangle_inner>" ), untangle_outer( "<type:untangle_outer>" );
00193     untangle_solver.use_global_patch();
00194 
00195     // untangle_inner.add_absolute_gradient_L2_norm( gradNorm );
00196     // untangle_inner.add_absolute_successive_improvement( successiveEps );
00197     // untangle_inner.add_relative_successive_improvement( 1.e-6 );
00198     // untangle_inner.add_untangled_mesh();
00199 
00200     untangle_inner.write_iterations( "untangle.gpt", err );
00201 
00202     // For parallel runs, we generally need to have the inner and outer TerminationCriterion
00203     //   have the same criteria else we can get an infinite loop (see VertexMover::loop_over_mesh)
00204     untangle_inner.add_absolute_quality_improvement( 0.0 );
00205     untangle_inner.add_iteration_limit( 20 );
00206 
00207     untangle_outer.add_absolute_quality_improvement( 0.0 );
00208     untangle_outer.add_iteration_limit( pmesh ? parallelIterations : 1 );
00209 
00210     untangle_solver.set_inner_termination_criterion( &untangle_inner );
00211     untangle_solver.set_outer_termination_criterion( &untangle_outer );
00212 
00213     // define shape improver
00214     IdealWeightInverseMeanRatio inverse_mean_ratio;
00215     inverse_mean_ratio.set_averaging_method( QualityMetric::LINEAR );
00216     LPtoPTemplate obj_func( 2, &inverse_mean_ratio );
00217 
00218     ConjugateGradient shape_solver( &obj_func );
00219     TerminationCriterion term_inner( "<type:shape_inner>" ), term_outer( "<type:shape_outer>" );
00220     term_inner.write_iterations( "shape.gpt", err );
00221 
00222     shape_solver.use_global_patch();
00223     qa->add_quality_assessment( &inverse_mean_ratio );
00224 
00225     // For parallel runs, we generally need to have the inner and outer TerminationCriterion
00226     //   have the same criteria else we can get an infinite loop (see VertexMover::loop_over_mesh)
00227     term_inner.add_absolute_gradient_L2_norm( gradNorm );
00228     term_inner.add_absolute_vertex_movement( 0.0 );
00229     term_inner.add_iteration_limit( innerIter );
00230 
00231     term_outer.add_absolute_gradient_L2_norm( gradNorm );
00232     term_outer.add_absolute_vertex_movement( 0.0 );
00233     term_outer.add_iteration_limit( pmesh ? parallelIterations : 1 );
00234 
00235     // term_outer.add_absolute_quality_improvement( 1.e-6 );
00236     //! term_outer.add_relative_successive_improvement( successiveEps );
00237 
00238     shape_solver.set_inner_termination_criterion( &term_inner );
00239     shape_solver.set_outer_termination_criterion( &term_outer );
00240 
00241     // Apply CPU time limit to untangler
00242     if( maxTime > 0.0 ) untangle_inner.add_cpu_time( maxTime );
00243 
00244     Timer totalTimer;
00245 
00246     // Run untangler
00247     std::cout << "\nP[" << rank << "] "
00248               << " ParShapeImprovementWrapper: running untangler...\n " << std::endl;
00249     bool use_untangle_wrapper = false;
00250     if( use_untangle_wrapper )
00251     {
00252         UntangleWrapper uw;
00253         // uw.set_untangle_metric(UntangleWrapper::BETA);
00254         uw.run_instructions( mesh_and_domain, err );
00255     }
00256     else
00257     {
00258         InstructionQueue q1;
00259         QualityAssessor qa_untangle( &untangle_metric );
00260         q1.add_quality_assessor( &qa_untangle, err );MSQ_ERRRTN( err );
00261         q1.set_master_quality_improver( &untangle_solver, err );MSQ_ERRRTN( err );
00262         q1.add_quality_assessor( &qa_untangle, err );MSQ_ERRRTN( err );
00263         q1.run_common( mesh_and_domain, pmesh, settings, err );
00264     }
00265     std::cout << "\nP[" << rank << "] "
00266               << " ParShapeImprovementWrapper: running untangler... done\n " << std::endl;
00267     std::cout << "\nP[" << rank << "] "
00268               << " ParShapeImprovementWrapper: MsqError after untangler: " << err << std::endl;
00269 
00270     bool check_quality_after_untangler = true;
00271     if( check_quality_after_untangler )
00272     {
00273         Mesh* mesh         = mesh_and_domain->get_mesh();
00274         MeshDomain* domain = mesh_and_domain->get_domain();
00275         int num_invalid    = count_invalid_elements( *mesh, domain );
00276         std::cout << "\nP[" << rank << "] "
00277                   << " ParShapeImprover num_invalid after untangler= " << num_invalid << " "
00278                   << ( num_invalid ? " ERROR still have invalid elements after MBMesquite untangle"
00279                                    : " SUCCESS: untangled invalid elements " )
00280                   << std::endl;
00281 
00282         if( check_untangle )
00283         {
00284             std::cout << "\nP[" << rank << "]  ParShapeImprover.... running QA with untangle_metric " << std::endl;
00285             InstructionQueue q1;
00286             QualityAssessor qa_untangle( &untangle_metric );
00287             q1.add_quality_assessor( &qa_untangle, err );MSQ_ERRRTN( err );
00288             q1.run_common( mesh_and_domain, pmesh, settings, err );
00289             std::cout << "\nP[" << rank << "]  ParShapeImprover.... running QA with untangle_metric... done "
00290                       << std::endl;
00291         }
00292 
00293         if( num_invalid ) return;
00294     }
00295     if( m_do_untangle_only ) return;MSQ_ERRRTN( err );
00296 
00297     // If limited by CPU time, limit next step to remaning time
00298     if( maxTime > 0.0 )
00299     {
00300         double remaining = maxTime - totalTimer.since_birth();
00301         if( remaining <= 0.0 )
00302         {
00303             MSQ_DBGOUT( 2 ) << "Optimization is terminating without perfoming shape improvement." << std::endl;
00304             remaining = 0.0;
00305         }
00306         term_inner.add_cpu_time( remaining );
00307     }
00308 
00309     // Run shape improver
00310     InstructionQueue q2;
00311     std::cout << "\nP[" << rank << "] "
00312               << " ParShapeImprovementWrapper: running shape improver... \n"
00313               << std::endl;
00314     q2.add_quality_assessor( qa, err );MSQ_ERRRTN( err );
00315     q2.set_master_quality_improver( &shape_solver, err );MSQ_ERRRTN( err );
00316     q2.add_quality_assessor( qa, err );MSQ_ERRRTN( err );
00317     q2.run_common( mesh_and_domain, pmesh, settings, err );
00318     std::cout << "\nP[" << rank << "] "
00319               << " ParShapeImprovementWrapper: running shape improver... done \n"
00320               << std::endl;MSQ_ERRRTN( err );
00321 }
00322 
00323 int ParShapeImprover::count_invalid_elements( Mesh& mesh, MeshDomain* domain )
00324 {
00325     MsqError err;
00326     InstructionQueue q;
00327 
00328     IdealWeightInverseMeanRatio metric;
00329     metric.set_averaging_method( QualityMetric::LINEAR );
00330 
00331     // Check for inverted elements in the mesh
00332     QualityAssessor inv_check( &metric );
00333     // inv_check.disable_printing_results();
00334     q.add_quality_assessor( &inv_check, err );
00335     MSQ_ERRZERO( err );
00336     Settings settings;
00337     // bug?  should we pass in pmesh?
00338     Mesh* mesh_ptr                  = &mesh;
00339     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( mesh_ptr, domain );
00340     q.run_common( &mesh_and_domain, 0, &settings, err );
00341     MSQ_ERRZERO( err );
00342     const QualityAssessor::Assessor* inv_b = inv_check.get_results( &metric );
00343     int num_invalid                        = inv_b->get_invalid_element_count();
00344     return num_invalid;
00345 }
00346 
00347 void ParShapeImprover::run( Mesh& mesh, MeshDomain* domain, MsqError&, bool always_smooth, int debug )
00348 {
00349     int rank, nprocs;
00350     MPI_Comm_rank( MPI_COMM_WORLD, &rank );
00351     MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
00352 
00353     if( debug )
00354     {
00355         MsqDebug::enable( 1 );
00356         if( debug > 1 ) MsqDebug::enable( 2 );
00357         if( debug > 2 ) MsqDebug::enable( 3 );
00358     }
00359 
00360     ParallelMesh* pmesh = dynamic_cast< ParallelMesh* >( &mesh );
00361     std::cout << "P[" << rank << "] "
00362               << " ParShapeImprover::run: pmesh= " << pmesh << std::endl;
00363 
00364     MsqError mErr;
00365     int num_invalid    = 0;
00366     bool check_quality = true;
00367     if( check_quality )
00368     {
00369         num_invalid = count_invalid_elements( mesh, domain );
00370         std::cout << "\nP[" << rank << "] "
00371                   << " ParShapeImprover num_invalid before= " << num_invalid
00372                   << ( num_invalid
00373                            ? " WARNING: invalid elements exist before MBMesquite smoothing"
00374                            : ( !always_smooth ? "WARNING: no smoothing requested since always_smooth=false" : " " ) )
00375                   << std::endl;
00376     }
00377 
00378     if( num_invalid || always_smooth )
00379     {
00380         bool use_canned_wrapper = false;
00381         if( use_canned_wrapper )
00382         {
00383             ShapeImprovementWrapper siw( mErr );
00384             if( pmesh )
00385                 siw.run_instructions( pmesh, domain, mErr );
00386             else
00387             {
00388                 MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &mesh, domain );
00389                 siw.run_instructions( &mesh_and_domain, mErr );
00390             }
00391         }
00392         else
00393         {
00394             // int  msq_debug             = debug; // 1,2,3 for more debug info
00395             // bool always_smooth_local   = false;
00396 
00397             bool do_untangle_only = false;
00398             ParShapeImprover::ParShapeImprovementWrapper siw( innerIter, 0.0, gradNorm, 100 );
00399             siw.m_do_untangle_only = do_untangle_only;
00400             if( pmesh )
00401                 siw.run_instructions( pmesh, domain, mErr );
00402             else
00403             {
00404                 MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &mesh, domain );
00405                 siw.run_instructions( &mesh_and_domain, mErr );
00406             }
00407         }
00408 
00409         std::cout << "\nP[" << rank << "] "
00410                   << " ParShapeImprover: MsqError after ShapeImprovementWrapper: " << mErr << std::endl;
00411 
00412         if( check_quality )
00413         {
00414             num_invalid = count_invalid_elements( mesh, domain );
00415             std::cout << "\nP[" << rank << "] "
00416                       << " ParShapeImprover num_invalid after= " << num_invalid << " "
00417                       << ( num_invalid ? " ERROR still have invalid elements after MBMesquite smoothing"
00418                                        : " SUCCESS: smoothed and removed invalid elements " )
00419                       << std::endl;
00420         }
00421 
00422         MSQ_ERRRTN( mErr );
00423     }
00424 }
00425 
00426 static int test( std::string filename_prefix, std::string mesh_topology_name, MeshDomain* domain = 0 )
00427 {
00428     int rank, nprocs;
00429     MPI_Comm_rank( MPI_COMM_WORLD, &rank );
00430     MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
00431 
00432     if( nprocs > 2 )
00433     {
00434         cerr << "parallel_untangle_shape::test(" << mesh_topology_name << " can only be run with 1 or 2 processors"
00435              << std::endl;
00436         return 0;
00437     }
00438 
00439     /* create processor-specific file names */
00440     ostringstream in_name, out_name, gold_name;
00441     in_name << filename_prefix << "par_untangle_original_" << mesh_topology_name << "_mesh." << nprocs << "." << rank
00442             << ".vtk";
00443     gold_name << filename_prefix << "par_untangle_smoothed_" << mesh_topology_name << "_mesh." << nprocs << "." << rank
00444               << ".vtk";
00445     out_name << "par_untangle_smoothed_" << mesh_topology_name << "_mesh." << nprocs << "." << rank << ".vtk";
00446 
00447     cout << "in_name= " << in_name.str() << " gold_name= " << gold_name.str() << " out_name= " << out_name.str()
00448          << std::endl;
00449 
00450     /* load different mesh files on each processor */
00451     MsqError err;
00452     MeshImpl mesh;
00453     mesh.read_vtk( in_name.str().c_str(), err );
00454     if( err )
00455     {
00456         cerr << err << endl;
00457         return 1;
00458     }
00459 
00460     /* create parallel mesh instance, specifying tags
00461      * containing parallel data */
00462     ParallelMeshImpl parallel_mesh( &mesh, "GLOBAL_ID", "PROCESSOR_ID" );
00463     ParallelHelperImpl helper;
00464     helper.set_communicator( MPI_COMM_WORLD );
00465     helper.set_parallel_mesh( &parallel_mesh );
00466     parallel_mesh.set_parallel_helper( &helper );
00467 
00468     /* do Laplacian smooth */
00469     // LaplaceWrapper optimizer;
00470     // optimizer.run_instructions(&parallel_mesh, err);
00471 
00472     int msq_debug      = 0;  // 1,2,3 for more debug info
00473     bool always_smooth = true;
00474     int innerIter      = 100;
00475     double gradNorm    = 1.e-9;
00476 
00477     ParShapeImprover si( innerIter, gradNorm );
00478     // Mesh *pmesh = &parallel_mesh;
00479     si.run( parallel_mesh, domain, err, always_smooth, msq_debug );
00480     if( err )
00481     {
00482         cerr << err << endl;
00483         return 1;
00484     }
00485 
00486     /* write mesh */
00487     mesh.write_vtk( out_name.str().c_str(), err );
00488     if( err )
00489     {
00490         cerr << err << endl;
00491         return 1;
00492     }
00493 
00494     // std::cout << "P[ " << rank <<"] reading gold..." << std::endl;
00495 
00496     /* compare mesh with gold copy */
00497     MeshImpl gold;
00498     gold.read_vtk( gold_name.str().c_str(), err );
00499     if( err )
00500     {
00501         cerr << err << endl;
00502         return 1;
00503     }
00504 
00505     // std::cout << "P[ " << rank <<"] read gold, checking mesh diff..." << std::endl;
00506     bool do_print = true;
00507     double tol    = 1.e-4;
00508     bool diff     = MeshUtil::meshes_are_different( mesh, gold, err, tol, do_print );
00509     if( err )
00510     {
00511         cerr << err << endl;
00512         return 1;
00513     }
00514     // std::cout << "P[ " << rank <<"] read gold, checking mesh diff...done" << std::endl;
00515 
00516     if( diff )
00517     {
00518         cerr << "Error, computed mesh different from gold copy" << std::endl;
00519         return 1;
00520     }
00521 
00522     print_timing_diagnostics( cout );
00523 
00524     return 0;
00525 }
00526 
00527 int main( int argc, char* argv[] )
00528 {
00529     /* init MPI */
00530     if( MPI_SUCCESS != MPI_Init( &argc, &argv ) )
00531     {
00532         cerr << "MPI_Init failed." << endl;
00533         return 2;
00534     }
00535 
00536     Vector3D pnt( 0, 0, 0 );
00537     Vector3D s_norm( 0, 0, 1 );
00538     PlanarDomain msq_geom( s_norm, pnt );
00539 
00540     int t1 = test( VTK_2D_DIR, "quad", &msq_geom );
00541     if( t1 ) return t1;
00542     t1 = test( VTK_3D_DIR, "hex" );
00543     if( t1 ) return t1;
00544 
00545     MPI_Finalize();
00546     return 0;
00547 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines