MOAB: Mesh Oriented datABase  (version 5.2.1)
VertexMover.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 /*!
00028   \file   VertexMover.cpp
00029   \brief
00030 
00031   The VertexMover Class is the base class for all the smoothing algorythms
00032 
00033   \author Thomas Leurent
00034   \date   2002-01-17
00035 */
00036 
00037 #include "VertexMover.hpp"
00038 #include "NonGradient.hpp"
00039 #include "QualityMetric.hpp"
00040 #include "TQualityMetric.hpp"
00041 #include "TMetricBarrier.hpp"
00042 #include "AWMetricBarrier.hpp"
00043 #include "ElementMaxQM.hpp"
00044 #include "ElementAvgQM.hpp"
00045 #include "ElementPMeanP.hpp"
00046 #include "PlanarDomain.hpp"
00047 #include "ObjectiveFunctionTemplate.hpp"
00048 #include "MaxTemplate.hpp"
00049 #include "MsqTimer.hpp"
00050 #include "MsqDebug.hpp"
00051 #include "PatchSet.hpp"
00052 #include "PatchData.hpp"
00053 #include "ParallelHelperInterface.hpp"
00054 #include <algorithm>
00055 
00056 namespace MBMesquite
00057 {
00058 
00059 extern int get_parallel_rank();
00060 extern int get_parallel_size();
00061 extern void parallel_barrier();
00062 
00063 VertexMover::VertexMover( ObjectiveFunction* OF ) : QualityImprover(), objFuncEval( OF ), jacobiOpt( false ) {}
00064 
00065 VertexMover::~VertexMover() {}
00066 
00067 /*
00068 
00069     +-----------+
00070     |Reset Outer|
00071     |Criterion  |
00072     +-----------+
00073           |
00074           V
00075           +
00076         /   \
00077        /Outer\  YES
00078 +--> <Criterion>-----> DONE
00079 |      \Done?/
00080 |       \   /
00081 |         +
00082 |         |NO
00083 |         V
00084 |   +-----------+
00085 1   |Reset Mesh |
00086 |   | Iteration |
00087 |   +-----------+
00088 |         |
00089 |         V
00090 |         +
00091 |       /  \
00092 |   NO /Next\
00093 +-----<Patch > <-----+
00094        \    /        |
00095         \  /         |
00096           +          |
00097        YES|          |
00098           V          |
00099     +-----------+    |
00100     |Reset Inner|    |
00101     |Criterion  |    2
00102     +-----------+    |
00103           |          |
00104           V          |
00105           +          |
00106         /   \        |
00107        /Inner\  YES  |
00108 +--> <Criterion>-----+    --------------+
00109 |      \Done?/                          |
00110 |       \   /                           |
00111 |         +                             |
00112 |         |NO                           |
00113 |         V                          Inside
00114 3   +-----------+                    Smoother
00115 |   |   Smooth  |                       |
00116 |   |   Patch   |                       |
00117 |   +-----------+                       |
00118 |         |                             |
00119 ----------+               --------------+
00120 
00121 */
00122 
00123 /*! \brief Improves the quality of the MeshSet, calling some
00124     methods specified in a class derived from VertexMover
00125 
00126     \param const MeshSet &: this MeshSet is looped over. Only the
00127     mutable data members are changed (such as currentVertexInd).
00128   */
00129 double VertexMover::loop_over_mesh( MeshDomainAssoc* mesh_and_domain, const Settings* settings, MsqError& err )
00130 {
00131     Mesh* mesh         = mesh_and_domain->get_mesh();
00132     MeshDomain* domain = mesh_and_domain->get_domain();
00133 
00134     TagHandle coord_tag      = 0;  // store uncommitted coords for jacobi optimization
00135     TagHandle* coord_tag_ptr = 0;
00136 
00137     TerminationCriterion* outer_crit = 0;
00138     TerminationCriterion* inner_crit = 0;
00139 
00140     // Clear culling flag, set hard fixed flag, etc on all vertices
00141     initialize_vertex_byte( mesh_and_domain, settings, err );
00142     MSQ_ERRZERO( err );
00143 
00144     // Get the patch data to use for the first iteration
00145     OFEvaluator& obj_func = get_objective_function_evaluator();
00146 
00147     // Check for impromer use of MaxTemplate
00148     ObjectiveFunctionTemplate* maxt_ptr = dynamic_cast< MaxTemplate* >( obj_func.get_objective_function() );
00149     if( maxt_ptr )
00150     {
00151         QualityImprover* ngqi_ptr = dynamic_cast< NonGradient* >( this );
00152         if( !ngqi_ptr )
00153             std::cout << "Warning: MaxTemplate results in non-differentiable objective function." << std::endl
00154                       << "   Therefore, it is best to use the NonGradient solver. Other Mesquite" << std::endl
00155                       << "   solvers require derivative information." << std::endl;
00156     }
00157 
00158     PatchData patch;
00159     patch.set_mesh( mesh );
00160     patch.set_domain( domain );
00161     if( settings ) patch.attach_settings( settings );
00162     bool one_patch = false, inner_crit_terminated, all_culled;
00163     std::vector< Mesh::VertexHandle > patch_vertices;
00164     std::vector< Mesh::ElementHandle > patch_elements;
00165     bool valid;
00166 
00167     PatchSet* patch_set = get_patch_set();
00168     if( !patch_set )
00169     {
00170         MSQ_SETERR( err )( "No PatchSet for QualityImprover!", MsqError::INVALID_STATE );
00171         return 0.0;
00172     }
00173     patch_set->set_mesh( mesh );
00174 
00175     std::vector< PatchSet::PatchHandle > patch_list;
00176     patch_set->get_patch_handles( patch_list, err );
00177     MSQ_ERRZERO( err );
00178 
00179     // check for inverted elements when using a barrietr target metric
00180     TQualityMetric* tqm_ptr        = NULL;
00181     TMetricBarrier* tm_ptr         = NULL;
00182     AWMetricBarrier* awm_ptr       = NULL;
00183     ElemSampleQM* sample_qm_ptr    = NULL;
00184     QualityMetric* qm_ptr          = NULL;
00185     TQualityMetric* pmeanp_ptr     = NULL;
00186     ElementMaxQM* elem_max_ptr     = NULL;
00187     ElementAvgQM* elem_avg_ptr     = NULL;
00188     ElementPMeanP* elem_pmeanp_ptr = NULL;
00189 
00190     ObjectiveFunctionTemplate* of_ptr = dynamic_cast< ObjectiveFunctionTemplate* >( obj_func.get_objective_function() );
00191     if( of_ptr ) qm_ptr = of_ptr->get_quality_metric();
00192     if( qm_ptr )
00193     {
00194         pmeanp_ptr      = dynamic_cast< TQualityMetric* >( qm_ptr );  // PMeanP case
00195         elem_max_ptr    = dynamic_cast< ElementMaxQM* >( qm_ptr );
00196         elem_avg_ptr    = dynamic_cast< ElementAvgQM* >( qm_ptr );
00197         elem_pmeanp_ptr = dynamic_cast< ElementPMeanP* >( qm_ptr );
00198     }
00199     if( elem_max_ptr )
00200         sample_qm_ptr = elem_max_ptr->get_quality_metric();
00201     else if( elem_pmeanp_ptr )
00202         sample_qm_ptr = elem_pmeanp_ptr->get_quality_metric();
00203     else if( elem_avg_ptr )
00204         sample_qm_ptr = elem_avg_ptr->get_quality_metric();
00205     else if( pmeanp_ptr )
00206     {
00207         tm_ptr  = dynamic_cast< TMetricBarrier* >( pmeanp_ptr->get_target_metric() );
00208         awm_ptr = dynamic_cast< AWMetricBarrier* >( pmeanp_ptr->get_target_metric() );
00209     }
00210 
00211     if( sample_qm_ptr || pmeanp_ptr )
00212     {
00213         if( !pmeanp_ptr )
00214         {
00215             tqm_ptr = dynamic_cast< TQualityMetric* >( sample_qm_ptr );
00216             if( tqm_ptr )
00217             {
00218                 tm_ptr  = dynamic_cast< TMetricBarrier* >( tqm_ptr->get_target_metric() );
00219                 awm_ptr = dynamic_cast< AWMetricBarrier* >( tqm_ptr->get_target_metric() );
00220             }
00221         }
00222         else
00223         {
00224             tqm_ptr = dynamic_cast< TQualityMetric* >( pmeanp_ptr );
00225         }
00226 
00227         if( tqm_ptr && ( tm_ptr || awm_ptr ) )
00228         {
00229             // check for inverted elements
00230             this->initialize( patch, err );
00231             std::vector< size_t > handles;
00232 
00233             // set up patch data
00234             std::vector< PatchSet::PatchHandle >::iterator patch_iter = patch_list.begin();
00235             while( patch_iter != patch_list.end() )
00236             {
00237                 do
00238                 {
00239                     patch_set->get_patch( *patch_iter, patch_elements, patch_vertices, err );
00240                     if( MSQ_CHKERR( err ) ) goto ERROR;
00241                     ++patch_iter;
00242                 } while( patch_elements.empty() && patch_iter != patch_list.end() );
00243 
00244                 patch.set_mesh_entities( patch_elements, patch_vertices, err );
00245                 if( MSQ_CHKERR( err ) ) goto ERROR;
00246 
00247                 qm_ptr->get_evaluations( patch, handles, true, err );  // MSQ_ERRFALSE(err);
00248 
00249                 // do actual check for inverted elements
00250                 std::vector< size_t >::const_iterator i;
00251                 double tvalue;
00252                 for( i = handles.begin(); i != handles.end(); ++i )
00253                 {
00254                     bool result = tqm_ptr->evaluate( patch, *i, tvalue, err );
00255                     if( MSQ_CHKERR( err ) || !result ) return false;  // inverted element detected
00256                 }
00257             }
00258         }
00259     }
00260 
00261     // Get termination criteria
00262     outer_crit = this->get_outer_termination_criterion();
00263     inner_crit = this->get_inner_termination_criterion();
00264     if( outer_crit == 0 )
00265     {
00266         MSQ_SETERR( err )( "Termination Criterion pointer is Null", MsqError::INVALID_STATE );
00267         return 0.;
00268     }
00269     if( inner_crit == 0 )
00270     {
00271         MSQ_SETERR( err )
00272         ( "Termination Criterion pointer for inner loop is Null", MsqError::INVALID_STATE );
00273         return 0.;
00274     }
00275 
00276     // Set Termination Criterion defaults if no other Criterion is set
00277     if( !outer_crit->criterion_is_set() ) outer_crit->add_iteration_limit( 1 );
00278     if( !inner_crit->criterion_is_set() ) inner_crit->add_iteration_limit( 10 );
00279 
00280     // If using a local patch, suppress output of inner termination criterion
00281     if( patch_list.size() > 1 )
00282         inner_crit->set_debug_output_level( 3 );
00283     else
00284         one_patch = true;
00285 
00286     if( jacobiOpt )
00287     {
00288         coord_tag = get_jacobi_coord_tag( mesh, err );
00289         MSQ_ERRZERO( err );
00290         coord_tag_ptr = &coord_tag;
00291     }
00292 
00293     // Initialize outer loop
00294 
00295     this->initialize( patch, err );
00296     if( MSQ_CHKERR( err ) ) goto ERROR;
00297 
00298     valid = obj_func.initialize( mesh_and_domain, settings, patch_set, err );
00299     if( MSQ_CHKERR( err ) ) goto ERROR;
00300     if( !valid )
00301     {
00302         MSQ_SETERR( err )
00303         ( "ObjectiveFunction initialization failed.  Mesh "
00304           "invalid at one or more sample points.",
00305           MsqError::INVALID_MESH );
00306         goto ERROR;
00307     }
00308 
00309     outer_crit->reset_outer( mesh, domain, obj_func, settings, err );
00310     if( MSQ_CHKERR( err ) ) goto ERROR;
00311 
00312     // if only one patch, get the patch now
00313     if( one_patch )
00314     {
00315         patch_set->get_patch( patch_list[0], patch_elements, patch_vertices, err );
00316         if( MSQ_CHKERR( err ) ) goto ERROR;
00317         patch.set_mesh_entities( patch_elements, patch_vertices, err );
00318         if( MSQ_CHKERR( err ) ) goto ERROR;
00319     }
00320 
00321     // Loop until outer termination criterion is met
00322     inner_crit_terminated = false;
00323     while( !outer_crit->terminate() )
00324     {
00325         if( inner_crit_terminated )
00326         {
00327             MSQ_SETERR( err )
00328             ( "Inner termiation criterion satisfied for all patches "
00329               "without meeting outer termination criterion.  This is "
00330               "an infinite loop.  Aborting.",
00331               MsqError::INVALID_STATE );
00332             break;
00333         }
00334         inner_crit_terminated = true;
00335         all_culled            = true;
00336 
00337         int num_patches = 0;
00338         // Loop over each patch
00339         std::vector< PatchSet::PatchHandle >::iterator p_iter = patch_list.begin();
00340         while( p_iter != patch_list.end() )
00341         {
00342             if( !one_patch )
00343             {  // if only one patch (global) re-use the previous one
00344                 // loop until we get a non-empty patch.  patch will be empty
00345                 // for culled vertices with element-on-vertex patches
00346                 do
00347                 {
00348                     patch_set->get_patch( *p_iter, patch_elements, patch_vertices, err );
00349                     if( MSQ_CHKERR( err ) ) goto ERROR;
00350                     ++p_iter;
00351                 } while( patch_elements.empty() && p_iter != patch_list.end() );
00352 
00353                 if( patch_elements.empty() )
00354                 {  // no more non-culled vertices
00355                     if( 0 ) std::cout << "P[" << get_parallel_rank() << "] tmp srk all vertices culled." << std::endl;
00356                     break;
00357                 }
00358 
00359                 all_culled = false;
00360                 patch.set_mesh_entities( patch_elements, patch_vertices, err );
00361                 if( MSQ_CHKERR( err ) ) goto ERROR;
00362             }
00363             else
00364             {
00365                 ++p_iter;
00366                 all_culled = false;
00367             }
00368 
00369             ++num_patches;
00370 
00371             // Initialize for inner iteration
00372 
00373             this->initialize_mesh_iteration( patch, err );
00374             if( MSQ_CHKERR( err ) ) goto ERROR;
00375 
00376             obj_func.reset();
00377 
00378             outer_crit->reset_patch( patch, err );
00379             if( MSQ_CHKERR( err ) ) goto ERROR;
00380 
00381             inner_crit->reset_inner( patch, obj_func, err );
00382             if( MSQ_CHKERR( err ) ) goto ERROR;
00383 
00384             inner_crit->reset_patch( patch, err );
00385             if( MSQ_CHKERR( err ) ) goto ERROR;
00386 
00387             // Don't even call optimizer if inner termination
00388             // criterion has already been met.
00389             if( !inner_crit->terminate() )
00390             {
00391                 inner_crit_terminated = false;
00392 
00393                 // Call optimizer - should loop on inner_crit->terminate()
00394                 this->optimize_vertex_positions( patch, err );
00395                 if( MSQ_CHKERR( err ) ) goto ERROR;
00396 
00397                 // Update for changes during inner iteration
00398                 // (during optimizer loop)
00399 
00400                 outer_crit->accumulate_patch( patch, err );
00401                 if( MSQ_CHKERR( err ) ) goto ERROR;
00402 
00403                 inner_crit->cull_vertices( patch, obj_func, err );
00404                 if( MSQ_CHKERR( err ) ) goto ERROR;
00405 
00406                 // FIXME
00407                 if( 0 )
00408                 {
00409                     inner_crit->cull_vertices_global( patch, mesh, domain, settings, obj_func, err );
00410                     if( MSQ_CHKERR( err ) ) goto ERROR;
00411                 }
00412 
00413                 patch.update_mesh( err, coord_tag_ptr );
00414                 if( MSQ_CHKERR( err ) ) goto ERROR;
00415             }
00416         }
00417 
00418         if( jacobiOpt ) commit_jacobi_coords( coord_tag, mesh, err );
00419 
00420         this->terminate_mesh_iteration( patch, err );
00421         if( MSQ_CHKERR( err ) ) goto ERROR;
00422 
00423         outer_crit->accumulate_outer( mesh, domain, obj_func, settings, err );
00424         if( MSQ_CHKERR( err ) ) goto ERROR;
00425 
00426         if( all_culled ) break;
00427     }
00428 
00429 ERROR:
00430     if( jacobiOpt ) mesh->tag_delete( coord_tag, err );
00431 
00432     // call the criteria's cleanup funtions.
00433     if( outer_crit ) outer_crit->cleanup( mesh, domain, err );
00434     if( inner_crit ) inner_crit->cleanup( mesh, domain, err );
00435     // call the optimization cleanup function.
00436     this->cleanup();
00437 
00438     return 0.;
00439 }
00440 
00441 static void checkpoint_bytes( Mesh* mesh, std::vector< unsigned char >& saved_bytes, MsqError& err )
00442 {
00443     std::vector< Mesh::VertexHandle > vertexHandlesArray;
00444     mesh->get_all_vertices( vertexHandlesArray, err );MSQ_ERRRTN( err );
00445     saved_bytes.resize( vertexHandlesArray.size() );
00446     mesh->vertices_get_byte( &vertexHandlesArray[0], &saved_bytes[0], vertexHandlesArray.size(), err );MSQ_ERRRTN( err );
00447 }
00448 
00449 static void restore_bytes( Mesh* mesh, std::vector< unsigned char >& saved_bytes, MsqError& err )
00450 {
00451     std::vector< Mesh::VertexHandle > vertexHandlesArray;
00452     mesh->get_all_vertices( vertexHandlesArray, err );MSQ_ERRRTN( err );
00453     mesh->vertices_set_byte( &vertexHandlesArray[0], &saved_bytes[0], vertexHandlesArray.size(), err );MSQ_ERRRTN( err );
00454 }
00455 
00456 static void save_or_restore_debug_state( bool save )
00457 {
00458     static bool debug[3] = { false, false, false };
00459     if( save )
00460     {
00461         debug[0] = MsqDebug::get( 1 );
00462         debug[1] = MsqDebug::get( 2 );
00463         debug[2] = MsqDebug::get( 3 );
00464     }
00465     else
00466     {
00467         if( debug[0] ) MsqDebug::enable( 1 );
00468         if( debug[1] ) MsqDebug::enable( 2 );
00469         if( debug[2] ) MsqDebug::enable( 3 );
00470     }
00471 }
00472 
00473 /*! \brief Improves the quality of the MeshSet, calling some
00474     methods specified in a class derived from VertexMover
00475 
00476     \param const MeshSet &: this MeshSet is looped over. Only the
00477     mutable data members are changed (such as currentVertexInd).
00478   */
00479 double VertexMover::loop_over_mesh( ParallelMesh* mesh, MeshDomain* domain, const Settings* settings, MsqError& err )
00480 {
00481     std::vector< size_t > junk;
00482     Mesh::VertexHandle vertex_handle;
00483     TagHandle coord_tag      = 0;  // store uncommitted coords for jacobi optimization
00484     TagHandle* coord_tag_ptr = 0;
00485     int outer_iter           = 0;
00486     int inner_iter           = 0;
00487     bool one_patch           = false;
00488 
00489     // Clear culling flag, set hard fixed flag, etc on all vertices
00490     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( (Mesh*)mesh, 0 );
00491     initialize_vertex_byte( &mesh_and_domain, settings, err );
00492     MSQ_ERRZERO( err );
00493 
00494     // Get the patch data to use for the first iteration
00495     OFEvaluator& obj_func = get_objective_function_evaluator();
00496 
00497     PatchData patch;
00498     patch.set_mesh( (Mesh*)mesh );
00499     patch.set_domain( domain );
00500     patch.attach_settings( settings );
00501 
00502     ParallelHelper* helper = mesh->get_parallel_helper();
00503     if( !helper )
00504     {
00505         MSQ_SETERR( err )( "No ParallelHelper instance", MsqError::INVALID_STATE );
00506         return 0;
00507     }
00508 
00509     helper->smoothing_init( err );
00510     MSQ_ERRZERO( err );
00511 
00512     bool inner_crit_terminated, all_culled;
00513     std::vector< Mesh::VertexHandle > patch_vertices;
00514     std::vector< Mesh::ElementHandle > patch_elements;
00515     std::vector< Mesh::VertexHandle > fixed_vertices;
00516     std::vector< Mesh::VertexHandle > free_vertices;
00517 
00518     // Get termination criteria
00519     TerminationCriterion* outer_crit = this->get_outer_termination_criterion();
00520     TerminationCriterion* inner_crit = this->get_inner_termination_criterion();
00521     if( outer_crit == 0 )
00522     {
00523         MSQ_SETERR( err )( "Termination Criterion pointer is Null", MsqError::INVALID_STATE );
00524         return 0.;
00525     }
00526     if( inner_crit == 0 )
00527     {
00528         MSQ_SETERR( err )
00529         ( "Termination Criterion pointer for inner loop is Null", MsqError::INVALID_STATE );
00530         return 0.;
00531     }
00532 
00533     PatchSet* patch_set = get_patch_set();
00534     if( !patch_set )
00535     {
00536         MSQ_SETERR( err )( "No PatchSet for QualityImprover!", MsqError::INVALID_STATE );
00537         return 0.0;
00538     }
00539     patch_set->set_mesh( (Mesh*)mesh );
00540 
00541     std::vector< PatchSet::PatchHandle > patch_list;
00542     patch_set->get_patch_handles( patch_list, err );
00543     MSQ_ERRZERO( err );
00544 
00545     if( patch_list.size() > 1 )
00546         inner_crit->set_debug_output_level( 3 );
00547     else
00548         one_patch = true;
00549 
00550     if( jacobiOpt )
00551     {
00552         coord_tag = get_jacobi_coord_tag( mesh, err );
00553         MSQ_ERRZERO( err );
00554         coord_tag_ptr = &coord_tag;
00555     }
00556 
00557     // parallel error checking
00558     MsqError perr;
00559     bool pdone = false;
00560 
00561 #define PERROR_COND continue
00562 
00563     // Initialize outer loop
00564 
00565     this->initialize( patch, err );
00566     if( MSQ_CHKERR( err ) ) { MSQ_SETERR( perr )( "initialize patch", MsqError::INVALID_STATE ); }  // goto ERROR;
00567 
00568     MeshDomainAssoc mesh_and_domain2 = MeshDomainAssoc( (Mesh*)mesh, domain );
00569     obj_func.initialize( &mesh_and_domain2, settings, patch_set, err );
00570     if( MSQ_CHKERR( err ) ) { MSQ_SETERR( perr )( "initialize obj_func", MsqError::INVALID_STATE ); }  // goto ERROR;
00571 
00572     outer_crit->reset_outer( (Mesh*)mesh, domain, obj_func, settings, err );
00573     if( MSQ_CHKERR( err ) ) { MSQ_SETERR( perr )( "reset_outer", MsqError::INVALID_STATE ); }  // goto ERROR;
00574 
00575     // Loop until outer termination criterion is met
00576     inner_crit_terminated = false;
00577     all_culled            = false;
00578     for( ;; )
00579     {
00580         if( 0 )
00581             std::cout << "P[" << get_parallel_rank() << "] tmp srk inner_iter= " << inner_iter
00582                       << " outer_iter= " << outer_iter << std::endl;
00583 
00584         // PERROR:
00585 
00586         if( MSQ_CHKERR( perr ) )
00587         {
00588             std::cout << "P[" << get_parallel_rank() << "] VertexMover::loop_over_mesh found parallel error: " << perr
00589                       << "\n quitting... pdone= " << pdone << std::endl;
00590             pdone = true;
00591         }
00592 
00593         helper->communicate_any_true( pdone, err );
00594 
00595         if( 0 )
00596             std::cout << "P[" << get_parallel_rank() << "] tmp srk inner_iter= " << inner_iter
00597                       << " outer_iter= " << outer_iter << " pdone= " << pdone << std::endl;
00598 
00599         if( pdone )
00600         {
00601             std::cout << "P[" << get_parallel_rank()
00602                       << "] VertexMover::loop_over_mesh found parallel error, quitting... pdone= " << pdone
00603                       << std::endl;
00604             MSQ_SETERR( err )( "PARALLEL ERROR", MsqError::PARALLEL_ERROR );
00605             break;
00606         }
00607 
00608         ++outer_iter;
00609 
00610         /// srkenno@sandia.gov 1/19/12: the logic here was changed so that all proc's must agree
00611         ///   on the values used for outer and inner termination before the iteration is stopped.
00612         ///   Previously, the ParallelHelper::communicate_all_true method returned true if any
00613         ///   proc sent it a true value, which seems to be a bug, at least in the name of the
00614         ///   method. The method has been changed to return true only if all proc's values are true.
00615         ///   In the previous version, this meant that if any proc hit its inner or outer
00616         ///   termination criterion, the loop was exited, and thus some parts of the mesh
00617         ///   are potentially left unconverged.  Also, if the outer criterion was satisfied on
00618         ///   part of the mesh (say a uniform part), the iterations were not executed at all.
00619         /// Also, changed name of "did_some" to "inner_crit_terminated", and flipped its boolean
00620         ///   value to be consistent with the name - for readability and for correctness since
00621         ///   we want to communicate a true value through the helper.
00622 
00623         bool outer_crit_terminated       = outer_crit->terminate();
00624         bool outer_crit_terminated_local = outer_crit_terminated;
00625         helper->communicate_all_true( outer_crit_terminated, err );
00626 
00627         bool inner_crit_terminated_local = inner_crit_terminated;
00628         helper->communicate_all_true( inner_crit_terminated, err );
00629 
00630         bool all_culled_local = all_culled;
00631         helper->communicate_all_true( all_culled, err );
00632 
00633         bool done = all_culled || outer_crit_terminated;
00634         if( inner_crit_terminated )
00635         {
00636             MSQ_SETERR( err )
00637             ( "Inner termination criterion satisfied for all patches "
00638               "without meeting outer termination criterion.  This is "
00639               "an infinite loop.  Aborting.",
00640               MsqError::INVALID_STATE );
00641             done = true;
00642             helper->communicate_any_true( done, err );
00643         }
00644 
00645         bool local_done = done;
00646 
00647         helper->communicate_all_true( done, err );
00648 
00649         if( 0 )
00650             std::cout << "P[" << get_parallel_rank() << "] tmp srk done= " << done << " local_done= " << local_done
00651                       << " all_culled= " << all_culled << " outer_crit->terminate()= " << outer_crit->terminate()
00652                       << " outer_term= " << outer_crit_terminated
00653                       << " outer_term_local= " << outer_crit_terminated_local
00654                       << " inner_crit_terminated = " << inner_crit_terminated
00655                       << " inner_crit_terminated_local = " << inner_crit_terminated_local
00656                       << " all_culled = " << all_culled << " all_culled_local = " << all_culled_local << std::endl;
00657 
00658         if( MSQ_CHKERR( err ) ) { MSQ_SETERR( perr )( "loop start", MsqError::INVALID_STATE ); }  // goto ERROR;
00659 
00660         if( done ) break;
00661 
00662         inner_crit_terminated = true;
00663         all_culled            = true;
00664 
00665         ///*** smooth the interior ***////
00666 
00667         // get the fixed vertices (i.e. the ones *not* part of the first independent set)
00668         helper->compute_first_independent_set( fixed_vertices );
00669 
00670         // sort the fixed vertices
00671         std::sort( fixed_vertices.begin(), fixed_vertices.end() );
00672 
00673         // Loop over each patch
00674         if( 0 && MSQ_DBG( 2 ) )
00675             std::cout << "P[" << get_parallel_rank() << "] tmp srk number of patches = " << patch_list.size()
00676                       << " inner_iter= " << inner_iter << " outer_iter= " << outer_iter
00677                       << " inner.globalInvertedCount = " << inner_crit->globalInvertedCount
00678                       << " outer.globalInvertedCount = " << outer_crit->globalInvertedCount
00679                       << " inner.patchInvertedCount = " << inner_crit->patchInvertedCount
00680                       << " outer.patchInvertedCount = " << outer_crit->patchInvertedCount << std::endl;
00681 
00682         save_or_restore_debug_state( true );
00683         // MsqDebug::disable_all();
00684 
00685         int num_patches = 0;
00686 
00687         std::vector< PatchSet::PatchHandle >::iterator p_iter = patch_list.begin();
00688         while( p_iter != patch_list.end() )
00689         {
00690 
00691             // loop until we get a non-empty patch.  patch will be empty
00692             // for culled vertices with element-on-vertex patches
00693             do
00694             {
00695                 patch_set->get_patch( *p_iter, patch_elements, patch_vertices, err );
00696                 if( MSQ_CHKERR( err ) )
00697                 {
00698                     MSQ_SETERR( perr )( "get_patch", MsqError::INVALID_STATE );
00699                     PERROR_COND;
00700                 }  // goto ERROR;
00701                 ++p_iter;
00702             } while( patch_elements.empty() && p_iter != patch_list.end() );
00703 
00704             if( patch_elements.empty() )
00705             {  // no more non-culled vertices
00706                 if( 0 ) std::cout << "P[" << get_parallel_rank() << "] tmp srk all vertices culled." << std::endl;
00707                 break;
00708             }
00709 
00710             if( patch_vertices.empty() )  // global patch hack (means all mesh vertices)
00711             { mesh->get_all_vertices( patch_vertices, err ); }
00712 
00713             free_vertices.clear();
00714 
00715             for( size_t i = 0; i < patch_vertices.size(); ++i )
00716                 if( !std::binary_search( fixed_vertices.begin(), fixed_vertices.end(), patch_vertices[i] ) )
00717                     free_vertices.push_back( patch_vertices[i] );
00718 
00719             if( free_vertices.empty() )
00720             {  // all vertices were fixed -> skip patch
00721                 continue;
00722             }
00723 
00724             ++num_patches;
00725             all_culled = false;
00726             patch.set_mesh_entities( patch_elements, free_vertices, err );
00727             if( MSQ_CHKERR( err ) )
00728             {
00729                 MSQ_SETERR( perr )( "set_mesh_entities", MsqError::INVALID_STATE );
00730                 PERROR_COND;
00731             }  // goto ERROR;
00732 
00733             // Initialize for inner iteration
00734 
00735             this->initialize_mesh_iteration( patch, err );
00736             if( MSQ_CHKERR( err ) )
00737             {
00738                 MSQ_SETERR( perr )( "initialize_mesh_iteration", MsqError::INVALID_STATE );
00739                 PERROR_COND;
00740             }  // goto ERROR;
00741 
00742             obj_func.reset();
00743 
00744             outer_crit->reset_patch( patch, err );
00745             if( MSQ_CHKERR( err ) )
00746             {
00747                 MSQ_SETERR( perr )( "reset_patch outer", MsqError::INVALID_STATE );
00748                 PERROR_COND;
00749             }  // goto ERROR;
00750 
00751             inner_crit->reset_inner( patch, obj_func, err );
00752             if( MSQ_CHKERR( err ) )
00753             {
00754                 MSQ_SETERR( perr )( "reset_inner", MsqError::INVALID_STATE );
00755                 PERROR_COND;
00756             }  // goto ERROR;
00757 
00758             inner_crit->reset_patch( patch, err );
00759             if( MSQ_CHKERR( err ) )
00760             {
00761                 MSQ_SETERR( perr )( "inner reset_patch", MsqError::INVALID_STATE );
00762                 PERROR_COND;
00763             }  // goto ERROR;
00764 
00765             // Don't even call optimizer if inner termination
00766             // criterion has already been met.
00767             if( !inner_crit->terminate() )
00768             {
00769                 inner_crit_terminated = false;
00770                 if( one_patch ) ++inner_iter;
00771 
00772                 // Call optimizer - should loop on inner_crit->terminate()
00773                 // size_t num_vert=patch.num_free_vertices();
00774                 // std::cout << "P[" << get_parallel_rank() << "] tmp srk VertexMover num_vert= " <<
00775                 // num_vert << std::endl;
00776 
00777                 this->optimize_vertex_positions( patch, err );
00778                 if( MSQ_CHKERR( err ) )
00779                 {
00780                     MSQ_SETERR( perr )( "optimize_vertex_positions", MsqError::INVALID_STATE );
00781                     PERROR_COND;
00782                 }  // goto ERROR;
00783 
00784                 // Update for changes during inner iteration
00785                 // (during optimizer loop)
00786 
00787                 outer_crit->accumulate_patch( patch, err );
00788                 if( MSQ_CHKERR( err ) )
00789                 {
00790                     MSQ_SETERR( perr )( "outer accumulate_patch", MsqError::INVALID_STATE );
00791                     PERROR_COND;
00792                 }  // goto ERROR;
00793 
00794                 inner_crit->cull_vertices( patch, obj_func, err );
00795                 if( MSQ_CHKERR( err ) )
00796                 {
00797                     MSQ_SETERR( perr )( "inner cull_vertices", MsqError::INVALID_STATE );
00798                     PERROR_COND;
00799                 }  // goto ERROR;
00800 
00801                 // experimental...
00802                 if( 0 )
00803                 {
00804                     inner_crit->cull_vertices_global( patch, mesh, domain, settings, obj_func, err );
00805                     if( MSQ_CHKERR( err ) )
00806                     {
00807                         MSQ_SETERR( perr )( "cull_vertices_global", MsqError::INVALID_STATE );
00808                         PERROR_COND;
00809                     }  // goto ERROR;
00810                 }
00811 
00812                 patch.update_mesh( err, coord_tag_ptr );
00813                 if( MSQ_CHKERR( err ) )
00814                 {
00815                     MSQ_SETERR( perr )( "update_mesh", MsqError::INVALID_STATE );
00816                     PERROR_COND;
00817                 }  // goto ERROR;
00818             }
00819         }  // while(p_iter....
00820 
00821         save_or_restore_debug_state( false );
00822 
00823         if( 1 )
00824         {
00825             bool pdone_inner = false;
00826 
00827             if( MSQ_CHKERR( perr ) )
00828             {
00829                 std::cout << "P[" << get_parallel_rank()
00830                           << "] VertexMover::loop_over_mesh found parallel error: " << perr
00831                           << "\n quitting... pdone_inner= " << pdone_inner << std::endl;
00832                 pdone_inner = true;
00833             }
00834 
00835             helper->communicate_any_true( pdone_inner, err );
00836 
00837             if( 0 )
00838                 std::cout << "P[" << get_parallel_rank() << "] tmp srk inner_iter= " << inner_iter
00839                           << " outer_iter= " << outer_iter << " pdone_inner= " << pdone_inner << std::endl;
00840 
00841             if( pdone_inner )
00842             {
00843                 if( 0 )
00844                     std::cout << "P[" << get_parallel_rank()
00845                               << "] tmp srk found parallel error, quitting... pdone_inner= " << pdone_inner
00846                               << std::endl;
00847                 MSQ_SETERR( err )( "PARALLEL ERROR", MsqError::PARALLEL_ERROR );
00848                 break;
00849             }
00850         }
00851 
00852         /// srkenno@sandia.gov save vertex bytes since boundary smoothing changes them
00853         std::vector< unsigned char > saved_bytes;
00854         checkpoint_bytes( mesh, saved_bytes, err );
00855         if( MSQ_CHKERR( err ) )
00856         {
00857             MSQ_SETERR( perr )( "checkpoint_bytes ", MsqError::INVALID_STATE );
00858             PERROR_COND;
00859         }  // goto ERROR;
00860 
00861         helper->communicate_first_independent_set( err );
00862         if( MSQ_CHKERR( err ) )
00863         {
00864             MSQ_SETERR( perr )( "communicate_first_independent_set ", MsqError::INVALID_STATE );
00865             PERROR_COND;
00866         }  // goto ERROR;
00867 
00868         ///*** smooth the boundary ***////
00869         save_or_restore_debug_state( true );
00870         MsqDebug::disable_all();
00871 
00872         while( helper->compute_next_independent_set() )
00873         {
00874             // Loop over all boundary elements
00875             while( helper->get_next_partition_boundary_vertex( vertex_handle ) )
00876             {
00877 
00878                 patch_vertices.clear();
00879                 patch_vertices.push_back( vertex_handle );
00880                 patch_elements.clear();
00881                 mesh->vertices_get_attached_elements( &vertex_handle, 1, patch_elements, junk, err );
00882 
00883                 all_culled = false;
00884                 patch.set_mesh_entities( patch_elements, patch_vertices, err );
00885 
00886                 if( MSQ_CHKERR( err ) )
00887                 {
00888                     MSQ_SETERR( perr )( "set_mesh_entities 2 ", MsqError::INVALID_STATE );
00889                     PERROR_COND;
00890                 }  // goto ERROR;
00891 
00892                 // Initialize for inner iteration
00893 
00894                 this->initialize_mesh_iteration( patch, err );
00895                 if( MSQ_CHKERR( err ) )
00896                 {
00897                     MSQ_SETERR( perr )( " initialize_mesh_iteration 2", MsqError::INVALID_STATE );
00898                     PERROR_COND;
00899                 }  // goto ERROR;
00900 
00901                 obj_func.reset();
00902 
00903                 outer_crit->reset_patch( patch, err );
00904                 if( MSQ_CHKERR( err ) )
00905                 {
00906                     MSQ_SETERR( perr )( "outer reset_patch 2 ", MsqError::INVALID_STATE );
00907                     PERROR_COND;
00908                 }  // goto ERROR;
00909 
00910                 inner_crit->reset_inner( patch, obj_func, err );
00911                 if( MSQ_CHKERR( err ) )
00912                 {
00913                     MSQ_SETERR( perr )( " inner reset_inner 2", MsqError::INVALID_STATE );
00914                     PERROR_COND;
00915                 }  // goto ERROR;
00916 
00917                 inner_crit->reset_patch( patch, err );
00918                 if( MSQ_CHKERR( err ) )
00919                 {
00920                     MSQ_SETERR( perr )( " inner_crit reset_patch 2", MsqError::INVALID_STATE );
00921                     PERROR_COND;
00922                 }  // goto ERROR;
00923 
00924                 // Don't even call optimizer if inner termination
00925                 // criterion has already been met.
00926                 if( !inner_crit->terminate() )
00927                 {
00928                     inner_crit_terminated = false;
00929 
00930                     // Call optimizer - should loop on inner_crit->terminate()
00931                     this->optimize_vertex_positions( patch, err );
00932                     if( MSQ_CHKERR( err ) )
00933                     {
00934                         MSQ_SETERR( perr )
00935                         ( " optimize_vertex_positions 2", MsqError::INVALID_STATE );
00936                         PERROR_COND;
00937                     }  // goto ERROR;
00938 
00939                     // Update for changes during inner iteration
00940                     // (during optimizer loop)
00941 
00942                     outer_crit->accumulate_patch( patch, err );
00943                     if( MSQ_CHKERR( err ) )
00944                     {
00945                         MSQ_SETERR( perr )( " outer accumulate_patch 2", MsqError::INVALID_STATE );
00946                         PERROR_COND;
00947                     }  // goto ERROR;
00948 
00949                     inner_crit->cull_vertices( patch, obj_func, err );
00950                     if( MSQ_CHKERR( err ) )
00951                     {
00952                         MSQ_SETERR( perr )( " inner cull_vertices", MsqError::INVALID_STATE );
00953                         PERROR_COND;
00954                     }  // goto ERROR;
00955 
00956                     // FIXME
00957                     if( 0 )
00958                     {
00959                         inner_crit->cull_vertices_global( patch, mesh, domain, settings, obj_func, err );
00960                         if( MSQ_CHKERR( err ) )
00961                         {
00962                             MSQ_SETERR( perr )
00963                             ( " cull_vertices_global 2 ", MsqError::INVALID_STATE );
00964                             PERROR_COND;
00965                         }  // goto ERROR;
00966                     }
00967 
00968                     patch.update_mesh( err, coord_tag_ptr );
00969                     if( MSQ_CHKERR( err ) )
00970                     {
00971                         MSQ_SETERR( perr )( " update_mesh 2", MsqError::INVALID_STATE );
00972                         PERROR_COND;
00973                     }  // goto ERROR;
00974                 }
00975             }
00976             helper->communicate_next_independent_set( err );
00977             if( MSQ_CHKERR( err ) )
00978             {
00979                 MSQ_SETERR( perr )
00980                 ( " communicate_next_independent_set 2", MsqError::INVALID_STATE );
00981                 PERROR_COND;
00982             }  // goto ERROR;
00983         }      // while(helper->compute_next_independent_set())
00984 
00985         save_or_restore_debug_state( false );
00986 
00987         // if (!get_parallel_rank())
00988         // std::cout << "P[" << get_parallel_rank() << "] tmp srk num_patches= " << num_patches <<
00989         // std::endl;
00990 
00991         /// srkenno@sandia.gov restore vertex bytes since boundary smoothing changes them
00992         restore_bytes( mesh, saved_bytes, err );
00993         if( MSQ_CHKERR( err ) )
00994         {
00995             MSQ_SETERR( perr )( " restore_bytes", MsqError::INVALID_STATE );
00996             PERROR_COND;
00997         }  // goto ERROR;
00998 
00999         if( jacobiOpt ) commit_jacobi_coords( coord_tag, mesh, err );
01000 
01001         this->terminate_mesh_iteration( patch, err );
01002         if( MSQ_CHKERR( err ) )
01003         {
01004             MSQ_SETERR( perr )( " terminate_mesh_iteration", MsqError::INVALID_STATE );
01005             PERROR_COND;
01006         }  // goto ERROR;
01007 
01008         outer_crit->accumulate_outer( mesh, domain, obj_func, settings, err );
01009         if( MSQ_CHKERR( err ) )
01010         {
01011             MSQ_SETERR( perr )( " outer_crit accumulate_outer", MsqError::INVALID_STATE );
01012             PERROR_COND;
01013         }  // goto ERROR;
01014     }
01015 
01016     // ERROR:
01017 
01018     if( MSQ_CHKERR( err ) )
01019     {
01020         std::cout << "P[" << get_parallel_rank() << "] VertexMover::loop_over_mesh error = " << err.error_message()
01021                   << std::endl;
01022     }
01023 
01024     if( jacobiOpt ) mesh->tag_delete( coord_tag, err );
01025 
01026     // call the criteria's cleanup funtions.
01027     outer_crit->cleanup( mesh, domain, err );MSQ_CHKERR( err );
01028     inner_crit->cleanup( mesh, domain, err );MSQ_CHKERR( err );
01029     // call the optimization cleanup function.
01030     this->cleanup();
01031     // close the helper
01032     helper->smoothing_close( err );MSQ_CHKERR( err );
01033 
01034     return 0.;
01035 }
01036 
01037 void VertexMover::initialize_queue( MeshDomainAssoc* mesh_and_domain, const Settings* settings, MsqError& err )
01038 {
01039     QualityImprover::initialize_queue( mesh_and_domain, settings, err );MSQ_ERRRTN( err );
01040     objFuncEval.initialize_queue( mesh_and_domain, settings, err );MSQ_ERRRTN( err );
01041 }
01042 
01043 TagHandle VertexMover::get_jacobi_coord_tag( Mesh* mesh, MsqError& err )
01044 {
01045     // Get tag handle
01046     const char tagname[] = "msq_jacobi_temp_coords";
01047     TagHandle tag        = mesh->tag_create( tagname, Mesh::DOUBLE, 3, 0, err );
01048     MSQ_ERRZERO( err );
01049     /* VertexMover will always delete the tag when it is done, so it is probably
01050        best not to accept an existing tag
01051   if (err.error_code() == TAG_ALREADY_EXISTS) {
01052     err.clear();
01053     tag = tag_get( tagname, err ); MSQ_ERRZERO(err);
01054     std::string name;
01055     Mesh::TagType type;
01056     unsigned length;
01057     mesh->tag_properties( tag, name, type, length, err ); MSQ_ERRZERO(err);
01058     if (type != Mesh::DOUBLE || length != 3) {
01059       MSQ_SETERR(err)(TAG_ALREADY_EXISTS,
01060                      "Tag \"%s\" already exists with invalid type",
01061                      tagname);
01062       return 0;
01063     }
01064   }
01065     */
01066 
01067     // Initialize tag value with initial vertex coordinates so that
01068     // vertices that never get moved end up with the correct coords.
01069     std::vector< Mesh::VertexHandle > vertices;
01070     mesh->get_all_vertices( vertices, err );
01071     MSQ_ERRZERO( err );
01072     // remove fixed vertices
01073     // to avoid really huge arrays (especially since we need to copy
01074     // coords out of annoying MsqVertex class to double array), work
01075     // in blocks
01076     const size_t blocksize = 4096;
01077     std::vector< bool > fixed( blocksize );
01078     MsqVertex coords[blocksize];
01079     double tag_data[3 * blocksize];
01080     for( size_t i = 0; i < vertices.size(); i += blocksize )
01081     {
01082         size_t count = std::min( blocksize, vertices.size() - i );
01083         // remove fixed vertices
01084         mesh->vertices_get_fixed_flag( &vertices[i], fixed, count, err );
01085         MSQ_ERRZERO( err );
01086         size_t w = 0;
01087         for( size_t j = 0; j < count; ++j )
01088             if( !fixed[j] ) vertices[i + w++] = vertices[i + j];
01089         count = w;
01090         // store tag data for free vertices
01091         mesh->vertices_get_coordinates( &vertices[i], coords, count, err );
01092         MSQ_ERRZERO( err );
01093         for( size_t j = 0; j < count; ++j )
01094         {
01095             tag_data[3 * j]     = coords[j][0];
01096             tag_data[3 * j + 1] = coords[j][1];
01097             tag_data[3 * j + 2] = coords[j][2];
01098         }
01099         mesh->tag_set_vertex_data( tag, count, &vertices[i], tag_data, err );
01100         MSQ_ERRZERO( err );
01101     }
01102 
01103     return tag;
01104 }
01105 
01106 void VertexMover::commit_jacobi_coords( TagHandle tag, Mesh* mesh, MsqError& err )
01107 {
01108     Vector3D coords;
01109     std::vector< Mesh::VertexHandle > vertices;
01110     mesh->get_all_vertices( vertices, err );MSQ_ERRRTN( err );
01111     std::vector< bool > fixed( vertices.size() );
01112     mesh->vertices_get_fixed_flag( &vertices[0], fixed, vertices.size(), err );
01113     for( size_t i = 0; i < vertices.size(); ++i )
01114     {
01115         if( !fixed[i] )
01116         {
01117             mesh->tag_get_vertex_data( tag, 1, &vertices[i], coords.to_array(), err );MSQ_ERRRTN( err );
01118             mesh->vertex_set_coordinates( vertices[i], coords, err );MSQ_ERRRTN( err );
01119         }
01120     }
01121 }
01122 
01123 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines