MOAB: Mesh Oriented datABase  (version 5.3.0)
TerminationCriterion.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2004 Sandia Corporation and Argonne National
00005     Laboratory.  Under the terms of Contract DE-AC04-94AL85000
00006     with Sandia Corporation, the U.S. Government retains certain
00007     rights in this software.
00008 
00009     This library is free software; you can redistribute it and/or
00010     modify it under the terms of the GNU Lesser General Public
00011     License as published by the Free Software Foundation; either
00012     version 2.1 of the License, or (at your option) any later version.
00013 
00014     This library is distributed in the hope that it will be useful,
00015     but WITHOUT ANY WARRANTY; without even the implied warranty of
00016     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00017     Lesser General Public License for more details.
00018 
00019     You should have received a copy of the GNU Lesser General Public License
00020     (lgpl.txt) along with this library; if not, write to the Free Software
00021     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00022 
00023     diachin2@llnl.gov, djmelan@sandia.gov, mbrewer@sandia.gov,
00024     pknupp@sandia.gov, tleurent@mcs.anl.gov, tmunson@mcs.anl.gov
00025 
00026   ***************************************************************** */
00027 // -*- Mode : c++; tab-width: 2; c-tab-always-indent: t; indent-tabs-mode: nil; c-basic-offset: 2
00028 // -*-
00029 //
00030 // DESCRIPTION:
00031 // ============
00032 /*! \file TerminationCriterion.cpp
00033 
00034     \brief  Member functions of the MBMesquite::TerminationCriterion class
00035 
00036     \author Michael Brewer
00037     \author Thomas Leurent
00038     \date   Feb. 14, 2003
00039  */
00040 
00041 #include "TerminationCriterion.hpp"
00042 #include "MsqVertex.hpp"
00043 #include "MsqInterrupt.hpp"
00044 #include "OFEvaluator.hpp"
00045 #include "MsqError.hpp"
00046 #include "MsqDebug.hpp"
00047 #include "PatchData.hpp"
00048 #include "MeshWriter.hpp"
00049 #include "MeshUtil.hpp"
00050 #include "SimpleStats.hpp"
00051 
00052 #include <sstream>
00053 #include <set>
00054 
00055 namespace MBMesquite
00056 {
00057 
00058 extern int get_parallel_rank();
00059 extern int get_parallel_size();
00060 extern double reduce_parallel_max( double value );
00061 
00062 #define MSQ_DBGOUT_P0_ONLY( flag ) \
00063     if( !get_parallel_rank() ) MSQ_DBGOUT( flag )
00064 
00065 #define RPM( val ) val
00066 
00067 // this causes race conditions - don't use it
00068 #define RPM1( val ) reduce_parallel_max( val )
00069 
00070 /*! \enum TCType  defines the termination criterion */
00071 enum TCType
00072 {
00073     NONE = 0,
00074     //! checks the gradient \f$\nabla f \f$ of objective function
00075     //! \f$f : I\!\!R^{3N} \rightarrow I\!\!R \f$ against a double \f$d\f$
00076     //! and stops when \f$\sqrt{\sum_{i=1}^{3N}\nabla f_i^2}<d\f$
00077     GRADIENT_L2_NORM_ABSOLUTE = 1 << 0,
00078     //! checks the gradient \f$\nabla f \f$ of objective function
00079     //! \f$f : I\!\!R^{3N} \rightarrow I\!\!R \f$ against a double \f$d\f$
00080     //! and stops when \f$ \max_{i=1}^{3N} \nabla f_i < d \f$
00081     GRADIENT_INF_NORM_ABSOLUTE = 1 << 1,
00082     //! terminates on the j_th iteration when
00083     //! \f$\sqrt{\sum_{i=1}^{3N}\nabla f_{i,j}^2}<d\sqrt{\sum_{i=1}^{3N}\nabla f_{i,0}^2}\f$
00084     //!  That is, terminates when the norm of the gradient is small
00085     //! than some scaling factor times the norm of the original gradient.
00086     GRADIENT_L2_NORM_RELATIVE = 1 << 2,
00087     //! terminates on the j_th iteration when
00088     //! \f$\max_{i=1 \cdots 3N}\nabla f_{i,j}<d \max_{i=1 \cdots 3N}\nabla f_{i,0}\f$
00089     //!  That is, terminates when the norm of the gradient is small
00090     //! than some scaling factor times the norm of the original gradient.
00091     //! (Using the infinity norm.)
00092     GRADIENT_INF_NORM_RELATIVE = 1 << 3,
00093     //! Not yet implemented.
00094     KKT = 1 << 4,
00095     //! Terminates when the objective function value is smaller than
00096     //! the given scalar value.
00097     QUALITY_IMPROVEMENT_ABSOLUTE = 1 << 5,
00098     //! Terminates when the objective function value is smaller than
00099     //! the given scalar value times the original objective function
00100     //! value.
00101     QUALITY_IMPROVEMENT_RELATIVE = 1 << 6,
00102     //! Terminates when the number of iterations exceeds a given integer.
00103     NUMBER_OF_ITERATES = 1 << 7,
00104     //! Terminates when the algorithm exceeds an allotted time limit
00105     //! (given in seconds).
00106     CPU_TIME = 1 << 8,
00107     //! Terminates when a the maximum distance moved by any vertex
00108     //! during the previous iteration is below the given value.
00109     VERTEX_MOVEMENT_ABSOLUTE = 1 << 9,
00110     //! Terminates when a the maximum distance moved by any vertex
00111     //! during the previous iteration is below a value calculated
00112     //! from the edge lengths of the initial mesh.
00113     VERTEX_MOVEMENT_ABS_EDGE_LENGTH = 1 << 10,
00114     //! Terminates when a the maximum distance moved by any vertex
00115     //! during the previous iteration is below the given value
00116     //! times the maximum distance moved by any vertex over the
00117     //! entire course of the optimization.
00118     VERTEX_MOVEMENT_RELATIVE = 1 << 11,
00119     //! Terminates when the decrease in the objective function value since
00120     //! the previous iteration is below the given value.
00121     SUCCESSIVE_IMPROVEMENTS_ABSOLUTE = 1 << 12,
00122     //! Terminates when the decrease in the objective function value since
00123     //! the previous iteration is below the given value times the
00124     //! decrease in the objective function value since the beginning
00125     //! of this optimization process.
00126     SUCCESSIVE_IMPROVEMENTS_RELATIVE = 1 << 13,
00127     //! Terminates when any vertex leaves the bounding box, defined
00128     //! by the given value, d.  That is, when the absolute value of
00129     //! a single coordinate of vertex's position exceeds d.
00130     BOUNDED_VERTEX_MOVEMENT = 1 << 14,
00131     //! Terminate when no elements are inverted
00132     UNTANGLED_MESH = 1 << 15
00133 };
00134 
00135 const unsigned long GRAD_FLAGS =
00136     GRADIENT_L2_NORM_ABSOLUTE | GRADIENT_INF_NORM_ABSOLUTE | GRADIENT_L2_NORM_RELATIVE | GRADIENT_INF_NORM_RELATIVE;
00137 const unsigned long OF_FLAGS = QUALITY_IMPROVEMENT_ABSOLUTE | QUALITY_IMPROVEMENT_RELATIVE |
00138                                SUCCESSIVE_IMPROVEMENTS_ABSOLUTE | SUCCESSIVE_IMPROVEMENTS_RELATIVE;
00139 
00140 const unsigned long MOVEMENT_FLAGS =
00141     VERTEX_MOVEMENT_ABSOLUTE | VERTEX_MOVEMENT_ABS_EDGE_LENGTH | VERTEX_MOVEMENT_RELATIVE;
00142 
00143 /*!Constructor initializes all of the data members which are not
00144   necessarily automatically initialized in their constructors.*/
00145 TerminationCriterion::TerminationCriterion( std::string name, InnerOuterType pinnerOuterType )
00146     : mGrad( 8 ), initialVerticesMemento( 0 ), previousVerticesMemento( 0 ), debugLevel( 2 ),
00147       timeStepFileType( NOTYPE ), moniker( name ), innerOuterType( pinnerOuterType )
00148 {
00149     terminationCriterionFlag = NONE;
00150     cullingMethodFlag        = NONE;
00151     cullingEps               = 0.0;
00152     cullingGlobalPatch       = false;
00153     initialOFValue           = 0.0;
00154     previousOFValue          = 0.0;
00155     currentOFValue           = 0.0;
00156     lowerOFBound             = 0.0;
00157     initialGradL2NormSquared = 0.0;
00158     initialGradInfNorm       = 0.0;
00159     // initial size of the gradient array
00160     gradL2NormAbsoluteEpsSquared      = 0.0;
00161     gradL2NormRelativeEpsSquared      = 0.0;
00162     gradInfNormAbsoluteEps            = 0.0;
00163     gradInfNormRelativeEps            = 0.0;
00164     qualityImprovementAbsoluteEps     = 0.0;
00165     qualityImprovementRelativeEps     = 0.0;
00166     iterationBound                    = 0;
00167     iterationCounter                  = 0;
00168     timeBound                         = 0.0;
00169     vertexMovementAbsoluteEps         = 0.0;
00170     vertexMovementRelativeEps         = 0.0;
00171     vertexMovementAvgBeta             = -1;
00172     vertexMovementAbsoluteAvgEdge     = -1;
00173     successiveImprovementsAbsoluteEps = 0.0;
00174     successiveImprovementsRelativeEps = 0.0;
00175     boundedVertexMovementEps          = 0.0;
00176     globalInvertedCount               = 0;
00177     patchInvertedCount                = 0;
00178 }
00179 
00180 std::string TerminationCriterion::par_string()
00181 {
00182     if( get_parallel_size() )
00183     {
00184         std::ostringstream str;
00185         str << "P[" << get_parallel_rank() << "] " + moniker + " ";
00186         return str.str();
00187     }
00188     return moniker + " ";
00189 }
00190 
00191 void TerminationCriterion::add_absolute_gradient_L2_norm( double eps )
00192 {
00193     terminationCriterionFlag |= GRADIENT_L2_NORM_ABSOLUTE;
00194     gradL2NormAbsoluteEpsSquared = eps * eps;
00195 }
00196 
00197 void TerminationCriterion::add_absolute_gradient_inf_norm( double eps )
00198 {
00199     terminationCriterionFlag |= GRADIENT_INF_NORM_ABSOLUTE;
00200     gradInfNormAbsoluteEps = eps;
00201 }
00202 
00203 void TerminationCriterion::add_relative_gradient_L2_norm( double eps )
00204 {
00205     terminationCriterionFlag |= GRADIENT_L2_NORM_RELATIVE;
00206     gradL2NormRelativeEpsSquared = eps * eps;
00207 }
00208 
00209 void TerminationCriterion::add_relative_gradient_inf_norm( double eps )
00210 {
00211     terminationCriterionFlag |= GRADIENT_INF_NORM_RELATIVE;
00212     gradInfNormRelativeEps = eps;
00213 }
00214 
00215 void TerminationCriterion::add_absolute_quality_improvement( double eps )
00216 {
00217     terminationCriterionFlag |= QUALITY_IMPROVEMENT_ABSOLUTE;
00218     qualityImprovementAbsoluteEps = eps;
00219 }
00220 
00221 void TerminationCriterion::add_relative_quality_improvement( double eps )
00222 {
00223     terminationCriterionFlag |= QUALITY_IMPROVEMENT_RELATIVE;
00224     qualityImprovementRelativeEps = eps;
00225 }
00226 
00227 void TerminationCriterion::add_absolute_vertex_movement( double eps )
00228 {
00229     terminationCriterionFlag |= VERTEX_MOVEMENT_ABSOLUTE;
00230     // we actually compare squared movement to squared epsilon
00231     vertexMovementAbsoluteEps = ( eps * eps );
00232 }
00233 
00234 void TerminationCriterion::add_absolute_vertex_movement_edge_length( double beta )
00235 {
00236     terminationCriterionFlag |= VERTEX_MOVEMENT_ABS_EDGE_LENGTH;
00237     vertexMovementAvgBeta = beta;
00238 }
00239 
00240 void TerminationCriterion::add_relative_vertex_movement( double eps )
00241 {
00242     terminationCriterionFlag |= VERTEX_MOVEMENT_RELATIVE;
00243     // we actually compare squared movement to squared epsilon
00244     vertexMovementRelativeEps = ( eps * eps );
00245 }
00246 
00247 void TerminationCriterion::add_absolute_successive_improvement( double eps )
00248 {
00249     terminationCriterionFlag |= SUCCESSIVE_IMPROVEMENTS_ABSOLUTE;
00250     successiveImprovementsAbsoluteEps = eps;
00251 }
00252 
00253 void TerminationCriterion::add_relative_successive_improvement( double eps )
00254 {
00255     terminationCriterionFlag |= SUCCESSIVE_IMPROVEMENTS_RELATIVE;
00256     successiveImprovementsRelativeEps = eps;
00257 }
00258 
00259 void TerminationCriterion::add_cpu_time( double seconds )
00260 {
00261     terminationCriterionFlag |= CPU_TIME;
00262     timeBound = seconds;
00263 }
00264 
00265 void TerminationCriterion::add_iteration_limit( unsigned int max_iterations )
00266 {
00267     terminationCriterionFlag |= NUMBER_OF_ITERATES;
00268     iterationBound = max_iterations;
00269 }
00270 
00271 void TerminationCriterion::add_bounded_vertex_movement( double eps )
00272 {
00273     terminationCriterionFlag |= BOUNDED_VERTEX_MOVEMENT;
00274     boundedVertexMovementEps = eps;
00275 }
00276 
00277 void TerminationCriterion::add_untangled_mesh()
00278 {
00279     terminationCriterionFlag |= UNTANGLED_MESH;
00280 }
00281 
00282 void TerminationCriterion::remove_all_criteria()
00283 {
00284     terminationCriterionFlag = 0;
00285 }
00286 
00287 void TerminationCriterion::cull_on_absolute_quality_improvement( double limit )
00288 {
00289     cullingMethodFlag = QUALITY_IMPROVEMENT_ABSOLUTE;
00290     cullingEps        = limit;
00291 }
00292 void TerminationCriterion::cull_on_relative_quality_improvement( double limit )
00293 {
00294     cullingMethodFlag = QUALITY_IMPROVEMENT_RELATIVE;
00295     cullingEps        = limit;
00296 }
00297 void TerminationCriterion::cull_on_absolute_vertex_movement( double limit )
00298 {
00299     cullingMethodFlag = VERTEX_MOVEMENT_ABSOLUTE;
00300     cullingEps        = limit;
00301 }
00302 void TerminationCriterion::cull_on_relative_vertex_movement( double limit )
00303 {
00304     cullingMethodFlag = VERTEX_MOVEMENT_RELATIVE;
00305     cullingEps        = limit;
00306 }
00307 void TerminationCriterion::cull_on_absolute_vertex_movement_edge_length( double limit )
00308 {
00309     cullingMethodFlag     = VERTEX_MOVEMENT_ABS_EDGE_LENGTH;
00310     vertexMovementAvgBeta = limit;
00311 }
00312 void TerminationCriterion::cull_on_absolute_successive_improvement( double limit )
00313 {
00314     cullingMethodFlag = SUCCESSIVE_IMPROVEMENTS_ABSOLUTE;
00315     cullingEps        = limit;
00316 }
00317 void TerminationCriterion::cull_on_relative_successive_improvement( double limit )
00318 {
00319     cullingMethodFlag = SUCCESSIVE_IMPROVEMENTS_RELATIVE;
00320     cullingEps        = limit;
00321 }
00322 void TerminationCriterion::cull_untangled_mesh()
00323 {
00324     cullingMethodFlag = UNTANGLED_MESH;
00325 }
00326 
00327 void TerminationCriterion::remove_culling()
00328 {
00329     cullingMethodFlag = NONE;
00330 }
00331 
00332 void TerminationCriterion::cull_for_global_patch( bool val )
00333 {
00334     cullingGlobalPatch = val;
00335 }
00336 
00337 /*!This version of reset is called using a MeshSet, which implies
00338   it is only called when this criterion is used as the 'outer' termination
00339   criterion.
00340  */
00341 void TerminationCriterion::reset_outer( Mesh* mesh, MeshDomain* domain, OFEvaluator& obj_eval, const Settings* settings,
00342                                         MsqError& err )
00343 {
00344     const unsigned long totalFlag = terminationCriterionFlag | cullingMethodFlag;
00345     PatchData global_patch;
00346     if( settings ) global_patch.attach_settings( settings );
00347 
00348     // if we need to fill out the global patch data object.
00349     if( ( totalFlag & ( GRAD_FLAGS | OF_FLAGS | VERTEX_MOVEMENT_RELATIVE | UNTANGLED_MESH ) ) || timeStepFileType )
00350     {
00351         global_patch.set_mesh( mesh );
00352         global_patch.set_domain( domain );
00353         global_patch.fill_global_patch( err );MSQ_ERRRTN( err );
00354     }
00355 
00356     // now call the other reset
00357     reset_inner( global_patch, obj_eval, err );MSQ_ERRRTN( err );
00358 }
00359 
00360 /*!Reset function using using a PatchData object.  This function is
00361   called for the inner-stopping criterion directly from the
00362   loop over mesh function in VertexMover.  For outer criterion,
00363   it is called from the reset function which takes a MeshSet object.
00364   This function prepares the object to be used by setting the initial
00365   values of some of the data members.  As examples, if needed, it resets
00366   the cpu timer to zero, the iteration counter to zero, and the
00367   initial and previous objective function values to the current
00368   objective function value for this patch.
00369   The return value for this function is similar to that of terminate().
00370   The function returns false if the checked criteria have not been
00371   satisfied, and true if they have been.  reset() only checks the
00372   GRADIENT_INF_NORM_ABSOLUTE, GRADIENT_L2_NORM_ABSOLUTE, and the
00373   QUALITY_IMPROVEMENT_ABSOLUTE criteria.  Checking these criteria
00374   allows the QualityImprover to skip the entire optimization if
00375   the initial mesh satisfies the appropriate conditions.
00376  */
00377 void TerminationCriterion::reset_inner( PatchData& pd, OFEvaluator& obj_eval, MsqError& err )
00378 {
00379     const unsigned long totalFlag = terminationCriterionFlag | cullingMethodFlag;
00380 
00381     // clear flag for BOUNDED_VERTEX_MOVEMENT
00382     vertexMovementExceedsBound = 0;
00383 
00384     // Use -1 to denote that this isn't initialized yet.
00385     // As all valid values must be >= 0.0, a negative
00386     // value indicates that it is uninitialized and is
00387     // always less than any valid value.
00388     maxSquaredMovement = -1;
00389 
00390     // Clear the iteration count.
00391     iterationCounter = 0;
00392 
00393     // reset the inner timer if needed
00394     if( totalFlag & CPU_TIME ) { mTimer.reset(); }
00395 
00396     // GRADIENT
00397     currentGradInfNorm = initialGradInfNorm = 0.0;
00398     currentGradL2NormSquared = initialGradL2NormSquared = 0.0;
00399     if( totalFlag & GRAD_FLAGS )
00400     {
00401         if( !obj_eval.have_objective_function() )
00402         {
00403             MSQ_SETERR( err )
00404             ( "Error termination criteria set which uses objective "
00405               "functions, but no objective function is available.",
00406               MsqError::INVALID_STATE );
00407             return;
00408         }
00409         int num_vertices = pd.num_free_vertices();
00410         mGrad.resize( num_vertices );
00411 
00412         // get gradient and make sure it is valid
00413         bool b = obj_eval.evaluate( pd, currentOFValue, mGrad, err );MSQ_ERRRTN( err );
00414         if( !b )
00415         {
00416             MSQ_SETERR( err )
00417             ( "Initial patch is invalid for gradient computation.", MsqError::INVALID_STATE );
00418             return;
00419         }
00420 
00421         // get the gradient norms
00422         if( totalFlag & ( GRADIENT_INF_NORM_ABSOLUTE | GRADIENT_INF_NORM_RELATIVE ) )
00423         {
00424             currentGradInfNorm = initialGradInfNorm = Linf( mGrad );
00425             MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o Initial gradient Inf norm: "
00426                                              << " " << RPM( initialGradInfNorm ) << std::endl;
00427         }
00428 
00429         if( totalFlag & ( GRADIENT_L2_NORM_ABSOLUTE | GRADIENT_L2_NORM_RELATIVE ) )
00430         {
00431             currentGradL2NormSquared = initialGradL2NormSquared = length_squared( mGrad );
00432             MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o Initial gradient L2 norm: "
00433                                              << " " << RPM( std::sqrt( initialGradL2NormSquared ) ) << std::endl;
00434         }
00435 
00436         // the OFvalue comes for free, so save it
00437         previousOFValue = currentOFValue;
00438         initialOFValue  = currentOFValue;
00439     }
00440     // find the initial objective function value if needed and not already
00441     // computed.  If we needed the gradient, we have the OF value for free.
00442     // Also, if possible, get initial OF value if writing plot file.  Solvers
00443     // often supply the OF value for subsequent iterations so by calculating
00444     // the initial value we can generate OF value plots.
00445     else if( ( totalFlag & OF_FLAGS ) ||
00446              ( plotFile.is_open() && pd.num_free_vertices() && obj_eval.have_objective_function() ) )
00447     {
00448         // ensure the obj_ptr is not null
00449         if( !obj_eval.have_objective_function() )
00450         {
00451             MSQ_SETERR( err )
00452             ( "Error termination criteria set which uses objective "
00453               "functions, but no objective function is available.",
00454               MsqError::INVALID_STATE );
00455             return;
00456         }
00457 
00458         bool b = obj_eval.evaluate( pd, currentOFValue, err );MSQ_ERRRTN( err );
00459         if( !b )
00460         {
00461             MSQ_SETERR( err )
00462             ( "Initial patch is invalid for evaluation.", MsqError::INVALID_STATE );
00463             return;
00464         }
00465         // std::cout<<"\nReseting initial of value = "<<initialOFValue;
00466         previousOFValue = currentOFValue;
00467         initialOFValue  = currentOFValue;
00468     }
00469 
00470     if( totalFlag & ( GRAD_FLAGS | OF_FLAGS ) )
00471         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o Initial OF value: "
00472                                          << " " << RPM( initialOFValue ) << std::endl;
00473 
00474     // Store current vertex locations now, because we'll
00475     // need them later to compare the current movement with.
00476     if( totalFlag & VERTEX_MOVEMENT_RELATIVE )
00477     {
00478         if( initialVerticesMemento ) { pd.recreate_vertices_memento( initialVerticesMemento, err ); }
00479         else
00480         {
00481             initialVerticesMemento = pd.create_vertices_memento( err );
00482         }
00483         MSQ_ERRRTN( err );
00484         maxSquaredInitialMovement = DBL_MAX;
00485     }
00486     else
00487     {
00488         maxSquaredInitialMovement = 0;
00489     }
00490 
00491     if( terminationCriterionFlag & UNTANGLED_MESH )
00492     {
00493         globalInvertedCount = count_inverted( pd, err );
00494         // if (innerOuterType==TYPE_OUTER) MSQ_DBGOUT_P0_ONLY(debugLevel) << par_string() << "  o
00495         // Num Inverted: " << " " << globalInvertedCount << std::endl;
00496         patchInvertedCount = 0;MSQ_ERRRTN( err );
00497     }
00498 
00499     if( timeStepFileType )
00500     {
00501         // If didn't already calculate gradient abive, calculate it now.
00502         if( !( totalFlag & GRAD_FLAGS ) )
00503         {
00504             mGrad.resize( pd.num_free_vertices() );
00505             obj_eval.evaluate( pd, currentOFValue, mGrad, err );
00506             err.clear();
00507         }
00508         write_timestep( pd, mGrad.empty() ? 0 : arrptr( mGrad ), err );
00509     }
00510 
00511     if( plotFile.is_open() )
00512     {
00513         // two newlines so GNU plot knows that we are starting a new data set
00514         plotFile << std::endl << std::endl;
00515         // write column headings as comment in data file
00516         plotFile << "#Iter\tCPU\tObjFunc\tGradL2\tGradInf\tMovement\tInverted" << std::endl;
00517         // write initial values
00518         plotFile << 0 << '\t' << mTimer.since_birth() << '\t' << initialOFValue << '\t'
00519                  << std::sqrt( currentGradL2NormSquared ) << '\t' << currentGradInfNorm << '\t' << 0.0 << '\t'
00520                  << globalInvertedCount << std::endl;
00521     }
00522 }
00523 
00524 void TerminationCriterion::reset_patch( PatchData& pd, MsqError& err )
00525 {
00526     const unsigned long totalFlag = terminationCriterionFlag | cullingMethodFlag;
00527     if( totalFlag & MOVEMENT_FLAGS )
00528     {
00529         if( previousVerticesMemento )
00530             pd.recreate_vertices_memento( previousVerticesMemento, err );
00531         else
00532             previousVerticesMemento = pd.create_vertices_memento( err );MSQ_ERRRTN( err );
00533     }
00534 
00535     if( totalFlag & UNTANGLED_MESH )
00536     {
00537         patchInvertedCount = count_inverted( pd, err );
00538         // MSQ_DBGOUT_P0_ONLY(debugLevel) << par_string() << "  o Num Patch Inverted: " << " " <<
00539         // patchInvertedCount << std::endl;MSQ_ERRRTN( err );
00540     }
00541 }
00542 
00543 void TerminationCriterion::accumulate_inner( PatchData& pd, OFEvaluator& of_eval, MsqError& err )
00544 {
00545     double of_value = 0;
00546 
00547     if( terminationCriterionFlag & GRAD_FLAGS )
00548     {
00549         mGrad.resize( pd.num_free_vertices() );
00550         bool b = of_eval.evaluate( pd, of_value, mGrad, err );MSQ_ERRRTN( err );
00551         if( !b )
00552         {
00553             MSQ_SETERR( err )
00554             ( "Initial patch is invalid for gradient compuation.", MsqError::INVALID_MESH );
00555             return;
00556         }
00557     }
00558     else if( terminationCriterionFlag & OF_FLAGS )
00559     {
00560         bool b = of_eval.evaluate( pd, of_value, err );MSQ_ERRRTN( err );
00561         if( !b )
00562         {
00563             MSQ_SETERR( err )
00564             ( "Invalid patch passed to TerminationCriterion.", MsqError::INVALID_MESH );
00565             return;
00566         }
00567     }
00568 
00569     accumulate_inner( pd, of_value, mGrad.empty() ? 0 : arrptr( mGrad ), err );MSQ_CHKERR( err );
00570 }
00571 
00572 void TerminationCriterion::accumulate_inner( PatchData& pd, double of_value, Vector3D* grad_array, MsqError& err )
00573 {
00574     // if terminating on the norm of the gradient
00575     // currentGradL2NormSquared = HUGE_VAL;
00576     if( terminationCriterionFlag & ( GRADIENT_L2_NORM_ABSOLUTE | GRADIENT_L2_NORM_RELATIVE ) )
00577     {
00578         currentGradL2NormSquared = length_squared( grad_array, pd.num_free_vertices() );  // get the L2 norm
00579         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o Info -- gradient L2 norm: "
00580                                          << " " << RPM( std::sqrt( currentGradL2NormSquared ) ) << std::endl;
00581     }
00582     // currentGradInfNorm = 10e6;
00583     if( terminationCriterionFlag & ( GRADIENT_INF_NORM_ABSOLUTE | GRADIENT_INF_NORM_RELATIVE ) )
00584     {
00585         currentGradInfNorm = Linf( grad_array, pd.num_free_vertices() );  // get the Linf norm
00586         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o Info -- gradient Inf norm: "
00587                                          << " " << RPM( currentGradInfNorm ) << std::endl;
00588     }
00589 
00590     if( terminationCriterionFlag & VERTEX_MOVEMENT_RELATIVE )
00591     {
00592         maxSquaredInitialMovement = pd.get_max_vertex_movement_squared( initialVerticesMemento, err );MSQ_ERRRTN( err );
00593         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o Info -- max initial vertex movement: "
00594                                          << " " << RPM( maxSquaredInitialMovement ) << std::endl;
00595     }
00596 
00597     previousOFValue = currentOFValue;
00598     currentOFValue  = of_value;
00599     if( terminationCriterionFlag & OF_FLAGS )
00600     {
00601         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o Info -- OF Value: "
00602                                          << " " << RPM( of_value ) << " iterationCounter= " << iterationCounter
00603                                          << std::endl;
00604     }
00605     else if( grad_array )
00606     {
00607         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o OF Value: "
00608                                          << " " << RPM( of_value ) << " iterationCounter= "
00609                                          << iterationCounter
00610                                          //<< " terminationCriterionFlag= " <<
00611                                          // terminationCriterionFlag << " OF_FLAGS = " << OF_FLAGS
00612                                          << std::endl;
00613     }
00614 
00615     ++iterationCounter;
00616     if( timeStepFileType ) write_timestep( pd, grad_array, err );
00617 
00618     if( plotFile.is_open() )
00619         plotFile << iterationCounter << '\t' << mTimer.since_birth() << '\t' << of_value << '\t'
00620                  << std::sqrt( currentGradL2NormSquared ) << '\t' << currentGradInfNorm << '\t'
00621                  << ( maxSquaredMovement > 0.0 ? std::sqrt( maxSquaredMovement ) : 0.0 ) << '\t' << globalInvertedCount
00622                  << std::endl;
00623 }
00624 
00625 void TerminationCriterion::accumulate_outer( Mesh* mesh, MeshDomain* domain, OFEvaluator& of_eval,
00626                                              const Settings* settings, MsqError& err )
00627 {
00628     PatchData global_patch;
00629     if( settings ) global_patch.attach_settings( settings );
00630 
00631     // if we need to fill out the global patch data object.
00632     if( ( terminationCriterionFlag & ( GRAD_FLAGS | OF_FLAGS | VERTEX_MOVEMENT_RELATIVE ) ) || timeStepFileType )
00633     {
00634         global_patch.set_mesh( mesh );
00635         global_patch.set_domain( domain );
00636         global_patch.fill_global_patch( err );MSQ_ERRRTN( err );
00637     }
00638 
00639     accumulate_inner( global_patch, of_eval, err );MSQ_ERRRTN( err );
00640 }
00641 
00642 void TerminationCriterion::accumulate_patch( PatchData& pd, MsqError& err )
00643 {
00644     if( terminationCriterionFlag & MOVEMENT_FLAGS )
00645     {
00646         double patch_max_dist = pd.get_max_vertex_movement_squared( previousVerticesMemento, err );
00647         if( patch_max_dist > maxSquaredMovement ) maxSquaredMovement = patch_max_dist;
00648         pd.recreate_vertices_memento( previousVerticesMemento, err );MSQ_ERRRTN( err );
00649     }
00650 
00651     // if terminating on bounded vertex movement (a bounding box for the mesh)
00652     if( terminationCriterionFlag & BOUNDED_VERTEX_MOVEMENT )
00653     {
00654         const MsqVertex* vert = pd.get_vertex_array( err );
00655         int num_vert          = pd.num_free_vertices();
00656         int i                 = 0;
00657         // for each vertex
00658         for( i = 0; i < num_vert; ++i )
00659         {
00660             // if any of the coordinates are greater than eps
00661             if( ( vert[i][0] > boundedVertexMovementEps ) || ( vert[i][1] > boundedVertexMovementEps ) ||
00662                 ( vert[i][2] > boundedVertexMovementEps ) )
00663             {
00664                 ++vertexMovementExceedsBound;
00665             }
00666         }
00667     }
00668 
00669     if( ( terminationCriterionFlag | cullingMethodFlag ) & UNTANGLED_MESH )
00670     {
00671         size_t new_count = count_inverted( pd, err );
00672         // be careful here because size_t is unsigned
00673         globalInvertedCount += new_count;
00674         globalInvertedCount -= patchInvertedCount;
00675         patchInvertedCount = new_count;
00676         // if (innerOuterType==TYPE_OUTER)
00677         //  MSQ_DBGOUT_P0_ONLY(debugLevel) << par_string() << "  o Num Patch Inverted: " << " " <<
00678         //  patchInvertedCount << " globalInvertedCount= " << globalInvertedCount << std::endl;
00679 
00680         MSQ_ERRRTN( err );
00681     }
00682 }
00683 
00684 /*!  This function evaluates the needed information and then evaluates
00685   the termination criteria.  If any of the selected criteria are satisfied,
00686   the function returns true.  Otherwise, the function returns false.
00687  */
00688 bool TerminationCriterion::terminate()
00689 {
00690     bool return_flag = false;
00691     // std::cout<<"\nInside terminate(pd,of,err):  flag = "<<terminationCriterionFlag << std::endl;
00692     int type = 0;
00693 
00694     // First check for an interrupt signal
00695     if( MsqInterrupt::interrupt() )
00696     {
00697         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- INTERRUPTED"
00698                                          << " " << std::endl;
00699         return true;
00700     }
00701 
00702     // if terminating on numbering of inner iterations
00703     if( NUMBER_OF_ITERATES & terminationCriterionFlag && iterationCounter >= iterationBound )
00704     {
00705         return_flag = true;
00706         type        = 1;
00707         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- Reached "
00708                                          << " " << iterationBound << " iterations." << std::endl;
00709     }
00710 
00711     if( CPU_TIME & terminationCriterionFlag && mTimer.since_birth() >= timeBound )
00712     {
00713         return_flag = true;
00714         type        = 2;
00715         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- Exceeded CPU time. "
00716                                          << " " << std::endl;
00717     }
00718 
00719     if( MOVEMENT_FLAGS & terminationCriterionFlag && maxSquaredMovement >= 0.0 )
00720     {
00721         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o Info -- Maximum vertex movement: "
00722                                          << " " << RPM( sqrt( maxSquaredMovement ) ) << std::endl;
00723 
00724         if( VERTEX_MOVEMENT_ABSOLUTE & terminationCriterionFlag && maxSquaredMovement <= vertexMovementAbsoluteEps )
00725         {
00726             MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- VERTEX_MOVEMENT_ABSOLUTE: "
00727                                              << " " << RPM( sqrt( maxSquaredMovement ) ) << std::endl;
00728             return_flag = true;
00729             type        = 3;
00730         }
00731 
00732         if( VERTEX_MOVEMENT_RELATIVE & terminationCriterionFlag &&
00733             maxSquaredMovement <= vertexMovementRelativeEps * maxSquaredInitialMovement )
00734         {
00735             MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- VERTEX_MOVEMENT_RELATIVE: "
00736                                              << " " << RPM( sqrt( maxSquaredMovement ) ) << std::endl;
00737             return_flag = true;
00738             type        = 4;
00739         }
00740 
00741         if( VERTEX_MOVEMENT_ABS_EDGE_LENGTH & terminationCriterionFlag )
00742         {
00743             assert( vertexMovementAbsoluteAvgEdge > -1e-12 );  // make sure value actually got calculated
00744             if( maxSquaredMovement <= vertexMovementAbsoluteAvgEdge )
00745             {
00746                 MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- VERTEX_MOVEMENT_ABS_EDGE_LENGTH: "
00747                                                  << " " << RPM( sqrt( maxSquaredMovement ) ) << std::endl;
00748                 return_flag = true;
00749                 type        = 5;
00750             }
00751         }
00752     }
00753 
00754     if( GRADIENT_L2_NORM_ABSOLUTE & terminationCriterionFlag &&
00755         currentGradL2NormSquared <= gradL2NormAbsoluteEpsSquared )
00756     {
00757         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- GRADIENT_L2_NORM_ABSOLUTE: "
00758                                          << " " << RPM( currentGradL2NormSquared ) << std::endl;
00759         return_flag = true;
00760         type        = 6;
00761     }
00762 
00763     if( GRADIENT_INF_NORM_ABSOLUTE & terminationCriterionFlag && currentGradInfNorm <= gradInfNormAbsoluteEps )
00764     {
00765         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- GRADIENT_INF_NORM_ABSOLUTE: "
00766                                          << " " << RPM( currentGradInfNorm ) << std::endl;
00767         return_flag = true;
00768         type        = 7;
00769     }
00770 
00771     if( GRADIENT_L2_NORM_RELATIVE & terminationCriterionFlag &&
00772         currentGradL2NormSquared <= ( gradL2NormRelativeEpsSquared * initialGradL2NormSquared ) )
00773     {
00774         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- GRADIENT_L2_NORM_RELATIVE: "
00775                                          << " " << RPM( currentGradL2NormSquared ) << std::endl;
00776         return_flag = true;
00777         type        = 8;
00778     }
00779 
00780     if( GRADIENT_INF_NORM_RELATIVE & terminationCriterionFlag &&
00781         currentGradInfNorm <= ( gradInfNormRelativeEps * initialGradInfNorm ) )
00782     {
00783         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- GRADIENT_INF_NORM_RELATIVE: "
00784                                          << " " << RPM( currentGradInfNorm ) << std::endl;
00785         return_flag = true;
00786         type        = 9;
00787     }
00788     // Quality Improvement and Successive Improvements are below.
00789     // The relative forms are only valid after the first iteration.
00790     if( ( QUALITY_IMPROVEMENT_ABSOLUTE & terminationCriterionFlag ) && currentOFValue <= qualityImprovementAbsoluteEps )
00791     {
00792         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- QUALITY_IMPROVEMENT_ABSOLUTE: "
00793                                          << " " << RPM( currentOFValue ) << std::endl;
00794         return_flag = true;
00795         type        = 10;
00796     }
00797 
00798     // only valid if an iteration has occurred, see above.
00799     if( iterationCounter > 0 )
00800     {
00801         if( SUCCESSIVE_IMPROVEMENTS_ABSOLUTE & terminationCriterionFlag &&
00802             ( previousOFValue - currentOFValue ) <= successiveImprovementsAbsoluteEps )
00803         {
00804             MSQ_DBGOUT_P0_ONLY( debugLevel )
00805                 << par_string() << "  o TermCrit -- SUCCESSIVE_IMPROVEMENTS_ABSOLUTE: previousOFValue= "
00806                 << " " << previousOFValue << " currentOFValue= " << RPM( currentOFValue )
00807                 << " successiveImprovementsAbsoluteEps= " << successiveImprovementsAbsoluteEps << std::endl;
00808             return_flag = true;
00809             type        = 11;
00810         }
00811         if( QUALITY_IMPROVEMENT_RELATIVE & terminationCriterionFlag &&
00812             ( currentOFValue - lowerOFBound ) <= qualityImprovementRelativeEps * ( initialOFValue - lowerOFBound ) )
00813         {
00814             MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- QUALITY_IMPROVEMENT_RELATIVE: "
00815                                              << " " << std::endl;
00816             return_flag = true;
00817             type        = 12;
00818         }
00819         if( SUCCESSIVE_IMPROVEMENTS_RELATIVE & terminationCriterionFlag &&
00820             ( previousOFValue - currentOFValue ) <=
00821                 successiveImprovementsRelativeEps * ( initialOFValue - currentOFValue ) )
00822         {
00823             MSQ_DBGOUT_P0_ONLY( debugLevel )
00824                 << par_string() << "  o TermCrit -- SUCCESSIVE_IMPROVEMENTS_RELATIVE: previousOFValue= "
00825                 << " " << previousOFValue << " currentOFValue= " << RPM( currentOFValue )
00826                 << " successiveImprovementsRelativeEps= " << successiveImprovementsRelativeEps << std::endl;
00827             return_flag = true;
00828             type        = 13;
00829         }
00830     }
00831 
00832     if( BOUNDED_VERTEX_MOVEMENT & terminationCriterionFlag && vertexMovementExceedsBound )
00833     {
00834         return_flag = true;
00835         type        = 14;
00836         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- "
00837                                          << " " << vertexMovementExceedsBound << " vertices out of bounds."
00838                                          << std::endl;
00839     }
00840 
00841     if( UNTANGLED_MESH & terminationCriterionFlag )
00842     {
00843         if( innerOuterType == TYPE_OUTER )
00844         {
00845             // MSQ_DBGOUT_P0_ONLY(debugLevel)
00846             MSQ_DBGOUT( debugLevel ) << par_string() << "  o Num Inverted: "
00847                                      << " " << globalInvertedCount << std::endl;
00848         }
00849         if( !globalInvertedCount )
00850         {
00851             MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- UNTANGLED_MESH: "
00852                                              << " " << std::endl;
00853             return_flag = true;
00854             type        = 15;
00855         }
00856     }
00857 
00858     // clear this value at the end of each iteration
00859     vertexMovementExceedsBound = 0;
00860     maxSquaredMovement         = -1.0;
00861 
00862     if( timeStepFileType == GNUPLOT && return_flag )
00863     {
00864         MsqError err;
00865         MeshWriter::write_gnuplot_overlay( iterationCounter, timeStepFileName.c_str(), err );
00866     }
00867 
00868     if( 0 && return_flag && MSQ_DBG( 2 ) )
00869         std::cout << "P[" << get_parallel_rank() << "] tmp TerminationCriterion::terminate: " << moniker
00870                   << " return_flag= " << return_flag << " type= " << type
00871                   << " terminationCriterionFlag= " << terminationCriterionFlag << " debugLevel= " << debugLevel
00872                   << std::endl;
00873 
00874     // if none of the criteria were satisfied
00875     return return_flag;
00876 }
00877 
00878 bool TerminationCriterion::criterion_is_set()
00879 {
00880     if( !terminationCriterionFlag )
00881         return false;
00882     else
00883         return true;
00884 }
00885 
00886 /*!This function checks the culling method criterion supplied to the object
00887   by the user.  If the user does not supply a culling method criterion,
00888   the default criterion is NONE, and in that case, no culling is performed.
00889   If the culling method criterion is satisfied, the interior vertices
00890   of the given patch are flagged as soft_fixed.  Otherwise, the soft_fixed
00891   flag is removed from each of the vertices in the patch (interior and
00892   boundary vertices).  Also, if the criterion was satisfied, then the
00893   function returns true.  Otherwise, the function returns false.
00894  */
00895 bool TerminationCriterion::cull_vertices( PatchData& pd, OFEvaluator& of_eval, MsqError& err )
00896 {
00897     // PRINT_INFO("CULLING_METHOD FLAG = %i",cullingMethodFlag);
00898 
00899     // cull_bool will be changed to true if the criterion is satisfied
00900     bool b, cull_bool = false;
00901     double prev_m, init_m;
00902     switch( cullingMethodFlag )
00903     {
00904             // if no culling is requested, always return false
00905         case NONE:
00906             return cull_bool;
00907             // if culling on quality improvement absolute
00908         case QUALITY_IMPROVEMENT_ABSOLUTE:
00909             // get objective function value
00910             b = of_eval.evaluate( pd, currentOFValue, err );
00911             if( MSQ_CHKERR( err ) ) return false;
00912             if( !b )
00913             {
00914                 MSQ_SETERR( err )( MsqError::INVALID_MESH );
00915                 return false;
00916             }
00917             // if the improvement was enough, cull
00918             if( currentOFValue <= cullingEps ) { cull_bool = true; }
00919             // PRINT_INFO("\ncurrentOFValue = %f, bool = %i\n",currentOFValue,cull_bool);
00920 
00921             break;
00922             // if culing on quality improvement relative
00923         case QUALITY_IMPROVEMENT_RELATIVE:
00924             // get objective function value
00925             b = of_eval.evaluate( pd, currentOFValue, err );
00926             if( MSQ_CHKERR( err ) ) return false;
00927             if( !b )
00928             {
00929                 MSQ_SETERR( err )( MsqError::INVALID_MESH );
00930                 return false;
00931             }
00932             // if the improvement was enough, cull
00933             if( ( currentOFValue - lowerOFBound ) <= ( cullingEps * ( initialOFValue - lowerOFBound ) ) )
00934             {
00935                 cull_bool = true;
00936             }
00937             break;
00938             // if culling on vertex movement absolute
00939         case VERTEX_MOVEMENT_ABSOLUTE:
00940         case VERTEX_MOVEMENT_ABS_EDGE_LENGTH:
00941             // if movement was enough, cull
00942             prev_m = pd.get_max_vertex_movement_squared( previousVerticesMemento, err );
00943             MSQ_ERRZERO( err );
00944             if( prev_m <= cullingEps * cullingEps ) { cull_bool = true; }
00945 
00946             break;
00947             // if culling on vertex movement relative
00948         case VERTEX_MOVEMENT_RELATIVE:
00949             // if movement was small enough, cull
00950             prev_m = pd.get_max_vertex_movement_squared( previousVerticesMemento, err );
00951             MSQ_ERRZERO( err );
00952             init_m = pd.get_max_vertex_movement_squared( initialVerticesMemento, err );
00953             MSQ_ERRZERO( err );
00954             if( prev_m <= ( cullingEps * cullingEps * init_m ) ) { cull_bool = true; }
00955             break;
00956         case UNTANGLED_MESH:
00957             if( !patchInvertedCount ) cull_bool = true;
00958             break;
00959         default:
00960             MSQ_SETERR( err )
00961             ( "Requested culling method not yet implemented.", MsqError::NOT_IMPLEMENTED );
00962             return false;
00963     };
00964     // Now actually have patch data cull vertices
00965     if( cull_bool )
00966     {
00967         pd.set_free_vertices_soft_fixed( err );
00968         MSQ_ERRZERO( err );
00969     }
00970     else
00971     {
00972         pd.set_all_vertices_soft_free( err );
00973         MSQ_ERRZERO( err );
00974     }
00975     return cull_bool;
00976 }
00977 
00978 /*!This function is activated when cullingGlobalPatch is true.  It supplies
00979   cull_vertices with a single vertex-based patch at a time.  If the patch
00980   satisfies the culling criterion, it's free vertices are then soft-fixed.
00981  */
00982 bool TerminationCriterion::cull_vertices_global( PatchData& global_patch, Mesh* mesh, MeshDomain* domain,
00983                                                  const Settings* settings, OFEvaluator& of_eval, MsqError& err )
00984 {
00985     if( !cullingGlobalPatch ) return false;
00986 
00987     // PRINT_INFO("CULLING_METHOD FLAG = %i",cullingMethodFlag);
00988 
00989     // cull_bool will be changed to true if the criterion is satisfied
00990     bool cull_bool = false;
00991 
00992     std::vector< Mesh::VertexHandle > mesh_vertices;
00993     // std::vector<Mesh::VertexHandle> patch_vertices;
00994     // std::vector<Mesh::ElementHandle> patch_elements;
00995     // std::vector<Mesh::VertexHandle> fixed_vertices;
00996     // std::vector<Mesh::VertexHandle> free_vertices;
00997 
00998     // FIXME, verify global_patch is a global patch... how, is this right?
00999     mesh->get_all_vertices( mesh_vertices, err );
01000     size_t mesh_num_nodes         = mesh_vertices.size();
01001     size_t global_patch_num_nodes = global_patch.num_nodes();
01002     if( 0 )
01003         std::cout << "tmp srk mesh_num_nodes= " << mesh_num_nodes
01004                   << " global_patch_num_nodes= " << global_patch_num_nodes << std::endl;
01005     if( mesh_num_nodes != global_patch_num_nodes )
01006     {
01007         std::cout << "tmp srk cull_vertices_global found non global patch" << std::endl;
01008         exit( 123 );
01009         return false;
01010     }
01011     PatchData patch;
01012     patch.set_mesh( (Mesh*)mesh );
01013     patch.set_domain( domain );
01014     patch.attach_settings( settings );
01015 
01016     // const MsqVertex* global_patch_vertex_array = global_patch.get_vertex_array( err );
01017     Mesh::VertexHandle* global_patch_vertex_handles = global_patch.get_vertex_handles_array();
01018 
01019     int num_culled = 0;
01020     for( unsigned iv = 0; iv < global_patch_num_nodes; iv++ )
01021     {
01022         // form a patch for this vertex; if it is culled, set it to be soft fixed
01023         Mesh::VertexHandle vert = global_patch_vertex_handles[iv];
01024         std::vector< Mesh::ElementHandle > elements;
01025         std::vector< size_t > offsets;
01026         mesh->vertices_get_attached_elements( &vert, 1, elements, offsets, err );
01027 
01028         std::set< Mesh::VertexHandle > patch_free_vertices_set;
01029 
01030         for( unsigned ie = 0; ie < elements.size(); ie++ )
01031         {
01032             std::vector< Mesh::VertexHandle > vert_handles;
01033             std::vector< size_t > v_offsets;
01034             mesh->elements_get_attached_vertices( &elements[ie], 1, vert_handles, v_offsets, err );
01035             for( unsigned jv = 0; jv < vert_handles.size(); jv++ )
01036             {
01037                 unsigned char bt;
01038                 mesh->vertex_get_byte( vert_handles[jv], &bt, err );
01039                 MsqVertex v;
01040                 v.set_flags( bt );
01041                 if( v.is_free_vertex() ) patch_free_vertices_set.insert( vert_handles[jv] );
01042             }
01043         }
01044 
01045         std::vector< Mesh::VertexHandle > patch_free_vertices_vector( patch_free_vertices_set.begin(),
01046                                                                       patch_free_vertices_set.end() );
01047         // std::vector<unsigned char> byte_vector(patch_vertices_vector.size());
01048         // mesh->vertices_get_byte(&vert_handles[0], &byte_vector[0], vert_handles.size(), err);
01049 
01050         patch.set_mesh_entities( elements, patch_free_vertices_vector, err );
01051         if( cull_vertices( patch, of_eval, err ) )
01052         {
01053             // std::cout << "tmp srk cull_vertices_global found culled patch" << std::endl;
01054             Mesh::VertexHandle* patch_vertex_handles = patch.get_vertex_handles_array();
01055             const MsqVertex* patch_vertex_array      = patch.get_vertex_array( err );
01056             for( unsigned jv = 0; jv < patch.num_nodes(); jv++ )
01057             {
01058                 if( patch_vertex_handles[jv] == global_patch_vertex_handles[iv] )
01059                 {
01060                     if( patch_vertex_array[jv].is_flag_set( MsqVertex::MSQ_CULLED ) )
01061                     {
01062                         global_patch.set_vertex_culled( iv );
01063                         ++num_culled;
01064                         cull_bool = true;
01065                         // std::cout << "tmp srk cull_vertices_global found culled vertex" <<
01066                         // std::endl;
01067                     }
01068                 }
01069             }
01070         }
01071     }
01072     if( 0 )
01073         std::cout << "tmp srk cull_vertices_global found " << num_culled << " culled vertices out of "
01074                   << global_patch_num_nodes << std::endl;
01075 
01076     return cull_bool;
01077 }
01078 
01079 size_t TerminationCriterion::count_inverted( PatchData& pd, MsqError& err )
01080 {
01081     size_t num_elem = pd.num_elements();
01082     size_t count    = 0;
01083     int inverted, samples;
01084     for( size_t i = 0; i < num_elem; i++ )
01085     {
01086         pd.element_by_index( i ).check_element_orientation( pd, inverted, samples, err );
01087         if( inverted ) ++count;
01088     }
01089     return count;
01090 }
01091 
01092 /*!
01093   Currently this only deletes the memento of the vertex positions and the
01094   mGrad vector if neccessary.
01095   When culling, we remove the soft fixed flags from all of the vertices.
01096  */
01097 void TerminationCriterion::cleanup( Mesh*, MeshDomain*, MsqError& )
01098 {
01099     delete previousVerticesMemento;
01100     delete initialVerticesMemento;
01101     previousVerticesMemento = 0;
01102     initialVerticesMemento  = 0;
01103 }
01104 
01105 void TerminationCriterion::write_timestep( PatchData& pd, const Vector3D* gradient, MsqError& err )
01106 {
01107     std::ostringstream str;
01108     if( timeStepFileType == VTK )
01109     {
01110         str << timeStepFileName << '_' << iterationCounter << ".vtk";
01111         MeshWriter::write_vtk( pd, str.str().c_str(), err, gradient );
01112     }
01113     else if( timeStepFileType == GNUPLOT )
01114     {
01115         str << timeStepFileName << '.' << iterationCounter;
01116         MeshWriter::write_gnuplot( pd.get_mesh(), str.str().c_str(), err );
01117     }
01118 }
01119 
01120 void TerminationCriterion::write_iterations( const char* filename, MsqError& err )
01121 {
01122     if( filename )
01123     {
01124         plotFile.open( filename, std::ios::trunc );
01125         if( !plotFile ) MSQ_SETERR( err )
01126         ( MsqError::FILE_ACCESS, "Failed to open plot data file: '%s'", filename );
01127     }
01128     else
01129     {
01130         plotFile.close();
01131     }
01132 }
01133 
01134 void TerminationCriterion::initialize_queue( MeshDomainAssoc* mesh_and_domain, const Settings*, MsqError& err )
01135 {
01136     if( VERTEX_MOVEMENT_ABS_EDGE_LENGTH & ( terminationCriterionFlag | cullingMethodFlag ) )
01137     {
01138         Mesh* mesh = mesh_and_domain->get_mesh();
01139         MeshUtil tool( mesh );
01140         SimpleStats stats;
01141         tool.edge_length_distribution( stats, err );MSQ_ERRRTN( err );
01142         double limit = vertexMovementAvgBeta * ( stats.average() - stats.standard_deviation() );
01143         // we actually calculate the square of the length
01144         vertexMovementAbsoluteAvgEdge = limit * limit;
01145         if( VERTEX_MOVEMENT_ABS_EDGE_LENGTH & cullingMethodFlag ) cullingEps = limit;
01146     }
01147 }
01148 
01149 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines