MOAB: Mesh Oriented datABase  (version 5.3.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             {
00712                 mesh->get_all_vertices( patch_vertices, err );
00713             }
00714 
00715             free_vertices.clear();
00716 
00717             for( size_t i = 0; i < patch_vertices.size(); ++i )
00718                 if( !std::binary_search( fixed_vertices.begin(), fixed_vertices.end(), patch_vertices[i] ) )
00719                     free_vertices.push_back( patch_vertices[i] );
00720 
00721             if( free_vertices.empty() )
00722             {  // all vertices were fixed -> skip patch
00723                 continue;
00724             }
00725 
00726             ++num_patches;
00727             all_culled = false;
00728             patch.set_mesh_entities( patch_elements, free_vertices, err );
00729             if( MSQ_CHKERR( err ) )
00730             {
00731                 MSQ_SETERR( perr )( "set_mesh_entities", MsqError::INVALID_STATE );
00732                 PERROR_COND;
00733             }  // goto ERROR;
00734 
00735             // Initialize for inner iteration
00736 
00737             this->initialize_mesh_iteration( patch, err );
00738             if( MSQ_CHKERR( err ) )
00739             {
00740                 MSQ_SETERR( perr )( "initialize_mesh_iteration", MsqError::INVALID_STATE );
00741                 PERROR_COND;
00742             }  // goto ERROR;
00743 
00744             obj_func.reset();
00745 
00746             outer_crit->reset_patch( patch, err );
00747             if( MSQ_CHKERR( err ) )
00748             {
00749                 MSQ_SETERR( perr )( "reset_patch outer", MsqError::INVALID_STATE );
00750                 PERROR_COND;
00751             }  // goto ERROR;
00752 
00753             inner_crit->reset_inner( patch, obj_func, err );
00754             if( MSQ_CHKERR( err ) )
00755             {
00756                 MSQ_SETERR( perr )( "reset_inner", MsqError::INVALID_STATE );
00757                 PERROR_COND;
00758             }  // goto ERROR;
00759 
00760             inner_crit->reset_patch( patch, err );
00761             if( MSQ_CHKERR( err ) )
00762             {
00763                 MSQ_SETERR( perr )( "inner reset_patch", MsqError::INVALID_STATE );
00764                 PERROR_COND;
00765             }  // goto ERROR;
00766 
00767             // Don't even call optimizer if inner termination
00768             // criterion has already been met.
00769             if( !inner_crit->terminate() )
00770             {
00771                 inner_crit_terminated = false;
00772                 if( one_patch ) ++inner_iter;
00773 
00774                 // Call optimizer - should loop on inner_crit->terminate()
00775                 // size_t num_vert=patch.num_free_vertices();
00776                 // std::cout << "P[" << get_parallel_rank() << "] tmp srk VertexMover num_vert= " <<
00777                 // num_vert << std::endl;
00778 
00779                 this->optimize_vertex_positions( patch, err );
00780                 if( MSQ_CHKERR( err ) )
00781                 {
00782                     MSQ_SETERR( perr )( "optimize_vertex_positions", MsqError::INVALID_STATE );
00783                     PERROR_COND;
00784                 }  // goto ERROR;
00785 
00786                 // Update for changes during inner iteration
00787                 // (during optimizer loop)
00788 
00789                 outer_crit->accumulate_patch( patch, err );
00790                 if( MSQ_CHKERR( err ) )
00791                 {
00792                     MSQ_SETERR( perr )( "outer accumulate_patch", MsqError::INVALID_STATE );
00793                     PERROR_COND;
00794                 }  // goto ERROR;
00795 
00796                 inner_crit->cull_vertices( patch, obj_func, err );
00797                 if( MSQ_CHKERR( err ) )
00798                 {
00799                     MSQ_SETERR( perr )( "inner cull_vertices", MsqError::INVALID_STATE );
00800                     PERROR_COND;
00801                 }  // goto ERROR;
00802 
00803                 // experimental...
00804                 if( 0 )
00805                 {
00806                     inner_crit->cull_vertices_global( patch, mesh, domain, settings, obj_func, err );
00807                     if( MSQ_CHKERR( err ) )
00808                     {
00809                         MSQ_SETERR( perr )( "cull_vertices_global", MsqError::INVALID_STATE );
00810                         PERROR_COND;
00811                     }  // goto ERROR;
00812                 }
00813 
00814                 patch.update_mesh( err, coord_tag_ptr );
00815                 if( MSQ_CHKERR( err ) )
00816                 {
00817                     MSQ_SETERR( perr )( "update_mesh", MsqError::INVALID_STATE );
00818                     PERROR_COND;
00819                 }  // goto ERROR;
00820             }
00821         }  // while(p_iter....
00822 
00823         save_or_restore_debug_state( false );
00824 
00825         if( 1 )
00826         {
00827             bool pdone_inner = false;
00828 
00829             if( MSQ_CHKERR( perr ) )
00830             {
00831                 std::cout << "P[" << get_parallel_rank()
00832                           << "] VertexMover::loop_over_mesh found parallel error: " << perr
00833                           << "\n quitting... pdone_inner= " << pdone_inner << std::endl;
00834                 pdone_inner = true;
00835             }
00836 
00837             helper->communicate_any_true( pdone_inner, err );
00838 
00839             if( 0 )
00840                 std::cout << "P[" << get_parallel_rank() << "] tmp srk inner_iter= " << inner_iter
00841                           << " outer_iter= " << outer_iter << " pdone_inner= " << pdone_inner << std::endl;
00842 
00843             if( pdone_inner )
00844             {
00845                 if( 0 )
00846                     std::cout << "P[" << get_parallel_rank()
00847                               << "] tmp srk found parallel error, quitting... pdone_inner= " << pdone_inner
00848                               << std::endl;
00849                 MSQ_SETERR( err )( "PARALLEL ERROR", MsqError::PARALLEL_ERROR );
00850                 break;
00851             }
00852         }
00853 
00854         /// srkenno@sandia.gov save vertex bytes since boundary smoothing changes them
00855         std::vector< unsigned char > saved_bytes;
00856         checkpoint_bytes( mesh, saved_bytes, err );
00857         if( MSQ_CHKERR( err ) )
00858         {
00859             MSQ_SETERR( perr )( "checkpoint_bytes ", MsqError::INVALID_STATE );
00860             PERROR_COND;
00861         }  // goto ERROR;
00862 
00863         helper->communicate_first_independent_set( err );
00864         if( MSQ_CHKERR( err ) )
00865         {
00866             MSQ_SETERR( perr )( "communicate_first_independent_set ", MsqError::INVALID_STATE );
00867             PERROR_COND;
00868         }  // goto ERROR;
00869 
00870         ///*** smooth the boundary ***////
00871         save_or_restore_debug_state( true );
00872         MsqDebug::disable_all();
00873 
00874         while( helper->compute_next_independent_set() )
00875         {
00876             // Loop over all boundary elements
00877             while( helper->get_next_partition_boundary_vertex( vertex_handle ) )
00878             {
00879 
00880                 patch_vertices.clear();
00881                 patch_vertices.push_back( vertex_handle );
00882                 patch_elements.clear();
00883                 mesh->vertices_get_attached_elements( &vertex_handle, 1, patch_elements, junk, err );
00884 
00885                 all_culled = false;
00886                 patch.set_mesh_entities( patch_elements, patch_vertices, err );
00887 
00888                 if( MSQ_CHKERR( err ) )
00889                 {
00890                     MSQ_SETERR( perr )( "set_mesh_entities 2 ", MsqError::INVALID_STATE );
00891                     PERROR_COND;
00892                 }  // goto ERROR;
00893 
00894                 // Initialize for inner iteration
00895 
00896                 this->initialize_mesh_iteration( patch, err );
00897                 if( MSQ_CHKERR( err ) )
00898                 {
00899                     MSQ_SETERR( perr )( " initialize_mesh_iteration 2", MsqError::INVALID_STATE );
00900                     PERROR_COND;
00901                 }  // goto ERROR;
00902 
00903                 obj_func.reset();
00904 
00905                 outer_crit->reset_patch( patch, err );
00906                 if( MSQ_CHKERR( err ) )
00907                 {
00908                     MSQ_SETERR( perr )( "outer reset_patch 2 ", MsqError::INVALID_STATE );
00909                     PERROR_COND;
00910                 }  // goto ERROR;
00911 
00912                 inner_crit->reset_inner( patch, obj_func, err );
00913                 if( MSQ_CHKERR( err ) )
00914                 {
00915                     MSQ_SETERR( perr )( " inner reset_inner 2", MsqError::INVALID_STATE );
00916                     PERROR_COND;
00917                 }  // goto ERROR;
00918 
00919                 inner_crit->reset_patch( patch, err );
00920                 if( MSQ_CHKERR( err ) )
00921                 {
00922                     MSQ_SETERR( perr )( " inner_crit reset_patch 2", MsqError::INVALID_STATE );
00923                     PERROR_COND;
00924                 }  // goto ERROR;
00925 
00926                 // Don't even call optimizer if inner termination
00927                 // criterion has already been met.
00928                 if( !inner_crit->terminate() )
00929                 {
00930                     inner_crit_terminated = false;
00931 
00932                     // Call optimizer - should loop on inner_crit->terminate()
00933                     this->optimize_vertex_positions( patch, err );
00934                     if( MSQ_CHKERR( err ) )
00935                     {
00936                         MSQ_SETERR( perr )
00937                         ( " optimize_vertex_positions 2", MsqError::INVALID_STATE );
00938                         PERROR_COND;
00939                     }  // goto ERROR;
00940 
00941                     // Update for changes during inner iteration
00942                     // (during optimizer loop)
00943 
00944                     outer_crit->accumulate_patch( patch, err );
00945                     if( MSQ_CHKERR( err ) )
00946                     {
00947                         MSQ_SETERR( perr )( " outer accumulate_patch 2", MsqError::INVALID_STATE );
00948                         PERROR_COND;
00949                     }  // goto ERROR;
00950 
00951                     inner_crit->cull_vertices( patch, obj_func, err );
00952                     if( MSQ_CHKERR( err ) )
00953                     {
00954                         MSQ_SETERR( perr )( " inner cull_vertices", MsqError::INVALID_STATE );
00955                         PERROR_COND;
00956                     }  // goto ERROR;
00957 
00958                     // FIXME
00959                     if( 0 )
00960                     {
00961                         inner_crit->cull_vertices_global( patch, mesh, domain, settings, obj_func, err );
00962                         if( MSQ_CHKERR( err ) )
00963                         {
00964                             MSQ_SETERR( perr )
00965                             ( " cull_vertices_global 2 ", MsqError::INVALID_STATE );
00966                             PERROR_COND;
00967                         }  // goto ERROR;
00968                     }
00969 
00970                     patch.update_mesh( err, coord_tag_ptr );
00971                     if( MSQ_CHKERR( err ) )
00972                     {
00973                         MSQ_SETERR( perr )( " update_mesh 2", MsqError::INVALID_STATE );
00974                         PERROR_COND;
00975                     }  // goto ERROR;
00976                 }
00977             }
00978             helper->communicate_next_independent_set( err );
00979             if( MSQ_CHKERR( err ) )
00980             {
00981                 MSQ_SETERR( perr )
00982                 ( " communicate_next_independent_set 2", MsqError::INVALID_STATE );
00983                 PERROR_COND;
00984             }  // goto ERROR;
00985         }      // while(helper->compute_next_independent_set())
00986 
00987         save_or_restore_debug_state( false );
00988 
00989         // if (!get_parallel_rank())
00990         // std::cout << "P[" << get_parallel_rank() << "] tmp srk num_patches= " << num_patches <<
00991         // std::endl;
00992 
00993         /// srkenno@sandia.gov restore vertex bytes since boundary smoothing changes them
00994         restore_bytes( mesh, saved_bytes, err );
00995         if( MSQ_CHKERR( err ) )
00996         {
00997             MSQ_SETERR( perr )( " restore_bytes", MsqError::INVALID_STATE );
00998             PERROR_COND;
00999         }  // goto ERROR;
01000 
01001         if( jacobiOpt ) commit_jacobi_coords( coord_tag, mesh, err );
01002 
01003         this->terminate_mesh_iteration( patch, err );
01004         if( MSQ_CHKERR( err ) )
01005         {
01006             MSQ_SETERR( perr )( " terminate_mesh_iteration", MsqError::INVALID_STATE );
01007             PERROR_COND;
01008         }  // goto ERROR;
01009 
01010         outer_crit->accumulate_outer( mesh, domain, obj_func, settings, err );
01011         if( MSQ_CHKERR( err ) )
01012         {
01013             MSQ_SETERR( perr )( " outer_crit accumulate_outer", MsqError::INVALID_STATE );
01014             PERROR_COND;
01015         }  // goto ERROR;
01016     }
01017 
01018     // ERROR:
01019 
01020     if( MSQ_CHKERR( err ) )
01021     {
01022         std::cout << "P[" << get_parallel_rank() << "] VertexMover::loop_over_mesh error = " << err.error_message()
01023                   << std::endl;
01024     }
01025 
01026     if( jacobiOpt ) mesh->tag_delete( coord_tag, err );
01027 
01028     // call the criteria's cleanup funtions.
01029     outer_crit->cleanup( mesh, domain, err );MSQ_CHKERR( err );
01030     inner_crit->cleanup( mesh, domain, err );MSQ_CHKERR( err );
01031     // call the optimization cleanup function.
01032     this->cleanup();
01033     // close the helper
01034     helper->smoothing_close( err );MSQ_CHKERR( err );
01035 
01036     return 0.;
01037 }
01038 
01039 void VertexMover::initialize_queue( MeshDomainAssoc* mesh_and_domain, const Settings* settings, MsqError& err )
01040 {
01041     QualityImprover::initialize_queue( mesh_and_domain, settings, err );MSQ_ERRRTN( err );
01042     objFuncEval.initialize_queue( mesh_and_domain, settings, err );MSQ_ERRRTN( err );
01043 }
01044 
01045 TagHandle VertexMover::get_jacobi_coord_tag( Mesh* mesh, MsqError& err )
01046 {
01047     // Get tag handle
01048     const char tagname[] = "msq_jacobi_temp_coords";
01049     TagHandle tag        = mesh->tag_create( tagname, Mesh::DOUBLE, 3, 0, err );
01050     MSQ_ERRZERO( err );
01051     /* VertexMover will always delete the tag when it is done, so it is probably
01052        best not to accept an existing tag
01053   if (err.error_code() == TAG_ALREADY_EXISTS) {
01054     err.clear();
01055     tag = tag_get( tagname, err ); MSQ_ERRZERO(err);
01056     std::string name;
01057     Mesh::TagType type;
01058     unsigned length;
01059     mesh->tag_properties( tag, name, type, length, err ); MSQ_ERRZERO(err);
01060     if (type != Mesh::DOUBLE || length != 3) {
01061       MSQ_SETERR(err)(TAG_ALREADY_EXISTS,
01062                      "Tag \"%s\" already exists with invalid type",
01063                      tagname);
01064       return 0;
01065     }
01066   }
01067     */
01068 
01069     // Initialize tag value with initial vertex coordinates so that
01070     // vertices that never get moved end up with the correct coords.
01071     std::vector< Mesh::VertexHandle > vertices;
01072     mesh->get_all_vertices( vertices, err );
01073     MSQ_ERRZERO( err );
01074     // remove fixed vertices
01075     // to avoid really huge arrays (especially since we need to copy
01076     // coords out of annoying MsqVertex class to double array), work
01077     // in blocks
01078     const size_t blocksize = 4096;
01079     std::vector< bool > fixed( blocksize );
01080     MsqVertex coords[blocksize];
01081     double tag_data[3 * blocksize];
01082     for( size_t i = 0; i < vertices.size(); i += blocksize )
01083     {
01084         size_t count = std::min( blocksize, vertices.size() - i );
01085         // remove fixed vertices
01086         mesh->vertices_get_fixed_flag( &vertices[i], fixed, count, err );
01087         MSQ_ERRZERO( err );
01088         size_t w = 0;
01089         for( size_t j = 0; j < count; ++j )
01090             if( !fixed[j] ) vertices[i + w++] = vertices[i + j];
01091         count = w;
01092         // store tag data for free vertices
01093         mesh->vertices_get_coordinates( &vertices[i], coords, count, err );
01094         MSQ_ERRZERO( err );
01095         for( size_t j = 0; j < count; ++j )
01096         {
01097             tag_data[3 * j]     = coords[j][0];
01098             tag_data[3 * j + 1] = coords[j][1];
01099             tag_data[3 * j + 2] = coords[j][2];
01100         }
01101         mesh->tag_set_vertex_data( tag, count, &vertices[i], tag_data, err );
01102         MSQ_ERRZERO( err );
01103     }
01104 
01105     return tag;
01106 }
01107 
01108 void VertexMover::commit_jacobi_coords( TagHandle tag, Mesh* mesh, MsqError& err )
01109 {
01110     Vector3D coords;
01111     std::vector< Mesh::VertexHandle > vertices;
01112     mesh->get_all_vertices( vertices, err );MSQ_ERRRTN( err );
01113     std::vector< bool > fixed( vertices.size() );
01114     mesh->vertices_get_fixed_flag( &vertices[0], fixed, vertices.size(), err );
01115     for( size_t i = 0; i < vertices.size(); ++i )
01116     {
01117         if( !fixed[i] )
01118         {
01119             mesh->tag_get_vertex_data( tag, 1, &vertices[i], coords.to_array(), err );MSQ_ERRRTN( err );
01120             mesh->vertex_set_coordinates( vertices[i], coords, err );MSQ_ERRRTN( err );
01121         }
01122     }
01123 }
01124 
01125 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines