MOAB: Mesh Oriented datABase  (version 5.2.1)
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             { ++vertexMovementExceedsBound; }
00664         }
00665     }
00666 
00667     if( ( terminationCriterionFlag | cullingMethodFlag ) & UNTANGLED_MESH )
00668     {
00669         size_t new_count = count_inverted( pd, err );
00670         // be careful here because size_t is unsigned
00671         globalInvertedCount += new_count;
00672         globalInvertedCount -= patchInvertedCount;
00673         patchInvertedCount = new_count;
00674         // if (innerOuterType==TYPE_OUTER)
00675         //  MSQ_DBGOUT_P0_ONLY(debugLevel) << par_string() << "  o Num Patch Inverted: " << " " <<
00676         //  patchInvertedCount << " globalInvertedCount= " << globalInvertedCount << std::endl;
00677 
00678         MSQ_ERRRTN( err );
00679     }
00680 }
00681 
00682 /*!  This function evaluates the needed information and then evaluates
00683   the termination criteria.  If any of the selected criteria are satisfied,
00684   the function returns true.  Otherwise, the function returns false.
00685  */
00686 bool TerminationCriterion::terminate()
00687 {
00688     bool return_flag = false;
00689     // std::cout<<"\nInside terminate(pd,of,err):  flag = "<<terminationCriterionFlag << std::endl;
00690     int type = 0;
00691 
00692     // First check for an interrupt signal
00693     if( MsqInterrupt::interrupt() )
00694     {
00695         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- INTERRUPTED"
00696                                          << " " << std::endl;
00697         return true;
00698     }
00699 
00700     // if terminating on numbering of inner iterations
00701     if( NUMBER_OF_ITERATES & terminationCriterionFlag && iterationCounter >= iterationBound )
00702     {
00703         return_flag = true;
00704         type        = 1;
00705         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- Reached "
00706                                          << " " << iterationBound << " iterations." << std::endl;
00707     }
00708 
00709     if( CPU_TIME & terminationCriterionFlag && mTimer.since_birth() >= timeBound )
00710     {
00711         return_flag = true;
00712         type        = 2;
00713         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- Exceeded CPU time. "
00714                                          << " " << std::endl;
00715     }
00716 
00717     if( MOVEMENT_FLAGS & terminationCriterionFlag && maxSquaredMovement >= 0.0 )
00718     {
00719         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o Info -- Maximum vertex movement: "
00720                                          << " " << RPM( sqrt( maxSquaredMovement ) ) << std::endl;
00721 
00722         if( VERTEX_MOVEMENT_ABSOLUTE & terminationCriterionFlag && maxSquaredMovement <= vertexMovementAbsoluteEps )
00723         {
00724             MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- VERTEX_MOVEMENT_ABSOLUTE: "
00725                                              << " " << RPM( sqrt( maxSquaredMovement ) ) << std::endl;
00726             return_flag = true;
00727             type        = 3;
00728         }
00729 
00730         if( VERTEX_MOVEMENT_RELATIVE & terminationCriterionFlag &&
00731             maxSquaredMovement <= vertexMovementRelativeEps * maxSquaredInitialMovement )
00732         {
00733             MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- VERTEX_MOVEMENT_RELATIVE: "
00734                                              << " " << RPM( sqrt( maxSquaredMovement ) ) << std::endl;
00735             return_flag = true;
00736             type        = 4;
00737         }
00738 
00739         if( VERTEX_MOVEMENT_ABS_EDGE_LENGTH & terminationCriterionFlag )
00740         {
00741             assert( vertexMovementAbsoluteAvgEdge > -1e-12 );  // make sure value actually got calculated
00742             if( maxSquaredMovement <= vertexMovementAbsoluteAvgEdge )
00743             {
00744                 MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- VERTEX_MOVEMENT_ABS_EDGE_LENGTH: "
00745                                                  << " " << RPM( sqrt( maxSquaredMovement ) ) << std::endl;
00746                 return_flag = true;
00747                 type        = 5;
00748             }
00749         }
00750     }
00751 
00752     if( GRADIENT_L2_NORM_ABSOLUTE & terminationCriterionFlag &&
00753         currentGradL2NormSquared <= gradL2NormAbsoluteEpsSquared )
00754     {
00755         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- GRADIENT_L2_NORM_ABSOLUTE: "
00756                                          << " " << RPM( currentGradL2NormSquared ) << std::endl;
00757         return_flag = true;
00758         type        = 6;
00759     }
00760 
00761     if( GRADIENT_INF_NORM_ABSOLUTE & terminationCriterionFlag && currentGradInfNorm <= gradInfNormAbsoluteEps )
00762     {
00763         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- GRADIENT_INF_NORM_ABSOLUTE: "
00764                                          << " " << RPM( currentGradInfNorm ) << std::endl;
00765         return_flag = true;
00766         type        = 7;
00767     }
00768 
00769     if( GRADIENT_L2_NORM_RELATIVE & terminationCriterionFlag &&
00770         currentGradL2NormSquared <= ( gradL2NormRelativeEpsSquared * initialGradL2NormSquared ) )
00771     {
00772         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- GRADIENT_L2_NORM_RELATIVE: "
00773                                          << " " << RPM( currentGradL2NormSquared ) << std::endl;
00774         return_flag = true;
00775         type        = 8;
00776     }
00777 
00778     if( GRADIENT_INF_NORM_RELATIVE & terminationCriterionFlag &&
00779         currentGradInfNorm <= ( gradInfNormRelativeEps * initialGradInfNorm ) )
00780     {
00781         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- GRADIENT_INF_NORM_RELATIVE: "
00782                                          << " " << RPM( currentGradInfNorm ) << std::endl;
00783         return_flag = true;
00784         type        = 9;
00785     }
00786     // Quality Improvement and Successive Improvements are below.
00787     // The relative forms are only valid after the first iteration.
00788     if( ( QUALITY_IMPROVEMENT_ABSOLUTE & terminationCriterionFlag ) && currentOFValue <= qualityImprovementAbsoluteEps )
00789     {
00790         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- QUALITY_IMPROVEMENT_ABSOLUTE: "
00791                                          << " " << RPM( currentOFValue ) << std::endl;
00792         return_flag = true;
00793         type        = 10;
00794     }
00795 
00796     // only valid if an iteration has occurred, see above.
00797     if( iterationCounter > 0 )
00798     {
00799         if( SUCCESSIVE_IMPROVEMENTS_ABSOLUTE & terminationCriterionFlag &&
00800             ( previousOFValue - currentOFValue ) <= successiveImprovementsAbsoluteEps )
00801         {
00802             MSQ_DBGOUT_P0_ONLY( debugLevel )
00803                 << par_string() << "  o TermCrit -- SUCCESSIVE_IMPROVEMENTS_ABSOLUTE: previousOFValue= "
00804                 << " " << previousOFValue << " currentOFValue= " << RPM( currentOFValue )
00805                 << " successiveImprovementsAbsoluteEps= " << successiveImprovementsAbsoluteEps << std::endl;
00806             return_flag = true;
00807             type        = 11;
00808         }
00809         if( QUALITY_IMPROVEMENT_RELATIVE & terminationCriterionFlag &&
00810             ( currentOFValue - lowerOFBound ) <= qualityImprovementRelativeEps * ( initialOFValue - lowerOFBound ) )
00811         {
00812             MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- QUALITY_IMPROVEMENT_RELATIVE: "
00813                                              << " " << std::endl;
00814             return_flag = true;
00815             type        = 12;
00816         }
00817         if( SUCCESSIVE_IMPROVEMENTS_RELATIVE & terminationCriterionFlag &&
00818             ( previousOFValue - currentOFValue ) <=
00819                 successiveImprovementsRelativeEps * ( initialOFValue - currentOFValue ) )
00820         {
00821             MSQ_DBGOUT_P0_ONLY( debugLevel )
00822                 << par_string() << "  o TermCrit -- SUCCESSIVE_IMPROVEMENTS_RELATIVE: previousOFValue= "
00823                 << " " << previousOFValue << " currentOFValue= " << RPM( currentOFValue )
00824                 << " successiveImprovementsRelativeEps= " << successiveImprovementsRelativeEps << std::endl;
00825             return_flag = true;
00826             type        = 13;
00827         }
00828     }
00829 
00830     if( BOUNDED_VERTEX_MOVEMENT & terminationCriterionFlag && vertexMovementExceedsBound )
00831     {
00832         return_flag = true;
00833         type        = 14;
00834         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- "
00835                                          << " " << vertexMovementExceedsBound << " vertices out of bounds."
00836                                          << std::endl;
00837     }
00838 
00839     if( UNTANGLED_MESH & terminationCriterionFlag )
00840     {
00841         if( innerOuterType == TYPE_OUTER )
00842         {
00843             // MSQ_DBGOUT_P0_ONLY(debugLevel)
00844             MSQ_DBGOUT( debugLevel ) << par_string() << "  o Num Inverted: "
00845                                      << " " << globalInvertedCount << std::endl;
00846         }
00847         if( !globalInvertedCount )
00848         {
00849             MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- UNTANGLED_MESH: "
00850                                              << " " << std::endl;
00851             return_flag = true;
00852             type        = 15;
00853         }
00854     }
00855 
00856     // clear this value at the end of each iteration
00857     vertexMovementExceedsBound = 0;
00858     maxSquaredMovement         = -1.0;
00859 
00860     if( timeStepFileType == GNUPLOT && return_flag )
00861     {
00862         MsqError err;
00863         MeshWriter::write_gnuplot_overlay( iterationCounter, timeStepFileName.c_str(), err );
00864     }
00865 
00866     if( 0 && return_flag && MSQ_DBG( 2 ) )
00867         std::cout << "P[" << get_parallel_rank() << "] tmp TerminationCriterion::terminate: " << moniker
00868                   << " return_flag= " << return_flag << " type= " << type
00869                   << " terminationCriterionFlag= " << terminationCriterionFlag << " debugLevel= " << debugLevel
00870                   << std::endl;
00871 
00872     // if none of the criteria were satisfied
00873     return return_flag;
00874 }
00875 
00876 bool TerminationCriterion::criterion_is_set()
00877 {
00878     if( !terminationCriterionFlag )
00879         return false;
00880     else
00881         return true;
00882 }
00883 
00884 /*!This function checks the culling method criterion supplied to the object
00885   by the user.  If the user does not supply a culling method criterion,
00886   the default criterion is NONE, and in that case, no culling is performed.
00887   If the culling method criterion is satisfied, the interior vertices
00888   of the given patch are flagged as soft_fixed.  Otherwise, the soft_fixed
00889   flag is removed from each of the vertices in the patch (interior and
00890   boundary vertices).  Also, if the criterion was satisfied, then the
00891   function returns true.  Otherwise, the function returns false.
00892  */
00893 bool TerminationCriterion::cull_vertices( PatchData& pd, OFEvaluator& of_eval, MsqError& err )
00894 {
00895     // PRINT_INFO("CULLING_METHOD FLAG = %i",cullingMethodFlag);
00896 
00897     // cull_bool will be changed to true if the criterion is satisfied
00898     bool b, cull_bool = false;
00899     double prev_m, init_m;
00900     switch( cullingMethodFlag )
00901     {
00902             // if no culling is requested, always return false
00903         case NONE:
00904             return cull_bool;
00905             // if culling on quality improvement absolute
00906         case QUALITY_IMPROVEMENT_ABSOLUTE:
00907             // get objective function value
00908             b = of_eval.evaluate( pd, currentOFValue, err );
00909             if( MSQ_CHKERR( err ) ) return false;
00910             if( !b )
00911             {
00912                 MSQ_SETERR( err )( MsqError::INVALID_MESH );
00913                 return false;
00914             }
00915             // if the improvement was enough, cull
00916             if( currentOFValue <= cullingEps ) { cull_bool = true; }
00917             // PRINT_INFO("\ncurrentOFValue = %f, bool = %i\n",currentOFValue,cull_bool);
00918 
00919             break;
00920             // if culing on quality improvement relative
00921         case QUALITY_IMPROVEMENT_RELATIVE:
00922             // get objective function value
00923             b = of_eval.evaluate( pd, currentOFValue, err );
00924             if( MSQ_CHKERR( err ) ) return false;
00925             if( !b )
00926             {
00927                 MSQ_SETERR( err )( MsqError::INVALID_MESH );
00928                 return false;
00929             }
00930             // if the improvement was enough, cull
00931             if( ( currentOFValue - lowerOFBound ) <= ( cullingEps * ( initialOFValue - lowerOFBound ) ) )
00932             { cull_bool = true; }
00933             break;
00934             // if culling on vertex movement absolute
00935         case VERTEX_MOVEMENT_ABSOLUTE:
00936         case VERTEX_MOVEMENT_ABS_EDGE_LENGTH:
00937             // if movement was enough, cull
00938             prev_m = pd.get_max_vertex_movement_squared( previousVerticesMemento, err );
00939             MSQ_ERRZERO( err );
00940             if( prev_m <= cullingEps * cullingEps ) { cull_bool = true; }
00941 
00942             break;
00943             // if culling on vertex movement relative
00944         case VERTEX_MOVEMENT_RELATIVE:
00945             // if movement was small enough, cull
00946             prev_m = pd.get_max_vertex_movement_squared( previousVerticesMemento, err );
00947             MSQ_ERRZERO( err );
00948             init_m = pd.get_max_vertex_movement_squared( initialVerticesMemento, err );
00949             MSQ_ERRZERO( err );
00950             if( prev_m <= ( cullingEps * cullingEps * init_m ) ) { cull_bool = true; }
00951             break;
00952         case UNTANGLED_MESH:
00953             if( !patchInvertedCount ) cull_bool = true;
00954             break;
00955         default:
00956             MSQ_SETERR( err )
00957             ( "Requested culling method not yet implemented.", MsqError::NOT_IMPLEMENTED );
00958             return false;
00959     };
00960     // Now actually have patch data cull vertices
00961     if( cull_bool )
00962     {
00963         pd.set_free_vertices_soft_fixed( err );
00964         MSQ_ERRZERO( err );
00965     }
00966     else
00967     {
00968         pd.set_all_vertices_soft_free( err );
00969         MSQ_ERRZERO( err );
00970     }
00971     return cull_bool;
00972 }
00973 
00974 /*!This function is activated when cullingGlobalPatch is true.  It supplies
00975   cull_vertices with a single vertex-based patch at a time.  If the patch
00976   satisfies the culling criterion, it's free vertices are then soft-fixed.
00977  */
00978 bool TerminationCriterion::cull_vertices_global( PatchData& global_patch, Mesh* mesh, MeshDomain* domain,
00979                                                  const Settings* settings, OFEvaluator& of_eval, MsqError& err )
00980 {
00981     if( !cullingGlobalPatch ) return false;
00982 
00983     // PRINT_INFO("CULLING_METHOD FLAG = %i",cullingMethodFlag);
00984 
00985     // cull_bool will be changed to true if the criterion is satisfied
00986     bool cull_bool = false;
00987 
00988     std::vector< Mesh::VertexHandle > mesh_vertices;
00989     // std::vector<Mesh::VertexHandle> patch_vertices;
00990     // std::vector<Mesh::ElementHandle> patch_elements;
00991     // std::vector<Mesh::VertexHandle> fixed_vertices;
00992     // std::vector<Mesh::VertexHandle> free_vertices;
00993 
00994     // FIXME, verify global_patch is a global patch... how, is this right?
00995     mesh->get_all_vertices( mesh_vertices, err );
00996     size_t mesh_num_nodes         = mesh_vertices.size();
00997     size_t global_patch_num_nodes = global_patch.num_nodes();
00998     if( 0 )
00999         std::cout << "tmp srk mesh_num_nodes= " << mesh_num_nodes
01000                   << " global_patch_num_nodes= " << global_patch_num_nodes << std::endl;
01001     if( mesh_num_nodes != global_patch_num_nodes )
01002     {
01003         std::cout << "tmp srk cull_vertices_global found non global patch" << std::endl;
01004         exit( 123 );
01005         return false;
01006     }
01007     PatchData patch;
01008     patch.set_mesh( (Mesh*)mesh );
01009     patch.set_domain( domain );
01010     patch.attach_settings( settings );
01011 
01012     // const MsqVertex* global_patch_vertex_array = global_patch.get_vertex_array( err );
01013     Mesh::VertexHandle* global_patch_vertex_handles = global_patch.get_vertex_handles_array();
01014 
01015     int num_culled = 0;
01016     for( unsigned iv = 0; iv < global_patch_num_nodes; iv++ )
01017     {
01018         // form a patch for this vertex; if it is culled, set it to be soft fixed
01019         Mesh::VertexHandle vert = global_patch_vertex_handles[iv];
01020         std::vector< Mesh::ElementHandle > elements;
01021         std::vector< size_t > offsets;
01022         mesh->vertices_get_attached_elements( &vert, 1, elements, offsets, err );
01023 
01024         std::set< Mesh::VertexHandle > patch_free_vertices_set;
01025 
01026         for( unsigned ie = 0; ie < elements.size(); ie++ )
01027         {
01028             std::vector< Mesh::VertexHandle > vert_handles;
01029             std::vector< size_t > v_offsets;
01030             mesh->elements_get_attached_vertices( &elements[ie], 1, vert_handles, v_offsets, err );
01031             for( unsigned jv = 0; jv < vert_handles.size(); jv++ )
01032             {
01033                 unsigned char bt;
01034                 mesh->vertex_get_byte( vert_handles[jv], &bt, err );
01035                 MsqVertex v;
01036                 v.set_flags( bt );
01037                 if( v.is_free_vertex() ) patch_free_vertices_set.insert( vert_handles[jv] );
01038             }
01039         }
01040 
01041         std::vector< Mesh::VertexHandle > patch_free_vertices_vector( patch_free_vertices_set.begin(),
01042                                                                       patch_free_vertices_set.end() );
01043         // std::vector<unsigned char> byte_vector(patch_vertices_vector.size());
01044         // mesh->vertices_get_byte(&vert_handles[0], &byte_vector[0], vert_handles.size(), err);
01045 
01046         patch.set_mesh_entities( elements, patch_free_vertices_vector, err );
01047         if( cull_vertices( patch, of_eval, err ) )
01048         {
01049             // std::cout << "tmp srk cull_vertices_global found culled patch" << std::endl;
01050             Mesh::VertexHandle* patch_vertex_handles = patch.get_vertex_handles_array();
01051             const MsqVertex* patch_vertex_array      = patch.get_vertex_array( err );
01052             for( unsigned jv = 0; jv < patch.num_nodes(); jv++ )
01053             {
01054                 if( patch_vertex_handles[jv] == global_patch_vertex_handles[iv] )
01055                 {
01056                     if( patch_vertex_array[jv].is_flag_set( MsqVertex::MSQ_CULLED ) )
01057                     {
01058                         global_patch.set_vertex_culled( iv );
01059                         ++num_culled;
01060                         cull_bool = true;
01061                         // std::cout << "tmp srk cull_vertices_global found culled vertex" <<
01062                         // std::endl;
01063                     }
01064                 }
01065             }
01066         }
01067     }
01068     if( 0 )
01069         std::cout << "tmp srk cull_vertices_global found " << num_culled << " culled vertices out of "
01070                   << global_patch_num_nodes << std::endl;
01071 
01072     return cull_bool;
01073 }
01074 
01075 size_t TerminationCriterion::count_inverted( PatchData& pd, MsqError& err )
01076 {
01077     size_t num_elem = pd.num_elements();
01078     size_t count    = 0;
01079     int inverted, samples;
01080     for( size_t i = 0; i < num_elem; i++ )
01081     {
01082         pd.element_by_index( i ).check_element_orientation( pd, inverted, samples, err );
01083         if( inverted ) ++count;
01084     }
01085     return count;
01086 }
01087 
01088 /*!
01089   Currently this only deletes the memento of the vertex positions and the
01090   mGrad vector if neccessary.
01091   When culling, we remove the soft fixed flags from all of the vertices.
01092  */
01093 void TerminationCriterion::cleanup( Mesh*, MeshDomain*, MsqError& )
01094 {
01095     delete previousVerticesMemento;
01096     delete initialVerticesMemento;
01097     previousVerticesMemento = 0;
01098     initialVerticesMemento  = 0;
01099 }
01100 
01101 void TerminationCriterion::write_timestep( PatchData& pd, const Vector3D* gradient, MsqError& err )
01102 {
01103     std::ostringstream str;
01104     if( timeStepFileType == VTK )
01105     {
01106         str << timeStepFileName << '_' << iterationCounter << ".vtk";
01107         MeshWriter::write_vtk( pd, str.str().c_str(), err, gradient );
01108     }
01109     else if( timeStepFileType == GNUPLOT )
01110     {
01111         str << timeStepFileName << '.' << iterationCounter;
01112         MeshWriter::write_gnuplot( pd.get_mesh(), str.str().c_str(), err );
01113     }
01114 }
01115 
01116 void TerminationCriterion::write_iterations( const char* filename, MsqError& err )
01117 {
01118     if( filename )
01119     {
01120         plotFile.open( filename, std::ios::trunc );
01121         if( !plotFile ) MSQ_SETERR( err )
01122         ( MsqError::FILE_ACCESS, "Failed to open plot data file: '%s'", filename );
01123     }
01124     else
01125     {
01126         plotFile.close();
01127     }
01128 }
01129 
01130 void TerminationCriterion::initialize_queue( MeshDomainAssoc* mesh_and_domain, const Settings*, MsqError& err )
01131 {
01132     if( VERTEX_MOVEMENT_ABS_EDGE_LENGTH & ( terminationCriterionFlag | cullingMethodFlag ) )
01133     {
01134         Mesh* mesh = mesh_and_domain->get_mesh();
01135         MeshUtil tool( mesh );
01136         SimpleStats stats;
01137         tool.edge_length_distribution( stats, err );MSQ_ERRRTN( err );
01138         double limit = vertexMovementAvgBeta * ( stats.average() - stats.standard_deviation() );
01139         // we actually calculate the square of the length
01140         vertexMovementAbsoluteAvgEdge = limit * limit;
01141         if( VERTEX_MOVEMENT_ABS_EDGE_LENGTH & cullingMethodFlag ) cullingEps = limit;
01142     }
01143 }
01144 
01145 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines